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 }