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.math4.legacy.analysis.interpolation;
18
19 import java.util.Arrays;
20 import java.util.function.DoubleBinaryOperator;
21 import java.util.function.Function;
22
23 import org.apache.commons.numbers.core.Sum;
24 import org.apache.commons.math4.legacy.analysis.BivariateFunction;
25 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
26 import org.apache.commons.math4.legacy.exception.NoDataException;
27 import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
28 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
29 import org.apache.commons.math4.legacy.core.MathArrays;
30
31 /**
32 * Function that implements the
33 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
34 * bicubic spline interpolation</a>.
35 *
36 * @since 3.4
37 */
38 public class BicubicInterpolatingFunction
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 BicubicFunction[][] splines;
71
72 /**
73 * @param x Sample values of the x-coordinate, in increasing order.
74 * @param y Sample values of the y-coordinate, in increasing order.
75 * @param f Values of the function on every grid point.
76 * @param dFdX Values of the partial derivative of function with respect
77 * to x on every grid point.
78 * @param dFdY Values of the partial derivative of function with respect
79 * to y on every grid point.
80 * @param d2FdXdY Values of the cross partial derivative of function on
81 * every grid point.
82 * @throws DimensionMismatchException if the various arrays do not contain
83 * the expected number of elements.
84 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
85 * not strictly increasing.
86 * @throws NoDataException if any of the arrays has zero length.
87 */
88 public BicubicInterpolatingFunction(double[] x,
89 double[] y,
90 double[][] f,
91 double[][] dFdX,
92 double[][] dFdY,
93 double[][] d2FdXdY)
94 throws DimensionMismatchException,
95 NoDataException,
96 NonMonotonicSequenceException {
97 this(x, y, f, dFdX, dFdY, d2FdXdY, false);
98 }
99
100 /**
101 * @param x Sample values of the x-coordinate, in increasing order.
102 * @param y Sample values of the y-coordinate, in increasing order.
103 * @param f Values of the function on every grid point.
104 * @param dFdX Values of the partial derivative of function with respect
105 * to x on every grid point.
106 * @param dFdY Values of the partial derivative of function with respect
107 * to y on every grid point.
108 * @param d2FdXdY Values of the cross partial derivative of function on
109 * every grid point.
110 * @param initializeDerivatives Whether to initialize the internal data
111 * needed for calling any of the methods that compute the partial derivatives
112 * this function.
113 * @throws DimensionMismatchException if the various arrays do not contain
114 * the expected number of elements.
115 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
116 * not strictly increasing.
117 * @throws NoDataException if any of the arrays has zero length.
118 */
119 public BicubicInterpolatingFunction(double[] x,
120 double[] y,
121 double[][] f,
122 double[][] dFdX,
123 double[][] dFdY,
124 double[][] d2FdXdY,
125 boolean initializeDerivatives)
126 throws DimensionMismatchException,
127 NoDataException,
128 NonMonotonicSequenceException {
129 final int xLen = x.length;
130 final int yLen = y.length;
131
132 if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
133 throw new NoDataException();
134 }
135 if (xLen != f.length) {
136 throw new DimensionMismatchException(xLen, f.length);
137 }
138 if (xLen != dFdX.length) {
139 throw new DimensionMismatchException(xLen, dFdX.length);
140 }
141 if (xLen != dFdY.length) {
142 throw new DimensionMismatchException(xLen, dFdY.length);
143 }
144 if (xLen != d2FdXdY.length) {
145 throw new DimensionMismatchException(xLen, d2FdXdY.length);
146 }
147
148 MathArrays.checkOrder(x);
149 MathArrays.checkOrder(y);
150
151 xval = x.clone();
152 yval = y.clone();
153
154 final int lastI = xLen - 1;
155 final int lastJ = yLen - 1;
156 splines = new BicubicFunction[lastI][lastJ];
157
158 for (int i = 0; i < lastI; i++) {
159 if (f[i].length != yLen) {
160 throw new DimensionMismatchException(f[i].length, yLen);
161 }
162 if (dFdX[i].length != yLen) {
163 throw new DimensionMismatchException(dFdX[i].length, yLen);
164 }
165 if (dFdY[i].length != yLen) {
166 throw new DimensionMismatchException(dFdY[i].length, yLen);
167 }
168 if (d2FdXdY[i].length != yLen) {
169 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
170 }
171 final int ip1 = i + 1;
172 final double xR = xval[ip1] - xval[i];
173 for (int j = 0; j < lastJ; j++) {
174 final int jp1 = j + 1;
175 final double yR = yval[jp1] - yval[j];
176 final double xRyR = xR * yR;
177 final double[] beta = new double[] {
178 f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
179 dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
180 dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
181 d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
182 };
183
184 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta),
185 xR,
186 yR,
187 initializeDerivatives);
188 }
189 }
190 }
191
192 /**
193 * {@inheritDoc}
194 */
195 @Override
196 public double value(double x, double y)
197 throws OutOfRangeException {
198 final int i = searchIndex(x, xval);
199 final int j = searchIndex(y, yval);
200
201 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
202 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
203
204 return splines[i][j].value(xN, yN);
205 }
206
207 /**
208 * Indicates whether a point is within the interpolation range.
209 *
210 * @param x First coordinate.
211 * @param y Second coordinate.
212 * @return {@code true} if (x, y) is a valid point.
213 */
214 public boolean isValidPoint(double x, double y) {
215 return !(x < xval[0] ||
216 x > xval[xval.length - 1] ||
217 y < yval[0] ||
218 y > yval[yval.length - 1]);
219 }
220
221 /**
222 * @return the first partial derivative respect to x.
223 * @throws NullPointerException if the internal data were not initialized
224 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
225 * double[][],double[][],double[][],boolean) constructor}).
226 */
227 public DoubleBinaryOperator partialDerivativeX() {
228 return partialDerivative(BicubicFunction::partialDerivativeX);
229 }
230
231 /**
232 * @return the first partial derivative respect to y.
233 * @throws NullPointerException if the internal data were not initialized
234 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
235 * double[][],double[][],double[][],boolean) constructor}).
236 */
237 public DoubleBinaryOperator partialDerivativeY() {
238 return partialDerivative(BicubicFunction::partialDerivativeY);
239 }
240
241 /**
242 * @return the second partial derivative respect to x.
243 * @throws NullPointerException if the internal data were not initialized
244 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
245 * double[][],double[][],double[][],boolean) constructor}).
246 */
247 public DoubleBinaryOperator partialDerivativeXX() {
248 return partialDerivative(BicubicFunction::partialDerivativeXX);
249 }
250
251 /**
252 * @return the second partial derivative respect to y.
253 * @throws NullPointerException if the internal data were not initialized
254 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
255 * double[][],double[][],double[][],boolean) constructor}).
256 */
257 public DoubleBinaryOperator partialDerivativeYY() {
258 return partialDerivative(BicubicFunction::partialDerivativeYY);
259 }
260
261 /**
262 * @return the second partial cross derivative.
263 * @throws NullPointerException if the internal data were not initialized
264 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
265 * double[][],double[][],double[][],boolean) constructor}).
266 */
267 public DoubleBinaryOperator partialDerivativeXY() {
268 return partialDerivative(BicubicFunction::partialDerivativeXY);
269 }
270
271 /**
272 * @param which derivative function to apply.
273 * @return the selected partial derivative.
274 * @throws NullPointerException if the internal data were not initialized
275 * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
276 * double[][],double[][],double[][],boolean) constructor}).
277 */
278 private DoubleBinaryOperator partialDerivative(Function<BicubicFunction, BivariateFunction> which) {
279 return (x, y) -> {
280 final int i = searchIndex(x, xval);
281 final int j = searchIndex(y, yval);
282
283 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
284 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
285
286 return which.apply(splines[i][j]).value(xN, yN);
287 };
288 }
289
290 /**
291 * @param c Coordinate.
292 * @param val Coordinate samples.
293 * @return the index in {@code val} corresponding to the interval
294 * containing {@code c}.
295 * @throws OutOfRangeException if {@code c} is out of the
296 * range defined by the boundary values of {@code val}.
297 */
298 private static int searchIndex(double c, double[] val) {
299 final int r = Arrays.binarySearch(val, c);
300
301 if (r == -1 ||
302 r == -val.length - 1) {
303 throw new OutOfRangeException(c, val[0], val[val.length - 1]);
304 }
305
306 if (r < 0) {
307 // "c" in within an interpolation sub-interval: Return the
308 // index of the sample at the lower end of the sub-interval.
309 return -r - 2;
310 }
311 final int last = val.length - 1;
312 if (r == last) {
313 // "c" is the last sample of the range: Return the index
314 // of the sample at the lower end of the last sub-interval.
315 return last - 1;
316 }
317
318 // "c" is another sample point.
319 return r;
320 }
321
322 /**
323 * Compute the spline coefficients from the list of function values and
324 * function partial derivatives values at the four corners of a grid
325 * element. They must be specified in the following order:
326 * <ul>
327 * <li>f(0,0)</li>
328 * <li>f(1,0)</li>
329 * <li>f(0,1)</li>
330 * <li>f(1,1)</li>
331 * <li>f<sub>x</sub>(0,0)</li>
332 * <li>f<sub>x</sub>(1,0)</li>
333 * <li>f<sub>x</sub>(0,1)</li>
334 * <li>f<sub>x</sub>(1,1)</li>
335 * <li>f<sub>y</sub>(0,0)</li>
336 * <li>f<sub>y</sub>(1,0)</li>
337 * <li>f<sub>y</sub>(0,1)</li>
338 * <li>f<sub>y</sub>(1,1)</li>
339 * <li>f<sub>xy</sub>(0,0)</li>
340 * <li>f<sub>xy</sub>(1,0)</li>
341 * <li>f<sub>xy</sub>(0,1)</li>
342 * <li>f<sub>xy</sub>(1,1)</li>
343 * </ul>
344 * where the subscripts indicate the partial derivative with respect to
345 * the corresponding variable(s).
346 *
347 * @param beta List of function values and function partial derivatives
348 * values.
349 * @return the spline coefficients.
350 */
351 private static double[] computeSplineCoefficients(double[] beta) {
352 final double[] a = new double[NUM_COEFF];
353
354 for (int i = 0; i < NUM_COEFF; i++) {
355 double result = 0;
356 final double[] row = AINV[i];
357 for (int j = 0; j < NUM_COEFF; j++) {
358 result += row[j] * beta[j];
359 }
360 a[i] = result;
361 }
362
363 return a;
364 }
365
366 /**
367 * Bicubic function.
368 */
369 static final class BicubicFunction implements BivariateFunction {
370 /** Number of points. */
371 private static final short N = 4;
372 /** Coefficients. */
373 private final double[][] a;
374 /** First partial derivative along x. */
375 private final BivariateFunction partialDerivativeX;
376 /** First partial derivative along y. */
377 private final BivariateFunction partialDerivativeY;
378 /** Second partial derivative along x. */
379 private final BivariateFunction partialDerivativeXX;
380 /** Second partial derivative along y. */
381 private final BivariateFunction partialDerivativeYY;
382 /** Second crossed partial derivative. */
383 private final BivariateFunction partialDerivativeXY;
384
385 /**
386 * Simple constructor.
387 *
388 * @param coeff Spline coefficients.
389 * @param xR x spacing.
390 * @param yR y spacing.
391 * @param initializeDerivatives Whether to initialize the internal data
392 * needed for calling any of the methods that compute the partial derivatives
393 * this function.
394 */
395 BicubicFunction(double[] coeff,
396 double xR,
397 double yR,
398 boolean initializeDerivatives) {
399 a = new double[N][N];
400 for (int j = 0; j < N; j++) {
401 final double[] aJ = a[j];
402 for (int i = 0; i < N; i++) {
403 aJ[i] = coeff[i * N + j];
404 }
405 }
406
407 if (initializeDerivatives) {
408 // Compute all partial derivatives functions.
409 final double[][] aX = new double[N][N];
410 final double[][] aY = new double[N][N];
411 final double[][] aXX = new double[N][N];
412 final double[][] aYY = new double[N][N];
413 final double[][] aXY = new double[N][N];
414
415 for (int i = 0; i < N; i++) {
416 for (int j = 0; j < N; j++) {
417 final double c = a[i][j];
418 aX[i][j] = i * c;
419 aY[i][j] = j * c;
420 aXX[i][j] = (i - 1) * aX[i][j];
421 aYY[i][j] = (j - 1) * aY[i][j];
422 aXY[i][j] = j * aX[i][j];
423 }
424 }
425
426 partialDerivativeX = (double x, double y) -> {
427 final double x2 = x * x;
428 final double[] pX = {0, 1, x, x2};
429
430 final double y2 = y * y;
431 final double y3 = y2 * y;
432 final double[] pY = {1, y, y2, y3};
433
434 return apply(pX, 1, pY, 0, aX) / xR;
435 };
436 partialDerivativeY = (double x, double y) -> {
437 final double x2 = x * x;
438 final double x3 = x2 * x;
439 final double[] pX = {1, x, x2, x3};
440
441 final double y2 = y * y;
442 final double[] pY = {0, 1, y, y2};
443
444 return apply(pX, 0, pY, 1, aY) / yR;
445 };
446 partialDerivativeXX = (double x, double y) -> {
447 final double[] pX = {0, 0, 1, x};
448
449 final double y2 = y * y;
450 final double y3 = y2 * y;
451 final double[] pY = {1, y, y2, y3};
452
453 return apply(pX, 2, pY, 0, aXX) / (xR * xR);
454 };
455 partialDerivativeYY = (double x, double y) -> {
456 final double x2 = x * x;
457 final double x3 = x2 * x;
458 final double[] pX = {1, x, x2, x3};
459
460 final double[] pY = {0, 0, 1, y};
461
462 return apply(pX, 0, pY, 2, aYY) / (yR * yR);
463 };
464 partialDerivativeXY = (double x, double y) -> {
465 final double x2 = x * x;
466 final double[] pX = {0, 1, x, x2};
467
468 final double y2 = y * y;
469 final double[] pY = {0, 1, y, y2};
470
471 return apply(pX, 1, pY, 1, aXY) / (xR * yR);
472 };
473 } else {
474 partialDerivativeX = null;
475 partialDerivativeY = null;
476 partialDerivativeXX = null;
477 partialDerivativeYY = null;
478 partialDerivativeXY = null;
479 }
480 }
481
482 /**
483 * {@inheritDoc}
484 */
485 @Override
486 public double value(double x, double y) {
487 if (x < 0 || x > 1) {
488 throw new OutOfRangeException(x, 0, 1);
489 }
490 if (y < 0 || y > 1) {
491 throw new OutOfRangeException(y, 0, 1);
492 }
493
494 final double x2 = x * x;
495 final double x3 = x2 * x;
496 final double[] pX = {1, x, x2, x3};
497
498 final double y2 = y * y;
499 final double y3 = y2 * y;
500 final double[] pY = {1, y, y2, y3};
501
502 return apply(pX, 0, pY, 0, a);
503 }
504
505 /**
506 * Compute the value of the bicubic polynomial.
507 *
508 * <p>Assumes the powers are zero below the provided index, and 1 at the provided
509 * index. This allows skipping some zero products and optimising multiplication
510 * by one.
511 *
512 * @param pX Powers of the x-coordinate.
513 * @param i Index of pX[i] == 1
514 * @param pY Powers of the y-coordinate.
515 * @param j Index of pX[j] == 1
516 * @param coeff Spline coefficients.
517 * @return the interpolated value.
518 */
519 private static double apply(double[] pX, int i, double[] pY, int j, double[][] coeff) {
520 // assert pX[i] == 1
521 double result = sumOfProducts(coeff[i], pY, j);
522 while (++i < N) {
523 final double r = sumOfProducts(coeff[i], pY, j);
524 result += r * pX[i];
525 }
526 return result;
527 }
528
529 /**
530 * Compute the sum of products starting from the provided index.
531 * Assumes that factor {@code b[j] == 1}.
532 *
533 * @param a Factors.
534 * @param b Factors.
535 * @param j Index to initialise the sum.
536 * @return the double
537 */
538 private static double sumOfProducts(double[] a, double[] b, int j) {
539 // assert b[j] == 1
540 final Sum sum = Sum.of(a[j]);
541 while (++j < N) {
542 sum.addProduct(a[j], b[j]);
543 }
544 return sum.getAsDouble();
545 }
546
547 /**
548 * @return the partial derivative wrt {@code x}.
549 */
550 BivariateFunction partialDerivativeX() {
551 return partialDerivativeX;
552 }
553
554 /**
555 * @return the partial derivative wrt {@code y}.
556 */
557 BivariateFunction partialDerivativeY() {
558 return partialDerivativeY;
559 }
560
561 /**
562 * @return the second partial derivative wrt {@code x}.
563 */
564 BivariateFunction partialDerivativeXX() {
565 return partialDerivativeXX;
566 }
567
568 /**
569 * @return the second partial derivative wrt {@code y}.
570 */
571 BivariateFunction partialDerivativeYY() {
572 return partialDerivativeYY;
573 }
574
575 /**
576 * @return the second partial cross-derivative.
577 */
578 BivariateFunction partialDerivativeXY() {
579 return partialDerivativeXY;
580 }
581 }
582 }
583