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.math3.analysis.interpolation;
018    
019    import org.apache.commons.math3.exception.DimensionMismatchException;
020    import org.apache.commons.math3.exception.NoDataException;
021    import org.apache.commons.math3.exception.NotPositiveException;
022    import org.apache.commons.math3.util.MathArrays;
023    import org.apache.commons.math3.util.Precision;
024    import org.apache.commons.math3.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer;
025    import org.apache.commons.math3.fitting.PolynomialFitter;
026    import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
027    import org.apache.commons.math3.optim.SimpleVectorValueChecker;
028    
029    /**
030     * Generates a bicubic interpolation function.
031     * Prior to generating the interpolating function, the input is smoothed using
032     * polynomial fitting.
033     *
034     * @version $Id: SmoothingPolynomialBicubicSplineInterpolator.java 1423558 2012-12-18 18:07:45Z erans $
035     * @since 2.2
036     */
037    public class SmoothingPolynomialBicubicSplineInterpolator
038        extends BicubicSplineInterpolator {
039        /** Fitter for x. */
040        private final PolynomialFitter xFitter;
041        /** Degree of the fitting polynomial. */
042        private final int xDegree;
043        /** Fitter for y. */
044        private final PolynomialFitter yFitter;
045        /** Degree of the fitting polynomial. */
046        private final int yDegree;
047    
048        /**
049         * Default constructor. The degree of the fitting polynomials is set to 3.
050         */
051        public SmoothingPolynomialBicubicSplineInterpolator() {
052            this(3);
053        }
054    
055        /**
056         * @param degree Degree of the polynomial fitting functions.
057         */
058        public SmoothingPolynomialBicubicSplineInterpolator(int degree) {
059            this(degree, degree);
060        }
061    
062        /**
063         * @param xDegree Degree of the polynomial fitting functions along the
064         * x-dimension.
065         * @param yDegree Degree of the polynomial fitting functions along the
066         * y-dimension.
067         */
068        public SmoothingPolynomialBicubicSplineInterpolator(int xDegree,
069                                                            int yDegree) {
070            if (xDegree < 0) {
071                throw new NotPositiveException(xDegree);
072            }
073            if (yDegree < 0) {
074                throw new NotPositiveException(yDegree);
075            }
076            this.xDegree = xDegree;
077            this.yDegree = yDegree;
078    
079            final double safeFactor = 1e2;
080            final SimpleVectorValueChecker checker
081                = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON,
082                                               safeFactor * Precision.SAFE_MIN);
083            xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
084            yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
085        }
086    
087        /**
088         * {@inheritDoc}
089         */
090        @Override
091        public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
092                                                              final double[] yval,
093                                                              final double[][] fval)
094            throws NoDataException,
095                   DimensionMismatchException {
096            if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
097                throw new NoDataException();
098            }
099            if (xval.length != fval.length) {
100                throw new DimensionMismatchException(xval.length, fval.length);
101            }
102    
103            final int xLen = xval.length;
104            final int yLen = yval.length;
105    
106            for (int i = 0; i < xLen; i++) {
107                if (fval[i].length != yLen) {
108                    throw new DimensionMismatchException(fval[i].length, yLen);
109                }
110            }
111    
112            MathArrays.checkOrder(xval);
113            MathArrays.checkOrder(yval);
114    
115            // For each line y[j] (0 <= j < yLen), construct a polynomial, with
116            // respect to variable x, fitting array fval[][j]
117            final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
118            for (int j = 0; j < yLen; j++) {
119                xFitter.clearObservations();
120                for (int i = 0; i < xLen; i++) {
121                    xFitter.addObservedPoint(1, xval[i], fval[i][j]);
122                }
123    
124                // Initial guess for the fit is zero for each coefficients (of which
125                // there are "xDegree" + 1).
126                yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1]));
127            }
128    
129            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
130            // values fval_1
131            final double[][] fval_1 = new double[xLen][yLen];
132            for (int j = 0; j < yLen; j++) {
133                final PolynomialFunction f = yPolyX[j];
134                for (int i = 0; i < xLen; i++) {
135                    fval_1[i][j] = f.value(xval[i]);
136                }
137            }
138    
139            // For each line x[i] (0 <= i < xLen), construct a polynomial, with
140            // respect to variable y, fitting array fval_1[i][]
141            final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
142            for (int i = 0; i < xLen; i++) {
143                yFitter.clearObservations();
144                for (int j = 0; j < yLen; j++) {
145                    yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
146                }
147    
148                // Initial guess for the fit is zero for each coefficients (of which
149                // there are "yDegree" + 1).
150                xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1]));
151            }
152    
153            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
154            // values fval_2
155            final double[][] fval_2 = new double[xLen][yLen];
156            for (int i = 0; i < xLen; i++) {
157                final PolynomialFunction f = xPolyY[i];
158                for (int j = 0; j < yLen; j++) {
159                    fval_2[i][j] = f.value(yval[j]);
160                }
161            }
162    
163            return super.interpolate(xval, yval, fval_2);
164        }
165    }