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>.
31   *
32   * @since 2.1
33   * @version $Id: BicubicSplineInterpolatingFunction.java 1512547 2013-08-10 01:13:38Z erans $
34   */
35  public class BicubicSplineInterpolatingFunction
36      implements BivariateFunction {
37      /** Number of coefficients. */
38      private static final int NUM_COEFF = 16;
39      /**
40       * Matrix to compute the spline coefficients from the function values
41       * and function derivatives values
42       */
43      private static final double[][] AINV = {
44          { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
45          { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
46          { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
47          { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
48          { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
49          { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
50          { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
51          { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
52          { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
53          { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
54          { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
55          { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
56          { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
57          { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
58          { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
59          { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
60      };
61  
62      /** Samples x-coordinates */
63      private final double[] xval;
64      /** Samples y-coordinates */
65      private final double[] yval;
66      /** Set of cubic splines patching the whole data grid */
67      private final BicubicSplineFunction[][] splines;
68      /**
69       * Partial derivatives
70       * The value of the first index determines the kind of derivatives:
71       * 0 = first partial derivatives wrt x
72       * 1 = first partial derivatives wrt y
73       * 2 = second partial derivatives wrt x
74       * 3 = second partial derivatives wrt y
75       * 4 = cross partial derivatives
76       */
77      private BivariateFunction[][][] partialDerivatives = null;
78  
79      /**
80       * @param x Sample values of the x-coordinate, in increasing order.
81       * @param y Sample values of the y-coordinate, in increasing order.
82       * @param f Values of the function on every grid point.
83       * @param dFdX Values of the partial derivative of function with respect
84       * to x on every grid point.
85       * @param dFdY Values of the partial derivative of function with respect
86       * to y on every grid point.
87       * @param d2FdXdY Values of the cross partial derivative of function on
88       * every grid point.
89       * @throws DimensionMismatchException if the various arrays do not contain
90       * the expected number of elements.
91       * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
92       * not strictly increasing.
93       * @throws NoDataException if any of the arrays has zero length.
94       */
95      public BicubicSplineInterpolatingFunction(double[] x,
96                                                double[] y,
97                                                double[][] f,
98                                                double[][] dFdX,
99                                                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  */
385 class 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 }