    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    }