1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.interpolation;
18
19import org.apache.commons.math.DimensionMismatchException;
20import org.apache.commons.math.FunctionEvaluationException;
21import org.apache.commons.math.analysis.BivariateRealFunction;
22import org.apache.commons.math.exception.NoDataException;
23import org.apache.commons.math.exception.OutOfRangeException;
24import org.apache.commons.math.util.MathUtils;
25
26/**
27 * Function that implements the
28 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
29 * bicubic spline interpolation</a>.
30 *
31 * @version $Revision$ $Date$
32 * @since 2.1
33 */
34public class BicubicSplineInterpolatingFunction
35    implements BivariateRealFunction {
36    /**
37     * Matrix to compute the spline coefficients from the function values
38     * and function derivatives values
39     */
40    private static final double[][] AINV = {
41        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
42        { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
43        { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
44        { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
45        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
46        { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
47        { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
48        { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
49        { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
50        { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
51        { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
52        { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
53        { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
54        { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
55        { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
56        { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
57    };
58
59    /** Samples x-coordinates */
60    private final double[] xval;
61    /** Samples y-coordinates */
62    private final double[] yval;
63    /** Set of cubic splines patching the whole data grid */
64    private final BicubicSplineFunction[][] splines;
65    /**
66     * Partial derivatives
67     * The value of the first index determines the kind of derivatives:
68     * 0 = first partial derivatives wrt x
69     * 1 = first partial derivatives wrt y
70     * 2 = second partial derivatives wrt x
71     * 3 = second partial derivatives wrt y
72     * 4 = cross partial derivatives
73     */
74    private BivariateRealFunction[][][] partialDerivatives = null;
75
76    /**
77     * @param x Sample values of the x-coordinate, in increasing order.
78     * @param y Sample values of the y-coordinate, in increasing order.
79     * @param f Values of the function on every grid point.
80     * @param dFdX Values of the partial derivative of function with respect
81     * to x on every grid point.
82     * @param dFdY Values of the partial derivative of function with respect
83     * to y on every grid point.
84     * @param d2FdXdY Values of the cross partial derivative of function on
85     * every grid point.
86     * @throws DimensionMismatchException if the various arrays do not contain
87     * the expected number of elements.
88     * @throws org.apache.commons.math.exception.NonMonotonousSequenceException
89     * if {@code x} or {@code y} are not strictly increasing.
90     * @throws NoDataException if any of the arrays has zero length.
91     */
92    public BicubicSplineInterpolatingFunction(double[] x,
93                                              double[] y,
94                                              double[][] f,
95                                              double[][] dFdX,
96                                              double[][] dFdY,
97                                              double[][] d2FdXdY)
98        throws DimensionMismatchException {
99        final int xLen = x.length;
100        final int yLen = y.length;
101
102        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
103            throw new NoDataException();
104        }
105        if (xLen != f.length) {
106            throw new DimensionMismatchException(xLen, f.length);
107        }
108        if (xLen != dFdX.length) {
109            throw new DimensionMismatchException(xLen, dFdX.length);
110        }
111        if (xLen != dFdY.length) {
112            throw new DimensionMismatchException(xLen, dFdY.length);
113        }
114        if (xLen != d2FdXdY.length) {
115            throw new DimensionMismatchException(xLen, d2FdXdY.length);
116        }
117
118        MathUtils.checkOrder(x);
119        MathUtils.checkOrder(y);
120
121        xval = x.clone();
122        yval = y.clone();
123
124        final int lastI = xLen - 1;
125        final int lastJ = yLen - 1;
126        splines = new BicubicSplineFunction[lastI][lastJ];
127
128        for (int i = 0; i < lastI; i++) {
129            if (f[i].length != yLen) {
130                throw new DimensionMismatchException(f[i].length, yLen);
131            }
132            if (dFdX[i].length != yLen) {
133                throw new DimensionMismatchException(dFdX[i].length, yLen);
134            }
135            if (dFdY[i].length != yLen) {
136                throw new DimensionMismatchException(dFdY[i].length, yLen);
137            }
138            if (d2FdXdY[i].length != yLen) {
139                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
140            }
141            final int ip1 = i + 1;
142            for (int j = 0; j < lastJ; j++) {
143                final int jp1 = j + 1;
144                final double[] beta = new double[] {
145                    f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
146                    dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
147                    dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
148                    d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
149                };
150
151                splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
152            }
153        }
154    }
155
156    /**
157     * {@inheritDoc}
158     */
159    public double value(double x, double y) {
160        final int i = searchIndex(x, xval);
161        if (i == -1) {
162            throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
163        }
164        final int j = searchIndex(y, yval);
165        if (j == -1) {
166            throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
167        }
168
169        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
170        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
171
172        return splines[i][j].value(xN, yN);
173    }
174
175    /**
176     * @param x x-coordinate.
177     * @param y y-coordinate.
178     * @return the value at point (x, y) of the first partial derivative with
179     * respect to x.
180     * @since 2.2
181     */
182    public double partialDerivativeX(double x, double y) {
183        return partialDerivative(0, x, y);
184    }
185    /**
186     * @param x x-coordinate.
187     * @param y y-coordinate.
188     * @return the value at point (x, y) of the first partial derivative with
189     * respect to y.
190     * @since 2.2
191     */
192    public double partialDerivativeY(double x, double y) {
193        return partialDerivative(1, x, y);
194    }
195    /**
196     * @param x x-coordinate.
197     * @param y y-coordinate.
198     * @return the value at point (x, y) of the second partial derivative with
199     * respect to x.
200     * @since 2.2
201     */
202    public double partialDerivativeXX(double x, double y) {
203        return partialDerivative(2, x, y);
204    }
205    /**
206     * @param x x-coordinate.
207     * @param y y-coordinate.
208     * @return the value at point (x, y) of the second partial derivative with
209     * respect to y.
210     * @since 2.2
211     */
212    public double partialDerivativeYY(double x, double y) {
213        return partialDerivative(3, x, y);
214    }
215    /**
216     * @param x x-coordinate.
217     * @param y y-coordinate.
218     * @return the value at point (x, y) of the second partial cross-derivative.
219     * @since 2.2
220     */
221    public double partialDerivativeXY(double x, double y) {
222        return partialDerivative(4, x, y);
223    }
224
225    /**
226     * @param which First index in {@link #partialDerivatives}.
227     * @param x x-coordinate.
228     * @param y y-coordinate.
229     * @return the value at point (x, y) of the selected partial derivative.
230     * @throws FunctionEvaluationException
231     */
232    private double partialDerivative(int which, double x, double y) {
233        if (partialDerivatives == null) {
234            computePartialDerivatives();
235        }
236
237        final int i = searchIndex(x, xval);
238        if (i == -1) {
239            throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
240        }
241        final int j = searchIndex(y, yval);
242        if (j == -1) {
243            throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
244        }
245
246        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
247        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
248
249        try {
250            return partialDerivatives[which][i][j].value(xN, yN);
251        } catch (FunctionEvaluationException fee) {
252            // this should never happen
253            throw new RuntimeException(fee);
254        }
255
256    }
257
258    /**
259     * Compute all partial derivatives.
260     */
261    private void computePartialDerivatives() {
262        final int lastI = xval.length - 1;
263        final int lastJ = yval.length - 1;
264        partialDerivatives = new BivariateRealFunction[5][lastI][lastJ];
265
266        for (int i = 0; i < lastI; i++) {
267            for (int j = 0; j < lastJ; j++) {
268                final BicubicSplineFunction f = splines[i][j];
269                partialDerivatives[0][i][j] = f.partialDerivativeX();
270                partialDerivatives[1][i][j] = f.partialDerivativeY();
271                partialDerivatives[2][i][j] = f.partialDerivativeXX();
272                partialDerivatives[3][i][j] = f.partialDerivativeYY();
273                partialDerivatives[4][i][j] = f.partialDerivativeXY();
274            }
275        }
276    }
277
278    /**
279     * @param c Coordinate.
280     * @param val Coordinate samples.
281     * @return the index in {@code val} corresponding to the interval
282     * containing {@code c}, or {@code -1} if {@code c} is out of the
283     * range defined by the end values of {@code val}.
284     */
285    private int searchIndex(double c, double[] val) {
286        if (c < val[0]) {
287            return -1;
288        }
289
290        final int max = val.length;
291        for (int i = 1; i < max; i++) {
292            if (c <= val[i]) {
293                return i - 1;
294            }
295        }
296
297        return -1;
298    }
299
300    /**
301     * Compute the spline coefficients from the list of function values and
302     * function partial derivatives values at the four corners of a grid
303     * element. They must be specified in the following order:
304     * <ul>
305     *  <li>f(0,0)</li>
306     *  <li>f(1,0)</li>
307     *  <li>f(0,1)</li>
308     *  <li>f(1,1)</li>
309     *  <li>f<sub>x</sub>(0,0)</li>
310     *  <li>f<sub>x</sub>(1,0)</li>
311     *  <li>f<sub>x</sub>(0,1)</li>
312     *  <li>f<sub>x</sub>(1,1)</li>
313     *  <li>f<sub>y</sub>(0,0)</li>
314     *  <li>f<sub>y</sub>(1,0)</li>
315     *  <li>f<sub>y</sub>(0,1)</li>
316     *  <li>f<sub>y</sub>(1,1)</li>
317     *  <li>f<sub>xy</sub>(0,0)</li>
318     *  <li>f<sub>xy</sub>(1,0)</li>
319     *  <li>f<sub>xy</sub>(0,1)</li>
320     *  <li>f<sub>xy</sub>(1,1)</li>
321     * </ul>
322     * where the subscripts indicate the partial derivative with respect to
323     * the corresponding variable(s).
324     *
325     * @param beta List of function values and function partial derivatives
326     * values.
327     * @return the spline coefficients.
328     */
329    private double[] computeSplineCoefficients(double[] beta) {
330        final double[] a = new double[16];
331
332        for (int i = 0; i < 16; i++) {
333            double result = 0;
334            final double[] row = AINV[i];
335            for (int j = 0; j < 16; j++) {
336                result += row[j] * beta[j];
337            }
338            a[i] = result;
339        }
340
341        return a;
342    }
343}
344
345/**
346 * 2D-spline function.
347 *
348 * @version $Revision$ $Date$
349 */
350class BicubicSplineFunction
351    implements BivariateRealFunction {
352
353    /** Number of points. */
354    private static final short N = 4;
355
356    /** Coefficients */
357    private final double[][] a;
358
359    /** First partial derivative along x. */
360    private BivariateRealFunction partialDerivativeX;
361
362    /** First partial derivative along y. */
363    private BivariateRealFunction partialDerivativeY;
364
365    /** Second partial derivative along x. */
366    private BivariateRealFunction partialDerivativeXX;
367
368    /** Second partial derivative along y. */
369    private BivariateRealFunction partialDerivativeYY;
370
371    /** Second crossed partial derivative. */
372    private BivariateRealFunction partialDerivativeXY;
373
374    /**
375     * Simple constructor.
376     * @param a Spline coefficients
377     */
378    public BicubicSplineFunction(double[] a) {
379        this.a = new double[N][N];
380        for (int i = 0; i < N; i++) {
381            for (int j = 0; j < N; j++) {
382                this.a[i][j] = a[i + N * j];
383            }
384        }
385    }
386
387    /**
388     * {@inheritDoc}
389     */
390    public double value(double x, double y) {
391        if (x < 0 || x > 1) {
392            throw new OutOfRangeException(x, 0, 1);
393        }
394        if (y < 0 || y > 1) {
395            throw new OutOfRangeException(y, 0, 1);
396        }
397
398        final double x2 = x * x;
399        final double x3 = x2 * x;
400        final double[] pX = {1, x, x2, x3};
401
402        final double y2 = y * y;
403        final double y3 = y2 * y;
404        final double[] pY = {1, y, y2, y3};
405
406        return apply(pX, pY, a);
407    }
408
409    /**
410     * Compute the value of the bicubic polynomial.
411     *
412     * @param pX Powers of the x-coordinate.
413     * @param pY Powers of the y-coordinate.
414     * @param coeff Spline coefficients.
415     * @return the interpolated value.
416     */
417    private double apply(double[] pX, double[] pY, double[][] coeff) {
418        double result = 0;
419        for (int i = 0; i < N; i++) {
420            for (int j = 0; j < N; j++) {
421                result += coeff[i][j] * pX[i] * pY[j];
422            }
423        }
424
425        return result;
426    }
427
428    /**
429     * @return the partial derivative wrt {@code x}.
430     */
431    public BivariateRealFunction partialDerivativeX() {
432        if (partialDerivativeX == null) {
433            computePartialDerivatives();
434        }
435
436        return partialDerivativeX;
437    }
438    /**
439     * @return the partial derivative wrt {@code y}.
440     */
441    public BivariateRealFunction partialDerivativeY() {
442        if (partialDerivativeY == null) {
443            computePartialDerivatives();
444        }
445
446        return partialDerivativeY;
447    }
448    /**
449     * @return the second partial derivative wrt {@code x}.
450     */
451    public BivariateRealFunction partialDerivativeXX() {
452        if (partialDerivativeXX == null) {
453            computePartialDerivatives();
454        }
455
456        return partialDerivativeXX;
457    }
458    /**
459     * @return the second partial derivative wrt {@code y}.
460     */
461    public BivariateRealFunction partialDerivativeYY() {
462        if (partialDerivativeYY == null) {
463            computePartialDerivatives();
464        }
465
466        return partialDerivativeYY;
467    }
468    /**
469     * @return the second partial cross-derivative.
470     */
471    public BivariateRealFunction partialDerivativeXY() {
472        if (partialDerivativeXY == null) {
473            computePartialDerivatives();
474        }
475
476        return partialDerivativeXY;
477    }
478
479    /**
480     * Compute all partial derivatives functions.
481     */
482    private void computePartialDerivatives() {
483        final double[][] aX = new double[N][N];
484        final double[][] aY = new double[N][N];
485        final double[][] aXX = new double[N][N];
486        final double[][] aYY = new double[N][N];
487        final double[][] aXY = new double[N][N];
488
489        for (int i = 0; i < N; i++) {
490            for (int j = 0; j < N; j++) {
491                final double c = a[i][j];
492                aX[i][j] = i * c;
493                aY[i][j] = j * c;
494                aXX[i][j] = (i - 1) * aX[i][j];
495                aYY[i][j] = (j - 1) * aY[i][j];
496                aXY[i][j] = j * aX[i][j];
497            }
498        }
499
500        partialDerivativeX = new BivariateRealFunction() {
501                public double value(double x, double y)  {
502                    final double x2 = x * x;
503                    final double[] pX = {0, 1, x, x2};
504
505                    final double y2 = y * y;
506                    final double y3 = y2 * y;
507                    final double[] pY = {1, y, y2, y3};
508
509                    return apply(pX, pY, aX);
510                }
511            };
512        partialDerivativeY = new BivariateRealFunction() {
513                public double value(double x, double y)  {
514                    final double x2 = x * x;
515                    final double x3 = x2 * x;
516                    final double[] pX = {1, x, x2, x3};
517
518                    final double y2 = y * y;
519                    final double[] pY = {0, 1, y, y2};
520
521                    return apply(pX, pY, aY);
522                }
523            };
524        partialDerivativeXX = new BivariateRealFunction() {
525                public double value(double x, double y)  {
526                    final double[] pX = {0, 0, 1, x};
527
528                    final double y2 = y * y;
529                    final double y3 = y2 * y;
530                    final double[] pY = {1, y, y2, y3};
531
532                    return apply(pX, pY, aXX);
533                }
534            };
535        partialDerivativeYY = new BivariateRealFunction() {
536                public double value(double x, double y)  {
537                    final double x2 = x * x;
538                    final double x3 = x2 * x;
539                    final double[] pX = {1, x, x2, x3};
540
541                    final double[] pY = {0, 0, 1, y};
542
543                    return apply(pX, pY, aYY);
544                }
545            };
546        partialDerivativeXY = new BivariateRealFunction() {
547                public double value(double x, double y)  {
548                    final double x2 = x * x;
549                    final double[] pX = {0, 1, x, x2};
550
551                    final double y2 = y * y;
552                    final double[] pY = {0, 1, y, y2};
553
554                    return apply(pX, pY, aXY);
555                }
556            };
557    }
558}
559