    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.analysis.UnivariateFunction;
    020    import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
    021    import org.apache.commons.math3.exception.DimensionMismatchException;
    022    import org.apache.commons.math3.exception.NoDataException;
    023    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
    024    import org.apache.commons.math3.util.MathArrays;
    025
    026    /**
    027     * Generates a bicubic interpolating function.
    028     *
    029     * @version $Id: BicubicSplineInterpolator.java 1385313 2012-09-16 16:35:23Z tn $
    030     * @since 2.2
    031     */
    032    public class BicubicSplineInterpolator
    033        implements BivariateGridInterpolator {
    034        /**
    035         * {@inheritDoc}
    036         */
    037        public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
    038                                                              final double[] yval,
    039                                                              final double[][] fval)
    040            throws NoDataException,
    041                   DimensionMismatchException,
    042                   NonMonotonicSequenceException {
    043            if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
    044                throw new NoDataException();
    045            }
    046            if (xval.length != fval.length) {
    047                throw new DimensionMismatchException(xval.length, fval.length);
    048            }
    049
    050            MathArrays.checkOrder(xval);
    051            MathArrays.checkOrder(yval);
    052
    053            final int xLen = xval.length;
    054            final int yLen = yval.length;
    055
    056            // Samples (first index is y-coordinate, i.e. subarray variable is x)
    057            // 0 <= i < xval.length
    058            // 0 <= j < yval.length
    059            // fX[j][i] = f(xval[i], yval[j])
    060            final double[][] fX = new double[yLen][xLen];
    061            for (int i = 0; i < xLen; i++) {
    062                if (fval[i].length != yLen) {
    063                    throw new DimensionMismatchException(fval[i].length, yLen);
    064                }
    065
    066                for (int j = 0; j < yLen; j++) {
    067                    fX[j][i] = fval[i][j];
    068                }
    069            }
    070
    071            final SplineInterpolator spInterpolator = new SplineInterpolator();
    072
    073            // For each line y[j] (0 <= j < yLen), construct a 1D spline with
    074            // respect to variable x
    075            final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
    076            for (int j = 0; j < yLen; j++) {
    077                ySplineX[j] = spInterpolator.interpolate(xval, fX[j]);
    078            }
    079
    080            // For each line x[i] (0 <= i < xLen), construct a 1D spline with
    081            // respect to variable y generated by array fY_1[i]
    082            final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
    083            for (int i = 0; i < xLen; i++) {
    084                xSplineY[i] = spInterpolator.interpolate(yval, fval[i]);
    085            }
    086
    087            // Partial derivatives with respect to x at the grid knots
    088            final double[][] dFdX = new double[xLen][yLen];
    089            for (int j = 0; j < yLen; j++) {
    090                final UnivariateFunction f = ySplineX[j].derivative();
    091                for (int i = 0; i < xLen; i++) {
    092                    dFdX[i][j] = f.value(xval[i]);
    093                }
    094            }
    095
    096            // Partial derivatives with respect to y at the grid knots
    097            final double[][] dFdY = new double[xLen][yLen];
    098            for (int i = 0; i < xLen; i++) {
    099                final UnivariateFunction f = xSplineY[i].derivative();
    100                for (int j = 0; j < yLen; j++) {
    101                    dFdY[i][j] = f.value(yval[j]);
    102                }
    103            }
    104
    105            // Cross partial derivatives
    106            final double[][] d2FdXdY = new double[xLen][yLen];
    107            for (int i = 0; i < xLen ; i++) {
    108                final int nI = nextIndex(i, xLen);
    109                final int pI = previousIndex(i);
    110                for (int j = 0; j < yLen; j++) {
    111                    final int nJ = nextIndex(j, yLen);
    112                    final int pJ = previousIndex(j);
    113                    d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] -
    114                                     fval[pI][nJ] + fval[pI][pJ]) /
    115                        ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
    116                }
    117            }
    118
    119            // Create the interpolating splines
    120            return new BicubicSplineInterpolatingFunction(xval, yval, fval,
    121                                                          dFdX, dFdY, d2FdXdY);
    122        }
    123
    124        /**
    125         * Computes the next index of an array, clipping if necessary.
    126         * It is assumed (but not checked) that {@code i >= 0}.
    127         *
    128         * @param i Index.
    129         * @param max Upper limit of the array.
    130         * @return the next index.
    131         */
    132        private int nextIndex(int i, int max) {
    133            final int index = i + 1;
    134            return index < max ? index : index - 1;
    135        }
    136        /**
    137         * Computes the previous index of an array, clipping if necessary.
    138         * It is assumed (but not checked) that {@code i} is smaller than the size
    139         * of the array.
    140         *
    141         * @param i Index.
    142         * @return the previous index.
    143         */
    144        private int previousIndex(int i) {
    145            final int index = i - 1;
    146            return index >= 0 ? index : 0;
    147        }
    148    }