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.math.analysis.interpolation; 18 19 import org.apache.commons.math.DimensionMismatchException; 20 import org.apache.commons.math.FunctionEvaluationException; 21 import org.apache.commons.math.analysis.BivariateRealFunction; 22 import org.apache.commons.math.exception.NoDataException; 23 import org.apache.commons.math.exception.OutOfRangeException; 24 import org.apache.commons.math.util.MathUtils; 25 26 /** 27 * Function that implements the 28 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation"> 29 * bicubic spline interpolation</a>. 30 * 31 * @version $Revision$ $Date$ 32 * @since 2.1 33 */ 34 public class BicubicSplineInterpolatingFunction 35 implements BivariateRealFunction { 36 /** 37 * Matrix to compute the spline coefficients from the function values 38 * and function derivatives values 39 */ 40 private static final double[][] AINV = { 41 { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, 42 { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, 43 { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, 44 { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, 45 { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, 46 { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, 47 { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, 48 { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, 49 { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, 50 { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, 51 { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, 52 { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, 53 { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, 54 { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, 55 { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, 56 { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } 57 }; 58 59 /** Samples x-coordinates */ 60 private final double[] xval; 61 /** Samples y-coordinates */ 62 private final double[] yval; 63 /** Set of cubic splines patching the whole data grid */ 64 private final BicubicSplineFunction[][] splines; 65 /** 66 * Partial derivatives 67 * The value of the first index determines the kind of derivatives: 68 * 0 = first partial derivatives wrt x 69 * 1 = first partial derivatives wrt y 70 * 2 = second partial derivatives wrt x 71 * 3 = second partial derivatives wrt y 72 * 4 = cross partial derivatives 73 */ 74 private BivariateRealFunction[][][] partialDerivatives = null; 75 76 /** 77 * @param x Sample values of the x-coordinate, in increasing order. 78 * @param y Sample values of the y-coordinate, in increasing order. 79 * @param f Values of the function on every grid point. 80 * @param dFdX Values of the partial derivative of function with respect 81 * to x on every grid point. 82 * @param dFdY Values of the partial derivative of function with respect 83 * to y on every grid point. 84 * @param d2FdXdY Values of the cross partial derivative of function on 85 * every grid point. 86 * @throws DimensionMismatchException if the various arrays do not contain 87 * the expected number of elements. 88 * @throws org.apache.commons.math.exception.NonMonotonousSequenceException 89 * if {@code x} or {@code y} are not strictly increasing. 90 * @throws NoDataException if any of the arrays has zero length. 91 */ BicubicSplineInterpolatingFunction(double[] x, double[] y, double[][] f, double[][] dFdX, double[][] dFdY, double[][] d2FdXdY)92 public BicubicSplineInterpolatingFunction(double[] x, 93 double[] y, 94 double[][] f, 95 double[][] dFdX, 96 double[][] dFdY, 97 double[][] d2FdXdY) 98 throws DimensionMismatchException { 99 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 */ value(double x, double y)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 */ partialDerivativeX(double x, double y)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 */ partialDerivativeY(double x, double y)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 */ partialDerivativeXX(double x, double y)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 */ partialDerivativeYY(double x, double y)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 */ partialDerivativeXY(double x, double y)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 */ partialDerivative(int which, double x, double y)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 */ computePartialDerivatives()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 */ searchIndex(double c, double[] val)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 */ computeSplineCoefficients(double[] beta)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 */ BicubicSplineFunction(double[] a)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 */ value(double x, double y)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 */ apply(double[] pX, double[] pY, double[][] coeff)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 */ partialDerivativeX()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 */ partialDerivativeY()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 */ partialDerivativeXX()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 */ partialDerivativeYY()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 */ partialDerivativeXY()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 */ computePartialDerivatives()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 } 559