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