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