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