View Javadoc
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   */
17  package org.apache.commons.math3.analysis.interpolation;
18  
19  import java.util.Arrays;
20  import org.apache.commons.math3.analysis.BivariateFunction;
21  import org.apache.commons.math3.exception.DimensionMismatchException;
22  import org.apache.commons.math3.exception.NoDataException;
23  import org.apache.commons.math3.exception.OutOfRangeException;
24  import org.apache.commons.math3.exception.NonMonotonicSequenceException;
25  import org.apache.commons.math3.util.MathArrays;
26  
27  /**
28   * Function that implements the
29   * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
30   * bicubic spline interpolation</a>. Due to numerical accuracy issues this should not
31   * be used.
32   *
33   * @since 2.1
34   * @deprecated as of 3.4 replaced by
35   * {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction}
36   */
37  @Deprecated
38  public class BicubicSplineInterpolatingFunction
39      implements BivariateFunction {
40      /** Number of coefficients. */
41      private static final int NUM_COEFF = 16;
42      /**
43       * Matrix to compute the spline coefficients from the function values
44       * and function derivatives values
45       */
46      private static final double[][] AINV = {
47          { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
48          { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
49          { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
50          { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
51          { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
52          { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
53          { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
54          { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
55          { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
56          { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
57          { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
58          { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
59          { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
60          { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
61          { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
62          { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
63      };
64  
65      /** Samples x-coordinates */
66      private final double[] xval;
67      /** Samples y-coordinates */
68      private final double[] yval;
69      /** Set of cubic splines patching the whole data grid */
70      private final BicubicSplineFunction[][] splines;
71      /**
72       * Partial derivatives.
73       * The value of the first index determines the kind of derivatives:
74       * 0 = first partial derivatives wrt x
75       * 1 = first partial derivatives wrt y
76       * 2 = second partial derivatives wrt x
77       * 3 = second partial derivatives wrt y
78       * 4 = cross partial derivatives
79       */
80      private final BivariateFunction[][][] partialDerivatives;
81  
82      /**
83       * @param x Sample values of the x-coordinate, in increasing order.
84       * @param y Sample values of the y-coordinate, in increasing order.
85       * @param f Values of the function on every grid point.
86       * @param dFdX Values of the partial derivative of function with respect
87       * to x on every grid point.
88       * @param dFdY Values of the partial derivative of function with respect
89       * to y on every grid point.
90       * @param d2FdXdY Values of the cross partial derivative of function on
91       * every grid point.
92       * @throws DimensionMismatchException if the various arrays do not contain
93       * the expected number of elements.
94       * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
95       * not strictly increasing.
96       * @throws NoDataException if any of the arrays has zero length.
97       */
98      public BicubicSplineInterpolatingFunction(double[] x,
99                                                double[] y,
100                                               double[][] f,
101                                               double[][] dFdX,
102                                               double[][] dFdY,
103                                               double[][] d2FdXdY)
104         throws DimensionMismatchException,
105                NoDataException,
106                NonMonotonicSequenceException {
107         this(x, y, f, dFdX, dFdY, d2FdXdY, false);
108     }
109 
110     /**
111      * @param x Sample values of the x-coordinate, in increasing order.
112      * @param y Sample values of the y-coordinate, in increasing order.
113      * @param f Values of the function on every grid point.
114      * @param dFdX Values of the partial derivative of function with respect
115      * to x on every grid point.
116      * @param dFdY Values of the partial derivative of function with respect
117      * to y on every grid point.
118      * @param d2FdXdY Values of the cross partial derivative of function on
119      * every grid point.
120      * @param initializeDerivatives Whether to initialize the internal data
121      * needed for calling any of the methods that compute the partial derivatives
122      * this function.
123      * @throws DimensionMismatchException if the various arrays do not contain
124      * the expected number of elements.
125      * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
126      * not strictly increasing.
127      * @throws NoDataException if any of the arrays has zero length.
128      *
129      * @see #partialDerivativeX(double,double)
130      * @see #partialDerivativeY(double,double)
131      * @see #partialDerivativeXX(double,double)
132      * @see #partialDerivativeYY(double,double)
133      * @see #partialDerivativeXY(double,double)
134      */
135     public BicubicSplineInterpolatingFunction(double[] x,
136                                               double[] y,
137                                               double[][] f,
138                                               double[][] dFdX,
139                                               double[][] dFdY,
140                                               double[][] d2FdXdY,
141                                               boolean initializeDerivatives)
142         throws DimensionMismatchException,
143                NoDataException,
144                NonMonotonicSequenceException {
145         final int xLen = x.length;
146         final int yLen = y.length;
147 
148         if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
149             throw new NoDataException();
150         }
151         if (xLen != f.length) {
152             throw new DimensionMismatchException(xLen, f.length);
153         }
154         if (xLen != dFdX.length) {
155             throw new DimensionMismatchException(xLen, dFdX.length);
156         }
157         if (xLen != dFdY.length) {
158             throw new DimensionMismatchException(xLen, dFdY.length);
159         }
160         if (xLen != d2FdXdY.length) {
161             throw new DimensionMismatchException(xLen, d2FdXdY.length);
162         }
163 
164         MathArrays.checkOrder(x);
165         MathArrays.checkOrder(y);
166 
167         xval = x.clone();
168         yval = y.clone();
169 
170         final int lastI = xLen - 1;
171         final int lastJ = yLen - 1;
172         splines = new BicubicSplineFunction[lastI][lastJ];
173 
174         for (int i = 0; i < lastI; i++) {
175             if (f[i].length != yLen) {
176                 throw new DimensionMismatchException(f[i].length, yLen);
177             }
178             if (dFdX[i].length != yLen) {
179                 throw new DimensionMismatchException(dFdX[i].length, yLen);
180             }
181             if (dFdY[i].length != yLen) {
182                 throw new DimensionMismatchException(dFdY[i].length, yLen);
183             }
184             if (d2FdXdY[i].length != yLen) {
185                 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
186             }
187             final int ip1 = i + 1;
188             for (int j = 0; j < lastJ; j++) {
189                 final int jp1 = j + 1;
190                 final double[] beta = new double[] {
191                     f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
192                     dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
193                     dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
194                     d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
195                 };
196 
197                 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
198                                                           initializeDerivatives);
199             }
200         }
201 
202         if (initializeDerivatives) {
203             // Compute all partial derivatives.
204             partialDerivatives = new BivariateFunction[5][lastI][lastJ];
205 
206             for (int i = 0; i < lastI; i++) {
207                 for (int j = 0; j < lastJ; j++) {
208                     final BicubicSplineFunction bcs = splines[i][j];
209                     partialDerivatives[0][i][j] = bcs.partialDerivativeX();
210                     partialDerivatives[1][i][j] = bcs.partialDerivativeY();
211                     partialDerivatives[2][i][j] = bcs.partialDerivativeXX();
212                     partialDerivatives[3][i][j] = bcs.partialDerivativeYY();
213                     partialDerivatives[4][i][j] = bcs.partialDerivativeXY();
214                 }
215             }
216         } else {
217             // Partial derivative methods cannot be used.
218             partialDerivatives = null;
219         }
220     }
221 
222     /**
223      * {@inheritDoc}
224      */
225     public double value(double x, double y)
226         throws OutOfRangeException {
227         final int i = searchIndex(x, xval);
228         final int j = searchIndex(y, yval);
229 
230         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
231         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
232 
233         return splines[i][j].value(xN, yN);
234     }
235 
236     /**
237      * Indicates whether a point is within the interpolation range.
238      *
239      * @param x First coordinate.
240      * @param y Second coordinate.
241      * @return {@code true} if (x, y) is a valid point.
242      * @since 3.3
243      */
244     public boolean isValidPoint(double x, double y) {
245         if (x < xval[0] ||
246             x > xval[xval.length - 1] ||
247             y < yval[0] ||
248             y > yval[yval.length - 1]) {
249             return false;
250         } else {
251             return true;
252         }
253     }
254 
255     /**
256      * @param x x-coordinate.
257      * @param y y-coordinate.
258      * @return the value at point (x, y) of the first partial derivative with
259      * respect to x.
260      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
261      * the range defined by the boundary values of {@code xval} (resp.
262      * {@code yval}).
263      * @throws NullPointerException if the internal data were not initialized
264      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
265      *             double[][],double[][],double[][],boolean) constructor}).
266      */
267     public double partialDerivativeX(double x, double y)
268         throws OutOfRangeException {
269         return partialDerivative(0, x, y);
270     }
271     /**
272      * @param x x-coordinate.
273      * @param y y-coordinate.
274      * @return the value at point (x, y) of the first partial derivative with
275      * respect to y.
276      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
277      * the range defined by the boundary values of {@code xval} (resp.
278      * {@code yval}).
279      * @throws NullPointerException if the internal data were not initialized
280      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
281      *             double[][],double[][],double[][],boolean) constructor}).
282      */
283     public double partialDerivativeY(double x, double y)
284         throws OutOfRangeException {
285         return partialDerivative(1, x, y);
286     }
287     /**
288      * @param x x-coordinate.
289      * @param y y-coordinate.
290      * @return the value at point (x, y) of the second partial derivative with
291      * respect to x.
292      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
293      * the range defined by the boundary values of {@code xval} (resp.
294      * {@code yval}).
295      * @throws NullPointerException if the internal data were not initialized
296      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
297      *             double[][],double[][],double[][],boolean) constructor}).
298      */
299     public double partialDerivativeXX(double x, double y)
300         throws OutOfRangeException {
301         return partialDerivative(2, x, y);
302     }
303     /**
304      * @param x x-coordinate.
305      * @param y y-coordinate.
306      * @return the value at point (x, y) of the second partial derivative with
307      * respect to y.
308      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
309      * the range defined by the boundary values of {@code xval} (resp.
310      * {@code yval}).
311      * @throws NullPointerException if the internal data were not initialized
312      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
313      *             double[][],double[][],double[][],boolean) constructor}).
314      */
315     public double partialDerivativeYY(double x, double y)
316         throws OutOfRangeException {
317         return partialDerivative(3, x, y);
318     }
319     /**
320      * @param x x-coordinate.
321      * @param y y-coordinate.
322      * @return the value at point (x, y) of the second partial cross-derivative.
323      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
324      * the range defined by the boundary values of {@code xval} (resp.
325      * {@code yval}).
326      * @throws NullPointerException if the internal data were not initialized
327      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
328      *             double[][],double[][],double[][],boolean) constructor}).
329      */
330     public double partialDerivativeXY(double x, double y)
331         throws OutOfRangeException {
332         return partialDerivative(4, x, y);
333     }
334 
335     /**
336      * @param which First index in {@link #partialDerivatives}.
337      * @param x x-coordinate.
338      * @param y y-coordinate.
339      * @return the value at point (x, y) of the selected partial derivative.
340      * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
341      * the range defined by the boundary values of {@code xval} (resp.
342      * {@code yval}).
343      * @throws NullPointerException if the internal data were not initialized
344      * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
345      *             double[][],double[][],double[][],boolean) constructor}).
346      */
347     private double partialDerivative(int which, double x, double y)
348         throws OutOfRangeException {
349         final int i = searchIndex(x, xval);
350         final int j = searchIndex(y, yval);
351 
352         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
353         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
354 
355         return partialDerivatives[which][i][j].value(xN, yN);
356     }
357 
358     /**
359      * @param c Coordinate.
360      * @param val Coordinate samples.
361      * @return the index in {@code val} corresponding to the interval
362      * containing {@code c}.
363      * @throws OutOfRangeException if {@code c} is out of the
364      * range defined by the boundary values of {@code val}.
365      */
366     private int searchIndex(double c, double[] val) {
367         final int r = Arrays.binarySearch(val, c);
368 
369         if (r == -1 ||
370             r == -val.length - 1) {
371             throw new OutOfRangeException(c, val[0], val[val.length - 1]);
372         }
373 
374         if (r < 0) {
375             // "c" in within an interpolation sub-interval: Return the
376             // index of the sample at the lower end of the sub-interval.
377             return -r - 2;
378         }
379         final int last = val.length - 1;
380         if (r == last) {
381             // "c" is the last sample of the range: Return the index
382             // of the sample at the lower end of the last sub-interval.
383             return last - 1;
384         }
385 
386         // "c" is another sample point.
387         return r;
388     }
389 
390     /**
391      * Compute the spline coefficients from the list of function values and
392      * function partial derivatives values at the four corners of a grid
393      * element. They must be specified in the following order:
394      * <ul>
395      *  <li>f(0,0)</li>
396      *  <li>f(1,0)</li>
397      *  <li>f(0,1)</li>
398      *  <li>f(1,1)</li>
399      *  <li>f<sub>x</sub>(0,0)</li>
400      *  <li>f<sub>x</sub>(1,0)</li>
401      *  <li>f<sub>x</sub>(0,1)</li>
402      *  <li>f<sub>x</sub>(1,1)</li>
403      *  <li>f<sub>y</sub>(0,0)</li>
404      *  <li>f<sub>y</sub>(1,0)</li>
405      *  <li>f<sub>y</sub>(0,1)</li>
406      *  <li>f<sub>y</sub>(1,1)</li>
407      *  <li>f<sub>xy</sub>(0,0)</li>
408      *  <li>f<sub>xy</sub>(1,0)</li>
409      *  <li>f<sub>xy</sub>(0,1)</li>
410      *  <li>f<sub>xy</sub>(1,1)</li>
411      * </ul>
412      * where the subscripts indicate the partial derivative with respect to
413      * the corresponding variable(s).
414      *
415      * @param beta List of function values and function partial derivatives
416      * values.
417      * @return the spline coefficients.
418      */
419     private double[] computeSplineCoefficients(double[] beta) {
420         final double[] a = new double[NUM_COEFF];
421 
422         for (int i = 0; i < NUM_COEFF; i++) {
423             double result = 0;
424             final double[] row = AINV[i];
425             for (int j = 0; j < NUM_COEFF; j++) {
426                 result += row[j] * beta[j];
427             }
428             a[i] = result;
429         }
430 
431         return a;
432     }
433 }
434 
435 /**
436  * 2D-spline function.
437  *
438  */
439 class BicubicSplineFunction implements BivariateFunction {
440     /** Number of points. */
441     private static final short N = 4;
442     /** Coefficients */
443     private final double[][] a;
444     /** First partial derivative along x. */
445     private final BivariateFunction partialDerivativeX;
446     /** First partial derivative along y. */
447     private final BivariateFunction partialDerivativeY;
448     /** Second partial derivative along x. */
449     private final BivariateFunction partialDerivativeXX;
450     /** Second partial derivative along y. */
451     private final BivariateFunction partialDerivativeYY;
452     /** Second crossed partial derivative. */
453     private final BivariateFunction partialDerivativeXY;
454 
455     /**
456      * Simple constructor.
457      *
458      * @param coeff Spline coefficients.
459      */
460     public BicubicSplineFunction(double[] coeff) {
461         this(coeff, false);
462     }
463 
464     /**
465      * Simple constructor.
466      *
467      * @param coeff Spline coefficients.
468      * @param initializeDerivatives Whether to initialize the internal data
469      * needed for calling any of the methods that compute the partial derivatives
470      * this function.
471      */
472     public BicubicSplineFunction(double[] coeff,
473                                  boolean initializeDerivatives) {
474         a = new double[N][N];
475         for (int i = 0; i < N; i++) {
476             for (int j = 0; j < N; j++) {
477                 a[i][j] = coeff[i * N + j];
478             }
479         }
480 
481         if (initializeDerivatives) {
482             // Compute all partial derivatives functions.
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 BivariateFunction() {
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 BivariateFunction() {
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 BivariateFunction() {
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 BivariateFunction() {
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 BivariateFunction() {
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         } else {
558             partialDerivativeX = null;
559             partialDerivativeY = null;
560             partialDerivativeXX = null;
561             partialDerivativeYY = null;
562             partialDerivativeXY = null;
563         }
564     }
565 
566     /**
567      * {@inheritDoc}
568      */
569     public double value(double x, double y) {
570         if (x < 0 || x > 1) {
571             throw new OutOfRangeException(x, 0, 1);
572         }
573         if (y < 0 || y > 1) {
574             throw new OutOfRangeException(y, 0, 1);
575         }
576 
577         final double x2 = x * x;
578         final double x3 = x2 * x;
579         final double[] pX = {1, x, x2, x3};
580 
581         final double y2 = y * y;
582         final double y3 = y2 * y;
583         final double[] pY = {1, y, y2, y3};
584 
585         return apply(pX, pY, a);
586     }
587 
588     /**
589      * Compute the value of the bicubic polynomial.
590      *
591      * @param pX Powers of the x-coordinate.
592      * @param pY Powers of the y-coordinate.
593      * @param coeff Spline coefficients.
594      * @return the interpolated value.
595      */
596     private double apply(double[] pX, double[] pY, double[][] coeff) {
597         double result = 0;
598         for (int i = 0; i < N; i++) {
599             for (int j = 0; j < N; j++) {
600                 result += coeff[i][j] * pX[i] * pY[j];
601             }
602         }
603 
604         return result;
605     }
606 
607     /**
608      * @return the partial derivative wrt {@code x}.
609      */
610     public BivariateFunction partialDerivativeX() {
611         return partialDerivativeX;
612     }
613     /**
614      * @return the partial derivative wrt {@code y}.
615      */
616     public BivariateFunction partialDerivativeY() {
617         return partialDerivativeY;
618     }
619     /**
620      * @return the second partial derivative wrt {@code x}.
621      */
622     public BivariateFunction partialDerivativeXX() {
623         return partialDerivativeXX;
624     }
625     /**
626      * @return the second partial derivative wrt {@code y}.
627      */
628     public BivariateFunction partialDerivativeYY() {
629         return partialDerivativeYY;
630     }
631     /**
632      * @return the second partial cross-derivative.
633      */
634     public BivariateFunction partialDerivativeXY() {
635         return partialDerivativeXY;
636     }
637 }