    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    }