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