    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.DimensionMismatchException;
    020    import org.apache.commons.math.MathRuntimeException;
    021    import org.apache.commons.math.MathException;
    022    import org.apache.commons.math.util.MathUtils;
    023    import org.apache.commons.math.util.MathUtils.OrderDirection;
    024    import org.apache.commons.math.analysis.BivariateRealFunction;
    025    import org.apache.commons.math.analysis.UnivariateRealFunction;
    026    import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
    027    import org.apache.commons.math.exception.util.LocalizedFormats;
    028
    029    /**
    030     * Generates a bicubic interpolation function.
    031     * Before interpolating, smoothing of the input data is performed using
    032     * splines.
    033     * See <b>Handbook on splines for the user</b>, ISBN 084939404X,
    034     * chapter 2.
    035     *
    036     * @version $Revision: 1059400 $ $Date: 2011-01-15 20:35:27 +0100 (sam. 15 janv. 2011) $
    037     * @since 2.1
    038     * @deprecated This class does not perform smoothing; the name is thus misleading.
    039     * Please use {@link org.apache.commons.math.analysis.interpolation.BicubicSplineInterpolator}
    040     * instead. If smoothing is desired, a tentative implementation is provided in class
    041     * {@link org.apache.commons.math.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator}.
    042     * This class will be removed in math 3.0.
    043     */
    044    @Deprecated
    045    public class SmoothingBicubicSplineInterpolator
    046        implements BivariateRealGridInterpolator {
    047        /**
    048         * {@inheritDoc}
    049         */
    050        public BivariateRealFunction interpolate(final double[] xval,
    051                                                              final double[] yval,
    052                                                              final double[][] zval)
    053            throws MathException, IllegalArgumentException {
    054            if (xval.length == 0 || yval.length == 0 || zval.length == 0) {
    055                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NO_DATA);
    056            }
    057            if (xval.length != zval.length) {
    058                throw new DimensionMismatchException(xval.length, zval.length);
    059            }
    060
    061            MathUtils.checkOrder(xval, OrderDirection.INCREASING, true);
    062            MathUtils.checkOrder(yval, OrderDirection.INCREASING, true);
    063
    064            final int xLen = xval.length;
    065            final int yLen = yval.length;
    066
    067            // Samples (first index is y-coordinate, i.e. subarray variable is x)
    068            // 0 <= i < xval.length
    069            // 0 <= j < yval.length
    070            // zX[j][i] = f(xval[i], yval[j])
    071            final double[][] zX = new double[yLen][xLen];
    072            for (int i = 0; i < xLen; i++) {
    073                if (zval[i].length != yLen) {
    074                    throw new DimensionMismatchException(zval[i].length, yLen);
    075                }
    076
    077                for (int j = 0; j < yLen; j++) {
    078                    zX[j][i] = zval[i][j];
    079                }
    080            }
    081
    082            final SplineInterpolator spInterpolator = new SplineInterpolator();
    083
    084            // For each line y[j] (0 <= j < yLen), construct a 1D spline with
    085            // respect to variable x
    086            final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
    087            for (int j = 0; j < yLen; j++) {
    088                ySplineX[j] = spInterpolator.interpolate(xval, zX[j]);
    089            }
    090
    091            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
    092            // values zY_1
    093            final double[][] zY_1 = new double[xLen][yLen];
    094            for (int j = 0; j < yLen; j++) {
    095                final PolynomialSplineFunction f = ySplineX[j];
    096                for (int i = 0; i < xLen; i++) {
    097                    zY_1[i][j] = f.value(xval[i]);
    098                }
    099            }
    100
    101            // For each line x[i] (0 <= i < xLen), construct a 1D spline with
    102            // respect to variable y generated by array zY_1[i]
    103            final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
    104            for (int i = 0; i < xLen; i++) {
    105                xSplineY[i] = spInterpolator.interpolate(yval, zY_1[i]);
    106            }
    107
    108            // For every knot (xval[i], yval[j]) of the grid, calculate corrected
    109            // values zY_2
    110            final double[][] zY_2 = new double[xLen][yLen];
    111            for (int i = 0; i < xLen; i++) {
    112                final PolynomialSplineFunction f = xSplineY[i];
    113                for (int j = 0; j < yLen; j++) {
    114                    zY_2[i][j] = f.value(yval[j]);
    115                }
    116            }
    117
    118            // Partial derivatives with respect to x at the grid knots
    119            final double[][] dZdX = new double[xLen][yLen];
    120            for (int j = 0; j < yLen; j++) {
    121                final UnivariateRealFunction f = ySplineX[j].derivative();
    122                for (int i = 0; i < xLen; i++) {
    123                    dZdX[i][j] = f.value(xval[i]);
    124                }
    125            }
    126
    127            // Partial derivatives with respect to y at the grid knots
    128            final double[][] dZdY = new double[xLen][yLen];
    129            for (int i = 0; i < xLen; i++) {
    130                final UnivariateRealFunction f = xSplineY[i].derivative();
    131                for (int j = 0; j < yLen; j++) {
    132                    dZdY[i][j] = f.value(yval[j]);
    133                }
    134            }
    135
    136            // Cross partial derivatives
    137            final double[][] dZdXdY = new double[xLen][yLen];
    138            for (int i = 0; i < xLen ; i++) {
    139                final int nI = nextIndex(i, xLen);
    140                final int pI = previousIndex(i);
    141                for (int j = 0; j < yLen; j++) {
    142                    final int nJ = nextIndex(j, yLen);
    143                    final int pJ = previousIndex(j);
    144                    dZdXdY[i][j] = (zY_2[nI][nJ] - zY_2[nI][pJ] -
    145                                    zY_2[pI][nJ] + zY_2[pI][pJ]) /
    146                        ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
    147                }
    148            }
    149
    150            // Create the interpolating splines
    151            return new BicubicSplineInterpolatingFunction(xval, yval, zY_2,
    152                                                          dZdX, dZdY, dZdXdY);
    153        }
    154
    155        /**
    156         * Compute the next index of an array, clipping if necessary.
    157         * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
    158         *
    159         * @param i Index
    160         * @param max Upper limit of the array
    161         * @return the next index
    162         */
    163        private int nextIndex(int i, int max) {
    164            final int index = i + 1;
    165            return index < max ? index : index - 1;
    166        }
    167        /**
    168         * Compute the previous index of an array, clipping if necessary.
    169         * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
    170         *
    171         * @param i Index
    172         * @return the previous index
    173         */
    174        private int previousIndex(int i) {
    175            final int index = i - 1;
    176            return index >= 0 ? index : 0;
    177        }
    178    }