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