001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.analysis.interpolation;
018    
019    import org.apache.commons.math.exception.DimensionMismatchException;
020    import org.apache.commons.math.exception.NoDataException;
021    import org.apache.commons.math.MathException;
022    import org.apache.commons.math.util.MathUtils;
023    import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
024    import org.apache.commons.math.optimization.fitting.PolynomialFitter;
025    import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
026    
027    /**
028     * Generates a bicubic interpolation function.
029     * Prior to generating the interpolating function, the input is smoothed using
030     * polynomial fitting.
031     *
032     * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
033     * @since 2.2
034     */
035    public class SmoothingPolynomialBicubicSplineInterpolator
036        extends BicubicSplineInterpolator {
037    
038        /** Fitter for x. */
039        private final PolynomialFitter xFitter;
040    
041        /** Fitter for y. */
042        private final PolynomialFitter yFitter;
043    
044        /**
045         * Default constructor. The degree of the fitting polynomials is set to 3.
046         */
047        public SmoothingPolynomialBicubicSplineInterpolator() {
048            this(3);
049        }
050    
051        /**
052         * @param degree Degree of the polynomial fitting functions.
053         */
054        public SmoothingPolynomialBicubicSplineInterpolator(int degree) {
055            this(degree, degree);
056        }
057    
058        /**
059         * @param xDegree Degree of the polynomial fitting functions along the
060         * x-dimension.
061         * @param yDegree Degree of the polynomial fitting functions along the
062         * y-dimension.
063         */
064        public SmoothingPolynomialBicubicSplineInterpolator(int xDegree,
065                                                            int yDegree) {
066            xFitter = new PolynomialFitter(xDegree, new GaussNewtonOptimizer(false));
067            yFitter = new PolynomialFitter(yDegree, new GaussNewtonOptimizer(false));
068        }
069    
070        /**
071         * {@inheritDoc}
072         */
073        @Override
074        public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
075                                                              final double[] yval,
076                                                              final double[][] fval)
077            throws MathException {
078            if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
079                throw new NoDataException();
080            }
081            if (xval.length != fval.length) {
082                throw new DimensionMismatchException(xval.length, fval.length);
083            }
084    
085            final int xLen = xval.length;
086            final int yLen = yval.length;
087    
088            for (int i = 0; i < xLen; i++) {
089                if (fval[i].length != yLen) {
090                    throw new DimensionMismatchException(fval[i].length, yLen);
091                }
092            }
093    
094            MathUtils.checkOrder(xval);
095            MathUtils.checkOrder(yval);
096    
097            // For each line y[j] (0 <= j < yLen), construct a polynomial, with
098            // respect to variable x, fitting array fval[][j]
099            final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
100            for (int j = 0; j < yLen; j++) {
101                xFitter.clearObservations();
102                for (int i = 0; i < xLen; i++) {
103                    xFitter.addObservedPoint(1, xval[i], fval[i][j]);
104                }
105    
106                yPolyX[j] = xFitter.fit();
107            }
108    
109            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
110            // values fval_1
111            final double[][] fval_1 = new double[xLen][yLen];
112            for (int j = 0; j < yLen; j++) {
113                final PolynomialFunction f = yPolyX[j];
114                for (int i = 0; i < xLen; i++) {
115                    fval_1[i][j] = f.value(xval[i]);
116                }
117            }
118    
119            // For each line x[i] (0 <= i < xLen), construct a polynomial, with
120            // respect to variable y, fitting array fval_1[i][]
121            final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
122            for (int i = 0; i < xLen; i++) {
123                yFitter.clearObservations();
124                for (int j = 0; j < yLen; j++) {
125                    yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
126                }
127    
128                xPolyY[i] = yFitter.fit();
129            }
130    
131            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
132            // values fval_2
133            final double[][] fval_2 = new double[xLen][yLen];
134            for (int i = 0; i < xLen; i++) {
135                final PolynomialFunction f = xPolyY[i];
136                for (int j = 0; j < yLen; j++) {
137                    fval_2[i][j] = f.value(yval[j]);
138                }
139            }
140    
141            return super.interpolate(xval, yval, fval_2);
142        }
143    }