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 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import java.util.Arrays;
020import org.apache.commons.math3.analysis.BivariateFunction;
021import org.apache.commons.math3.exception.DimensionMismatchException;
022import org.apache.commons.math3.exception.NoDataException;
023import org.apache.commons.math3.exception.OutOfRangeException;
024import org.apache.commons.math3.exception.NonMonotonicSequenceException;
025import org.apache.commons.math3.util.MathArrays;
026
027/**
028 * Function that implements the
029 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
030 * bicubic spline interpolation</a>.
031 *
032 * @since 2.1
033 * @version $Id: BicubicSplineInterpolatingFunction.java 1512547 2013-08-10 01:13:38Z erans $
034 */
035public class BicubicSplineInterpolatingFunction
036    implements BivariateFunction {
037    /** Number of coefficients. */
038    private static final int NUM_COEFF = 16;
039    /**
040     * Matrix to compute the spline coefficients from the function values
041     * and function derivatives values
042     */
043    private static final double[][] AINV = {
044        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
045        { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
046        { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
047        { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
048        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
049        { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
050        { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
051        { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
052        { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
053        { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
054        { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
055        { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
056        { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
057        { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
058        { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
059        { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
060    };
061
062    /** Samples x-coordinates */
063    private final double[] xval;
064    /** Samples y-coordinates */
065    private final double[] yval;
066    /** Set of cubic splines patching the whole data grid */
067    private final BicubicSplineFunction[][] splines;
068    /**
069     * Partial derivatives
070     * The value of the first index determines the kind of derivatives:
071     * 0 = first partial derivatives wrt x
072     * 1 = first partial derivatives wrt y
073     * 2 = second partial derivatives wrt x
074     * 3 = second partial derivatives wrt y
075     * 4 = cross partial derivatives
076     */
077    private BivariateFunction[][][] partialDerivatives = null;
078
079    /**
080     * @param x Sample values of the x-coordinate, in increasing order.
081     * @param y Sample values of the y-coordinate, in increasing order.
082     * @param f Values of the function on every grid point.
083     * @param dFdX Values of the partial derivative of function with respect
084     * to x on every grid point.
085     * @param dFdY Values of the partial derivative of function with respect
086     * to y on every grid point.
087     * @param d2FdXdY Values of the cross partial derivative of function on
088     * every grid point.
089     * @throws DimensionMismatchException if the various arrays do not contain
090     * the expected number of elements.
091     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
092     * not strictly increasing.
093     * @throws NoDataException if any of the arrays has zero length.
094     */
095    public BicubicSplineInterpolatingFunction(double[] x,
096                                              double[] y,
097                                              double[][] f,
098                                              double[][] dFdX,
099                                              double[][] dFdY,
100                                              double[][] d2FdXdY)
101        throws DimensionMismatchException,
102               NoDataException,
103               NonMonotonicSequenceException {
104        final int xLen = x.length;
105        final int yLen = y.length;
106
107        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
108            throw new NoDataException();
109        }
110        if (xLen != f.length) {
111            throw new DimensionMismatchException(xLen, f.length);
112        }
113        if (xLen != dFdX.length) {
114            throw new DimensionMismatchException(xLen, dFdX.length);
115        }
116        if (xLen != dFdY.length) {
117            throw new DimensionMismatchException(xLen, dFdY.length);
118        }
119        if (xLen != d2FdXdY.length) {
120            throw new DimensionMismatchException(xLen, d2FdXdY.length);
121        }
122
123        MathArrays.checkOrder(x);
124        MathArrays.checkOrder(y);
125
126        xval = x.clone();
127        yval = y.clone();
128
129        final int lastI = xLen - 1;
130        final int lastJ = yLen - 1;
131        splines = new BicubicSplineFunction[lastI][lastJ];
132
133        for (int i = 0; i < lastI; i++) {
134            if (f[i].length != yLen) {
135                throw new DimensionMismatchException(f[i].length, yLen);
136            }
137            if (dFdX[i].length != yLen) {
138                throw new DimensionMismatchException(dFdX[i].length, yLen);
139            }
140            if (dFdY[i].length != yLen) {
141                throw new DimensionMismatchException(dFdY[i].length, yLen);
142            }
143            if (d2FdXdY[i].length != yLen) {
144                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
145            }
146            final int ip1 = i + 1;
147            for (int j = 0; j < lastJ; j++) {
148                final int jp1 = j + 1;
149                final double[] beta = new double[] {
150                    f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
151                    dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
152                    dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
153                    d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
154                };
155
156                splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
157            }
158        }
159    }
160
161    /**
162     * {@inheritDoc}
163     */
164    public double value(double x, double y)
165        throws OutOfRangeException {
166        final int i = searchIndex(x, xval);
167        final int j = searchIndex(y, yval);
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     * Indicates whether a point is within the interpolation range.
177     *
178     * @param x First coordinate.
179     * @param y Second coordinate.
180     * @return {@code true} if (x, y) is a valid point.
181     * @since 3.3
182     */
183    public boolean isValidPoint(double x, double y) {
184        if (x < xval[0] ||
185            x > xval[xval.length - 1] ||
186            y < yval[0] ||
187            y > yval[yval.length - 1]) {
188            return false;
189        } else {
190            return true;
191        }
192    }
193
194    /**
195     * @param x x-coordinate.
196     * @param y y-coordinate.
197     * @return the value at point (x, y) of the first partial derivative with
198     * respect to x.
199     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
200     * the range defined by the boundary values of {@code xval} (resp.
201     * {@code yval}).
202     */
203    public double partialDerivativeX(double x, double y)
204        throws OutOfRangeException {
205        return partialDerivative(0, x, y);
206    }
207    /**
208     * @param x x-coordinate.
209     * @param y y-coordinate.
210     * @return the value at point (x, y) of the first partial derivative with
211     * respect to y.
212     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
213     * the range defined by the boundary values of {@code xval} (resp.
214     * {@code yval}).
215     */
216    public double partialDerivativeY(double x, double y)
217        throws OutOfRangeException {
218        return partialDerivative(1, x, y);
219    }
220    /**
221     * @param x x-coordinate.
222     * @param y y-coordinate.
223     * @return the value at point (x, y) of the second partial derivative with
224     * respect to x.
225     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
226     * the range defined by the boundary values of {@code xval} (resp.
227     * {@code yval}).
228     */
229    public double partialDerivativeXX(double x, double y)
230        throws OutOfRangeException {
231        return partialDerivative(2, x, y);
232    }
233    /**
234     * @param x x-coordinate.
235     * @param y y-coordinate.
236     * @return the value at point (x, y) of the second partial derivative with
237     * respect to y.
238     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
239     * the range defined by the boundary values of {@code xval} (resp.
240     * {@code yval}).
241     */
242    public double partialDerivativeYY(double x, double y)
243        throws OutOfRangeException {
244        return partialDerivative(3, x, y);
245    }
246    /**
247     * @param x x-coordinate.
248     * @param y y-coordinate.
249     * @return the value at point (x, y) of the second partial cross-derivative.
250     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
251     * the range defined by the boundary values of {@code xval} (resp.
252     * {@code yval}).
253     */
254    public double partialDerivativeXY(double x, double y)
255        throws OutOfRangeException {
256        return partialDerivative(4, x, y);
257    }
258
259    /**
260     * @param which First index in {@link #partialDerivatives}.
261     * @param x x-coordinate.
262     * @param y y-coordinate.
263     * @return the value at point (x, y) of the selected partial derivative.
264     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
265     * the range defined by the boundary values of {@code xval} (resp.
266     * {@code yval}).
267     */
268    private double partialDerivative(int which, double x, double y)
269        throws OutOfRangeException {
270        if (partialDerivatives == null) {
271            computePartialDerivatives();
272        }
273
274        final int i = searchIndex(x, xval);
275        final int j = searchIndex(y, yval);
276
277        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
278        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
279
280        return partialDerivatives[which][i][j].value(xN, yN);
281    }
282
283    /**
284     * Compute all partial derivatives.
285     */
286    private void computePartialDerivatives() {
287        final int lastI = xval.length - 1;
288        final int lastJ = yval.length - 1;
289        partialDerivatives = new BivariateFunction[5][lastI][lastJ];
290
291        for (int i = 0; i < lastI; i++) {
292            for (int j = 0; j < lastJ; j++) {
293                final BicubicSplineFunction f = splines[i][j];
294                partialDerivatives[0][i][j] = f.partialDerivativeX();
295                partialDerivatives[1][i][j] = f.partialDerivativeY();
296                partialDerivatives[2][i][j] = f.partialDerivativeXX();
297                partialDerivatives[3][i][j] = f.partialDerivativeYY();
298                partialDerivatives[4][i][j] = f.partialDerivativeXY();
299            }
300        }
301    }
302
303    /**
304     * @param c Coordinate.
305     * @param val Coordinate samples.
306     * @return the index in {@code val} corresponding to the interval
307     * containing {@code c}.
308     * @throws OutOfRangeException if {@code c} is out of the
309     * range defined by the boundary values of {@code val}.
310     */
311    private int searchIndex(double c, double[] val) {
312        final int r = Arrays.binarySearch(val, c);
313
314        if (r == -1 ||
315            r == -val.length - 1) {
316            throw new OutOfRangeException(c, val[0], val[val.length - 1]);
317        }
318
319        if (r < 0) {
320            // "c" in within an interpolation sub-interval: Return the
321            // index of the sample at the lower end of the sub-interval.
322            return -r - 2;
323        }
324        final int last = val.length - 1;
325        if (r == last) {
326            // "c" is the last sample of the range: Return the index
327            // of the sample at the lower end of the last sub-interval.
328            return last - 1;
329        }
330
331        // "c" is another sample point.
332        return r;
333    }
334
335    /**
336     * Compute the spline coefficients from the list of function values and
337     * function partial derivatives values at the four corners of a grid
338     * element. They must be specified in the following order:
339     * <ul>
340     *  <li>f(0,0)</li>
341     *  <li>f(1,0)</li>
342     *  <li>f(0,1)</li>
343     *  <li>f(1,1)</li>
344     *  <li>f<sub>x</sub>(0,0)</li>
345     *  <li>f<sub>x</sub>(1,0)</li>
346     *  <li>f<sub>x</sub>(0,1)</li>
347     *  <li>f<sub>x</sub>(1,1)</li>
348     *  <li>f<sub>y</sub>(0,0)</li>
349     *  <li>f<sub>y</sub>(1,0)</li>
350     *  <li>f<sub>y</sub>(0,1)</li>
351     *  <li>f<sub>y</sub>(1,1)</li>
352     *  <li>f<sub>xy</sub>(0,0)</li>
353     *  <li>f<sub>xy</sub>(1,0)</li>
354     *  <li>f<sub>xy</sub>(0,1)</li>
355     *  <li>f<sub>xy</sub>(1,1)</li>
356     * </ul>
357     * where the subscripts indicate the partial derivative with respect to
358     * the corresponding variable(s).
359     *
360     * @param beta List of function values and function partial derivatives
361     * values.
362     * @return the spline coefficients.
363     */
364    private double[] computeSplineCoefficients(double[] beta) {
365        final double[] a = new double[NUM_COEFF];
366
367        for (int i = 0; i < NUM_COEFF; i++) {
368            double result = 0;
369            final double[] row = AINV[i];
370            for (int j = 0; j < NUM_COEFF; j++) {
371                result += row[j] * beta[j];
372            }
373            a[i] = result;
374        }
375
376        return a;
377    }
378}
379
380/**
381 * 2D-spline function.
382 *
383 * @version $Id: BicubicSplineInterpolatingFunction.java 1512547 2013-08-10 01:13:38Z erans $
384 */
385class BicubicSplineFunction
386    implements BivariateFunction {
387
388    /** Number of points. */
389    private static final short N = 4;
390
391    /** Coefficients */
392    private final double[][] a;
393
394    /** First partial derivative along x. */
395    private BivariateFunction partialDerivativeX;
396
397    /** First partial derivative along y. */
398    private BivariateFunction partialDerivativeY;
399
400    /** Second partial derivative along x. */
401    private BivariateFunction partialDerivativeXX;
402
403    /** Second partial derivative along y. */
404    private BivariateFunction partialDerivativeYY;
405
406    /** Second crossed partial derivative. */
407    private BivariateFunction partialDerivativeXY;
408
409    /**
410     * Simple constructor.
411     * @param a Spline coefficients
412     */
413    public BicubicSplineFunction(double[] a) {
414        this.a = new double[N][N];
415        for (int i = 0; i < N; i++) {
416            for (int j = 0; j < N; j++) {
417                this.a[i][j] = a[i * N + j];
418            }
419        }
420    }
421
422    /**
423     * {@inheritDoc}
424     */
425    public double value(double x, double y) {
426        if (x < 0 || x > 1) {
427            throw new OutOfRangeException(x, 0, 1);
428        }
429        if (y < 0 || y > 1) {
430            throw new OutOfRangeException(y, 0, 1);
431        }
432
433        final double x2 = x * x;
434        final double x3 = x2 * x;
435        final double[] pX = {1, x, x2, x3};
436
437        final double y2 = y * y;
438        final double y3 = y2 * y;
439        final double[] pY = {1, y, y2, y3};
440
441        return apply(pX, pY, a);
442    }
443
444    /**
445     * Compute the value of the bicubic polynomial.
446     *
447     * @param pX Powers of the x-coordinate.
448     * @param pY Powers of the y-coordinate.
449     * @param coeff Spline coefficients.
450     * @return the interpolated value.
451     */
452    private double apply(double[] pX, double[] pY, double[][] coeff) {
453        double result = 0;
454        for (int i = 0; i < N; i++) {
455            for (int j = 0; j < N; j++) {
456                result += coeff[i][j] * pX[i] * pY[j];
457            }
458        }
459
460        return result;
461    }
462
463    /**
464     * @return the partial derivative wrt {@code x}.
465     */
466    public BivariateFunction partialDerivativeX() {
467        if (partialDerivativeX == null) {
468            computePartialDerivatives();
469        }
470
471        return partialDerivativeX;
472    }
473    /**
474     * @return the partial derivative wrt {@code y}.
475     */
476    public BivariateFunction partialDerivativeY() {
477        if (partialDerivativeY == null) {
478            computePartialDerivatives();
479        }
480
481        return partialDerivativeY;
482    }
483    /**
484     * @return the second partial derivative wrt {@code x}.
485     */
486    public BivariateFunction partialDerivativeXX() {
487        if (partialDerivativeXX == null) {
488            computePartialDerivatives();
489        }
490
491        return partialDerivativeXX;
492    }
493    /**
494     * @return the second partial derivative wrt {@code y}.
495     */
496    public BivariateFunction partialDerivativeYY() {
497        if (partialDerivativeYY == null) {
498            computePartialDerivatives();
499        }
500
501        return partialDerivativeYY;
502    }
503    /**
504     * @return the second partial cross-derivative.
505     */
506    public BivariateFunction partialDerivativeXY() {
507        if (partialDerivativeXY == null) {
508            computePartialDerivatives();
509        }
510
511        return partialDerivativeXY;
512    }
513
514    /**
515     * Compute all partial derivatives functions.
516     */
517    private void computePartialDerivatives() {
518        final double[][] aX = new double[N][N];
519        final double[][] aY = new double[N][N];
520        final double[][] aXX = new double[N][N];
521        final double[][] aYY = new double[N][N];
522        final double[][] aXY = new double[N][N];
523
524        for (int i = 0; i < N; i++) {
525            for (int j = 0; j < N; j++) {
526                final double c = a[i][j];
527                aX[i][j] = i * c;
528                aY[i][j] = j * c;
529                aXX[i][j] = (i - 1) * aX[i][j];
530                aYY[i][j] = (j - 1) * aY[i][j];
531                aXY[i][j] = j * aX[i][j];
532            }
533        }
534
535        partialDerivativeX = new BivariateFunction() {
536                public double value(double x, double y)  {
537                    final double x2 = x * x;
538                    final double[] pX = {0, 1, x, x2};
539
540                    final double y2 = y * y;
541                    final double y3 = y2 * y;
542                    final double[] pY = {1, y, y2, y3};
543
544                    return apply(pX, pY, aX);
545                }
546            };
547        partialDerivativeY = new BivariateFunction() {
548                public double value(double x, double y)  {
549                    final double x2 = x * x;
550                    final double x3 = x2 * x;
551                    final double[] pX = {1, x, x2, x3};
552
553                    final double y2 = y * y;
554                    final double[] pY = {0, 1, y, y2};
555
556                    return apply(pX, pY, aY);
557                }
558            };
559        partialDerivativeXX = new BivariateFunction() {
560                public double value(double x, double y)  {
561                    final double[] pX = {0, 0, 1, x};
562
563                    final double y2 = y * y;
564                    final double y3 = y2 * y;
565                    final double[] pY = {1, y, y2, y3};
566
567                    return apply(pX, pY, aXX);
568                }
569            };
570        partialDerivativeYY = new BivariateFunction() {
571                public double value(double x, double y)  {
572                    final double x2 = x * x;
573                    final double x3 = x2 * x;
574                    final double[] pX = {1, x, x2, x3};
575
576                    final double[] pY = {0, 0, 1, y};
577
578                    return apply(pX, pY, aYY);
579                }
580            };
581        partialDerivativeXY = new BivariateFunction() {
582                public double value(double x, double y)  {
583                    final double x2 = x * x;
584                    final double[] pX = {0, 1, x, x2};
585
586                    final double y2 = y * y;
587                    final double[] pY = {0, 1, y, y2};
588
589                    return apply(pX, pY, aXY);
590                }
591            };
592    }
593}