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.math.analysis.interpolation;
018
019 import org.apache.commons.math.DimensionMismatchException;
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.analysis.BivariateRealFunction;
022 import org.apache.commons.math.exception.NoDataException;
023 import org.apache.commons.math.exception.OutOfRangeException;
024 import org.apache.commons.math.util.MathUtils;
025
026 /**
027 * Function that implements the
028 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
029 * bicubic spline interpolation</a>.
030 *
031 * @version $Revision$ $Date$
032 * @since 2.1
033 */
034 public class BicubicSplineInterpolatingFunction
035 implements BivariateRealFunction {
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 BivariateRealFunction[][][] 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 org.apache.commons.math.exception.NonMonotonousSequenceException
089 * if {@code x} or {@code y} are 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 final int xLen = x.length;
100 final int yLen = y.length;
101
102 if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
103 throw new NoDataException();
104 }
105 if (xLen != f.length) {
106 throw new DimensionMismatchException(xLen, f.length);
107 }
108 if (xLen != dFdX.length) {
109 throw new DimensionMismatchException(xLen, dFdX.length);
110 }
111 if (xLen != dFdY.length) {
112 throw new DimensionMismatchException(xLen, dFdY.length);
113 }
114 if (xLen != d2FdXdY.length) {
115 throw new DimensionMismatchException(xLen, d2FdXdY.length);
116 }
117
118 MathUtils.checkOrder(x);
119 MathUtils.checkOrder(y);
120
121 xval = x.clone();
122 yval = y.clone();
123
124 final int lastI = xLen - 1;
125 final int lastJ = yLen - 1;
126 splines = new BicubicSplineFunction[lastI][lastJ];
127
128 for (int i = 0; i < lastI; i++) {
129 if (f[i].length != yLen) {
130 throw new DimensionMismatchException(f[i].length, yLen);
131 }
132 if (dFdX[i].length != yLen) {
133 throw new DimensionMismatchException(dFdX[i].length, yLen);
134 }
135 if (dFdY[i].length != yLen) {
136 throw new DimensionMismatchException(dFdY[i].length, yLen);
137 }
138 if (d2FdXdY[i].length != yLen) {
139 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
140 }
141 final int ip1 = i + 1;
142 for (int j = 0; j < lastJ; j++) {
143 final int jp1 = j + 1;
144 final double[] beta = new double[] {
145 f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
146 dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
147 dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
148 d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
149 };
150
151 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
152 }
153 }
154 }
155
156 /**
157 * {@inheritDoc}
158 */
159 public double value(double x, double y) {
160 final int i = searchIndex(x, xval);
161 if (i == -1) {
162 throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
163 }
164 final int j = searchIndex(y, yval);
165 if (j == -1) {
166 throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
167 }
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 * @param x x-coordinate.
177 * @param y y-coordinate.
178 * @return the value at point (x, y) of the first partial derivative with
179 * respect to x.
180 * @since 2.2
181 */
182 public double partialDerivativeX(double x, double y) {
183 return partialDerivative(0, x, y);
184 }
185 /**
186 * @param x x-coordinate.
187 * @param y y-coordinate.
188 * @return the value at point (x, y) of the first partial derivative with
189 * respect to y.
190 * @since 2.2
191 */
192 public double partialDerivativeY(double x, double y) {
193 return partialDerivative(1, x, y);
194 }
195 /**
196 * @param x x-coordinate.
197 * @param y y-coordinate.
198 * @return the value at point (x, y) of the second partial derivative with
199 * respect to x.
200 * @since 2.2
201 */
202 public double partialDerivativeXX(double x, double y) {
203 return partialDerivative(2, x, y);
204 }
205 /**
206 * @param x x-coordinate.
207 * @param y y-coordinate.
208 * @return the value at point (x, y) of the second partial derivative with
209 * respect to y.
210 * @since 2.2
211 */
212 public double partialDerivativeYY(double x, double y) {
213 return partialDerivative(3, x, y);
214 }
215 /**
216 * @param x x-coordinate.
217 * @param y y-coordinate.
218 * @return the value at point (x, y) of the second partial cross-derivative.
219 * @since 2.2
220 */
221 public double partialDerivativeXY(double x, double y) {
222 return partialDerivative(4, x, y);
223 }
224
225 /**
226 * @param which First index in {@link #partialDerivatives}.
227 * @param x x-coordinate.
228 * @param y y-coordinate.
229 * @return the value at point (x, y) of the selected partial derivative.
230 * @throws FunctionEvaluationException
231 */
232 private double partialDerivative(int which, double x, double y) {
233 if (partialDerivatives == null) {
234 computePartialDerivatives();
235 }
236
237 final int i = searchIndex(x, xval);
238 if (i == -1) {
239 throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
240 }
241 final int j = searchIndex(y, yval);
242 if (j == -1) {
243 throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
244 }
245
246 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
247 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
248
249 try {
250 return partialDerivatives[which][i][j].value(xN, yN);
251 } catch (FunctionEvaluationException fee) {
252 // this should never happen
253 throw new RuntimeException(fee);
254 }
255
256 }
257
258 /**
259 * Compute all partial derivatives.
260 */
261 private void computePartialDerivatives() {
262 final int lastI = xval.length - 1;
263 final int lastJ = yval.length - 1;
264 partialDerivatives = new BivariateRealFunction[5][lastI][lastJ];
265
266 for (int i = 0; i < lastI; i++) {
267 for (int j = 0; j < lastJ; j++) {
268 final BicubicSplineFunction f = splines[i][j];
269 partialDerivatives[0][i][j] = f.partialDerivativeX();
270 partialDerivatives[1][i][j] = f.partialDerivativeY();
271 partialDerivatives[2][i][j] = f.partialDerivativeXX();
272 partialDerivatives[3][i][j] = f.partialDerivativeYY();
273 partialDerivatives[4][i][j] = f.partialDerivativeXY();
274 }
275 }
276 }
277
278 /**
279 * @param c Coordinate.
280 * @param val Coordinate samples.
281 * @return the index in {@code val} corresponding to the interval
282 * containing {@code c}, or {@code -1} if {@code c} is out of the
283 * range defined by the end values of {@code val}.
284 */
285 private int searchIndex(double c, double[] val) {
286 if (c < val[0]) {
287 return -1;
288 }
289
290 final int max = val.length;
291 for (int i = 1; i < max; i++) {
292 if (c <= val[i]) {
293 return i - 1;
294 }
295 }
296
297 return -1;
298 }
299
300 /**
301 * Compute the spline coefficients from the list of function values and
302 * function partial derivatives values at the four corners of a grid
303 * element. They must be specified in the following order:
304 * <ul>
305 * <li>f(0,0)</li>
306 * <li>f(1,0)</li>
307 * <li>f(0,1)</li>
308 * <li>f(1,1)</li>
309 * <li>f<sub>x</sub>(0,0)</li>
310 * <li>f<sub>x</sub>(1,0)</li>
311 * <li>f<sub>x</sub>(0,1)</li>
312 * <li>f<sub>x</sub>(1,1)</li>
313 * <li>f<sub>y</sub>(0,0)</li>
314 * <li>f<sub>y</sub>(1,0)</li>
315 * <li>f<sub>y</sub>(0,1)</li>
316 * <li>f<sub>y</sub>(1,1)</li>
317 * <li>f<sub>xy</sub>(0,0)</li>
318 * <li>f<sub>xy</sub>(1,0)</li>
319 * <li>f<sub>xy</sub>(0,1)</li>
320 * <li>f<sub>xy</sub>(1,1)</li>
321 * </ul>
322 * where the subscripts indicate the partial derivative with respect to
323 * the corresponding variable(s).
324 *
325 * @param beta List of function values and function partial derivatives
326 * values.
327 * @return the spline coefficients.
328 */
329 private double[] computeSplineCoefficients(double[] beta) {
330 final double[] a = new double[16];
331
332 for (int i = 0; i < 16; i++) {
333 double result = 0;
334 final double[] row = AINV[i];
335 for (int j = 0; j < 16; j++) {
336 result += row[j] * beta[j];
337 }
338 a[i] = result;
339 }
340
341 return a;
342 }
343 }
344
345 /**
346 * 2D-spline function.
347 *
348 * @version $Revision$ $Date$
349 */
350 class BicubicSplineFunction
351 implements BivariateRealFunction {
352
353 /** Number of points. */
354 private static final short N = 4;
355
356 /** Coefficients */
357 private final double[][] a;
358
359 /** First partial derivative along x. */
360 private BivariateRealFunction partialDerivativeX;
361
362 /** First partial derivative along y. */
363 private BivariateRealFunction partialDerivativeY;
364
365 /** Second partial derivative along x. */
366 private BivariateRealFunction partialDerivativeXX;
367
368 /** Second partial derivative along y. */
369 private BivariateRealFunction partialDerivativeYY;
370
371 /** Second crossed partial derivative. */
372 private BivariateRealFunction partialDerivativeXY;
373
374 /**
375 * Simple constructor.
376 * @param a Spline coefficients
377 */
378 public BicubicSplineFunction(double[] a) {
379 this.a = new double[N][N];
380 for (int i = 0; i < N; i++) {
381 for (int j = 0; j < N; j++) {
382 this.a[i][j] = a[i + N * j];
383 }
384 }
385 }
386
387 /**
388 * {@inheritDoc}
389 */
390 public double value(double x, double y) {
391 if (x < 0 || x > 1) {
392 throw new OutOfRangeException(x, 0, 1);
393 }
394 if (y < 0 || y > 1) {
395 throw new OutOfRangeException(y, 0, 1);
396 }
397
398 final double x2 = x * x;
399 final double x3 = x2 * x;
400 final double[] pX = {1, x, x2, x3};
401
402 final double y2 = y * y;
403 final double y3 = y2 * y;
404 final double[] pY = {1, y, y2, y3};
405
406 return apply(pX, pY, a);
407 }
408
409 /**
410 * Compute the value of the bicubic polynomial.
411 *
412 * @param pX Powers of the x-coordinate.
413 * @param pY Powers of the y-coordinate.
414 * @param coeff Spline coefficients.
415 * @return the interpolated value.
416 */
417 private double apply(double[] pX, double[] pY, double[][] coeff) {
418 double result = 0;
419 for (int i = 0; i < N; i++) {
420 for (int j = 0; j < N; j++) {
421 result += coeff[i][j] * pX[i] * pY[j];
422 }
423 }
424
425 return result;
426 }
427
428 /**
429 * @return the partial derivative wrt {@code x}.
430 */
431 public BivariateRealFunction partialDerivativeX() {
432 if (partialDerivativeX == null) {
433 computePartialDerivatives();
434 }
435
436 return partialDerivativeX;
437 }
438 /**
439 * @return the partial derivative wrt {@code y}.
440 */
441 public BivariateRealFunction partialDerivativeY() {
442 if (partialDerivativeY == null) {
443 computePartialDerivatives();
444 }
445
446 return partialDerivativeY;
447 }
448 /**
449 * @return the second partial derivative wrt {@code x}.
450 */
451 public BivariateRealFunction partialDerivativeXX() {
452 if (partialDerivativeXX == null) {
453 computePartialDerivatives();
454 }
455
456 return partialDerivativeXX;
457 }
458 /**
459 * @return the second partial derivative wrt {@code y}.
460 */
461 public BivariateRealFunction partialDerivativeYY() {
462 if (partialDerivativeYY == null) {
463 computePartialDerivatives();
464 }
465
466 return partialDerivativeYY;
467 }
468 /**
469 * @return the second partial cross-derivative.
470 */
471 public BivariateRealFunction partialDerivativeXY() {
472 if (partialDerivativeXY == null) {
473 computePartialDerivatives();
474 }
475
476 return partialDerivativeXY;
477 }
478
479 /**
480 * Compute all partial derivatives functions.
481 */
482 private void computePartialDerivatives() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 }
558 }