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