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.util.MathUtils;
020    import org.apache.commons.math.MathRuntimeException;
021    import org.apache.commons.math.DimensionMismatchException;
022    import org.apache.commons.math.analysis.BivariateRealFunction;
023    
024    /**
025     * Function that implements the
026     * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
027     * bicubic spline interpolation</a>.
028     *
029     * @version $Revision$ $Date$
030     * @since 2.1
031     */
032    public class BicubicSplineInterpolatingFunction
033        implements BivariateRealFunction {
034        /**
035         * Matrix to compute the spline coefficients from the function values
036         * and function derivatives values
037         */
038        private static final double[][] AINV = {
039            { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
040            { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
041            { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
042            { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
043            { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
044            { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
045            { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
046            { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
047            { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
048            { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
049            { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
050            { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
051            { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
052            { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
053            { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
054            { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
055        };
056    
057        /** Samples x-coordinates */
058        private final double[] xval;
059        /** Samples y-coordinates */
060        private final double[] yval;
061        /** Set of cubic splines pacthing the whole data grid */
062        private final BicubicSplineFunction[][] splines;
063    
064        /**
065         * @param x Sample values of the x-coordinate, in increasing order
066         * @param y Sample values of the y-coordinate, in increasing order
067         * @param z Values of the function on every grid point
068         * @param dZdX Values of the partial derivative of function with respect
069         * to x on every grid point
070         * @param dZdY Values of the partial derivative of function with respect
071         * to y on every grid point
072         * @param dZdXdY Values of the cross partial derivative of function on
073         * every grid point
074         * @throws DimensionMismatchException if the various arrays do not contain
075         * the expected number of elements.
076         * @throws IllegalArgumentException if {@code x} or {@code y} are not strictly
077         * increasing.
078         */
079        public BicubicSplineInterpolatingFunction(double[] x,
080                                                  double[] y,
081                                                  double[][] z,
082                                                  double[][] dZdX,
083                                                  double[][] dZdY,
084                                                  double[][] dZdXdY)
085            throws DimensionMismatchException {
086            final int xLen = x.length;
087            final int yLen = y.length;
088    
089            if (xLen == 0 || yLen == 0 || z.length == 0 || z[0].length == 0) {
090                throw MathRuntimeException.createIllegalArgumentException("no data");
091            }
092            if (xLen != z.length) {
093                throw new DimensionMismatchException(xLen, z.length);
094            }
095            if (xLen != dZdX.length) {
096                throw new DimensionMismatchException(xLen, dZdX.length);
097            }
098            if (xLen != dZdY.length) {
099                throw new DimensionMismatchException(xLen, dZdY.length);
100            }
101            if (xLen != dZdXdY.length) {
102                throw new DimensionMismatchException(xLen, dZdXdY.length);
103            }
104    
105            MathUtils.checkOrder(x, 1, true);
106            MathUtils.checkOrder(y, 1, true);
107    
108            xval = x.clone();
109            yval = y.clone();
110    
111            final int lastI = xLen - 1;
112            final int lastJ = yLen - 1;
113            splines = new BicubicSplineFunction[lastI][lastJ];
114    
115            for (int i = 0; i < lastI; i++) {
116                if (z[i].length != yLen) {
117                    throw new DimensionMismatchException(z[i].length, yLen);
118                }
119                if (dZdX[i].length != yLen) {
120                    throw new DimensionMismatchException(dZdX[i].length, yLen);
121                }
122                if (dZdY[i].length != yLen) {
123                    throw new DimensionMismatchException(dZdY[i].length, yLen);
124                }
125                if (dZdXdY[i].length != yLen) {
126                    throw new DimensionMismatchException(dZdXdY[i].length, yLen);
127                }
128                final int ip1 = i + 1;
129                for (int j = 0; j < lastJ; j++) {
130                    final int jp1 = j + 1;
131                    final double[] beta = new double[] {
132                        z[i][j],      z[ip1][j],      z[i][jp1],      z[ip1][jp1],
133                        dZdX[i][j],   dZdX[ip1][j],   dZdX[i][jp1],   dZdX[ip1][jp1],
134                        dZdY[i][j],   dZdY[ip1][j],   dZdY[i][jp1],   dZdY[ip1][jp1],
135                        dZdXdY[i][j], dZdXdY[ip1][j], dZdXdY[i][jp1], dZdXdY[ip1][jp1]
136                    };
137    
138                    splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
139                }
140            }
141        }
142    
143        /**
144         * {@inheritDoc}
145         */
146        public double value(double x, double y) {
147            final int i = searchIndex(x, xval);
148            if (i == -1) {
149                throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
150                                                                          x, xval[0], xval[xval.length - 1]);
151            }
152            final int j = searchIndex(y, yval);
153            if (j == -1) {
154                throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
155                                                                          y, yval[0], yval[yval.length - 1]);
156            }
157    
158            final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
159            final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
160    
161            return splines[i][j].value(xN, yN);
162        }
163    
164        /**
165         * @param c coordinate
166         * @param val coordinate samples
167         * @return the index in {@code val} corresponding to the interval
168         * containing {@code c}, or {@code -1} if {@code c} is out of the
169         * range defined by the end values of {@code val}
170         */
171        private int searchIndex(double c, double[] val) {
172            if (c < val[0]) {
173                return -1;
174            }
175    
176            int max = val.length;
177            for (int i = 1; i < max; i++) {
178                if (c <= val[i]) {
179                    return i - 1;
180                }
181            }
182    
183            return -1;
184        }
185    
186        /**
187         * Compute the spline coefficients from the list of function values and
188         * function partial derivatives values at the four corners of a grid
189         * element. They must be specified in the following order:
190         * <ul>
191         *  <li>f(0,0)</li>
192         *  <li>f(1,0)</li>
193         *  <li>f(0,1)</li>
194         *  <li>f(1,1)</li>
195         *  <li>fx(0,0)</li>
196         *  <li>fx(1,0)</li>
197         *  <li>fx(0,1)</li>
198         *  <li>fx(1,1)</li>
199         *  <li>fy(0,0)</li>
200         *  <li>fy(1,0)</li>
201         *  <li>fy(0,1)</li>
202         *  <li>fy(1,1)</li>
203         *  <li>fxy(0,0)</li>
204         *  <li>fxy(1,0)</li>
205         *  <li>fxy(0,1)</li>
206         *  <li>fxy(1,1)</li>
207         * </ul>
208         * @param beta List of function values and function partial derivatives
209         * values
210         * @return the spline coefficients
211         */
212        private double[] computeSplineCoefficients(double[] beta) {
213            final double[] a = new double[16];
214    
215            for (int i = 0; i < 16; i++) {
216                double result = 0;
217                final double[] row = AINV[i];
218                for (int j = 0; j < 16; j++) {
219                    result += row[j] * beta[j];
220                }
221                a[i] = result;
222            }
223    
224            return a;
225        }
226    }
227    /**
228     * 2D-spline function.
229     *
230     * @version $Revision$ $Date$
231     */
232    class BicubicSplineFunction
233        implements BivariateRealFunction {
234    //CHECKSTYLE: stop MultipleVariableDeclarations
235        /** Coefficients */
236        private final double
237            a00, a01, a02, a03,
238            a10, a11, a12, a13,
239            a20, a21, a22, a23,
240            a30, a31, a32, a33;
241    //CHECKSTYLE: resume MultipleVariableDeclarations
242    
243        /**
244         * @param a Spline coefficients
245         */
246        public BicubicSplineFunction(double[] a) {
247            this.a00 = a[0];
248            this.a10 = a[1];
249            this.a20 = a[2];
250            this.a30 = a[3];
251            this.a01 = a[4];
252            this.a11 = a[5];
253            this.a21 = a[6];
254            this.a31 = a[7];
255            this.a02 = a[8];
256            this.a12 = a[9];
257            this.a22 = a[10];
258            this.a32 = a[11];
259            this.a03 = a[12];
260            this.a13 = a[13];
261            this.a23 = a[14];
262            this.a33 = a[15];
263        }
264    
265        /**
266         * @param x x-coordinate of the interpolation point
267         * @param y y-coordinate of the interpolation point
268         * @return the interpolated value.
269         */
270        public double value(double x, double y) {
271            if (x < 0 || x > 1) {
272                throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
273                                                                          x, 0, 1);
274            }
275            if (y < 0 || y > 1) {
276                throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
277                                                                          y, 0, 1);
278            }
279    
280            final double x2 = x * x;
281            final double x3 = x2 * x;
282            final double y2 = y * y;
283            final double y3 = y2 * y;
284    
285            return a00 + a01 * y + a02 * y2 + a03 * y3 +
286                a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 +
287                a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 +
288                a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * y3;
289        }
290    }