BicubicInterpolatingFunction.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.legacy.analysis.interpolation;
import java.util.Arrays;
import java.util.function.DoubleBinaryOperator;
import java.util.function.Function;
import org.apache.commons.numbers.core.Sum;
import org.apache.commons.math4.legacy.analysis.BivariateFunction;
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
import org.apache.commons.math4.legacy.exception.NoDataException;
import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.legacy.exception.OutOfRangeException;
import org.apache.commons.math4.legacy.core.MathArrays;
/**
* Function that implements the
* <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
* bicubic spline interpolation</a>.
*
* @since 3.4
*/
public class BicubicInterpolatingFunction
implements BivariateFunction {
/** Number of coefficients. */
private static final int NUM_COEFF = 16;
/**
* Matrix to compute the spline coefficients from the function values
* and function derivatives values.
*/
private static final double[][] AINV = {
{ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
{ -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
{ 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
{ 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
{ 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
{ -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
{ 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
{ 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
{ -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
{ 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
{ 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
{ -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
{ 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
};
/** Samples x-coordinates. */
private final double[] xval;
/** Samples y-coordinates. */
private final double[] yval;
/** Set of cubic splines patching the whole data grid. */
private final BicubicFunction[][] splines;
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param f Values of the function on every grid point.
* @param dFdX Values of the partial derivative of function with respect
* to x on every grid point.
* @param dFdY Values of the partial derivative of function with respect
* to y on every grid point.
* @param d2FdXdY Values of the cross partial derivative of function on
* every grid point.
* @throws DimensionMismatchException if the various arrays do not contain
* the expected number of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are
* not strictly increasing.
* @throws NoDataException if any of the arrays has zero length.
*/
public BicubicInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
this(x, y, f, dFdX, dFdY, d2FdXdY, false);
}
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param f Values of the function on every grid point.
* @param dFdX Values of the partial derivative of function with respect
* to x on every grid point.
* @param dFdY Values of the partial derivative of function with respect
* to y on every grid point.
* @param d2FdXdY Values of the cross partial derivative of function on
* every grid point.
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* this function.
* @throws DimensionMismatchException if the various arrays do not contain
* the expected number of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are
* not strictly increasing.
* @throws NoDataException if any of the arrays has zero length.
*/
public BicubicInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY,
boolean initializeDerivatives)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
final int xLen = x.length;
final int yLen = y.length;
if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
throw new NoDataException();
}
if (xLen != f.length) {
throw new DimensionMismatchException(xLen, f.length);
}
if (xLen != dFdX.length) {
throw new DimensionMismatchException(xLen, dFdX.length);
}
if (xLen != dFdY.length) {
throw new DimensionMismatchException(xLen, dFdY.length);
}
if (xLen != d2FdXdY.length) {
throw new DimensionMismatchException(xLen, d2FdXdY.length);
}
MathArrays.checkOrder(x);
MathArrays.checkOrder(y);
xval = x.clone();
yval = y.clone();
final int lastI = xLen - 1;
final int lastJ = yLen - 1;
splines = new BicubicFunction[lastI][lastJ];
for (int i = 0; i < lastI; i++) {
if (f[i].length != yLen) {
throw new DimensionMismatchException(f[i].length, yLen);
}
if (dFdX[i].length != yLen) {
throw new DimensionMismatchException(dFdX[i].length, yLen);
}
if (dFdY[i].length != yLen) {
throw new DimensionMismatchException(dFdY[i].length, yLen);
}
if (d2FdXdY[i].length != yLen) {
throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
}
final int ip1 = i + 1;
final double xR = xval[ip1] - xval[i];
for (int j = 0; j < lastJ; j++) {
final int jp1 = j + 1;
final double yR = yval[jp1] - yval[j];
final double xRyR = xR * yR;
final double[] beta = new double[] {
f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
};
splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta),
xR,
yR,
initializeDerivatives);
}
}
}
/**
* {@inheritDoc}
*/
@Override
public double value(double x, double y)
throws OutOfRangeException {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
return splines[i][j].value(xN, yN);
}
/**
* Indicates whether a point is within the interpolation range.
*
* @param x First coordinate.
* @param y Second coordinate.
* @return {@code true} if (x, y) is a valid point.
*/
public boolean isValidPoint(double x, double y) {
return !(x < xval[0] ||
x > xval[xval.length - 1] ||
y < yval[0] ||
y > yval[yval.length - 1]);
}
/**
* @return the first partial derivative respect to x.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public DoubleBinaryOperator partialDerivativeX() {
return partialDerivative(BicubicFunction::partialDerivativeX);
}
/**
* @return the first partial derivative respect to y.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public DoubleBinaryOperator partialDerivativeY() {
return partialDerivative(BicubicFunction::partialDerivativeY);
}
/**
* @return the second partial derivative respect to x.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public DoubleBinaryOperator partialDerivativeXX() {
return partialDerivative(BicubicFunction::partialDerivativeXX);
}
/**
* @return the second partial derivative respect to y.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public DoubleBinaryOperator partialDerivativeYY() {
return partialDerivative(BicubicFunction::partialDerivativeYY);
}
/**
* @return the second partial cross derivative.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public DoubleBinaryOperator partialDerivativeXY() {
return partialDerivative(BicubicFunction::partialDerivativeXY);
}
/**
* @param which derivative function to apply.
* @return the selected partial derivative.
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
private DoubleBinaryOperator partialDerivative(Function<BicubicFunction, BivariateFunction> which) {
return (x, y) -> {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
return which.apply(splines[i][j]).value(xN, yN);
};
}
/**
* @param c Coordinate.
* @param val Coordinate samples.
* @return the index in {@code val} corresponding to the interval
* containing {@code c}.
* @throws OutOfRangeException if {@code c} is out of the
* range defined by the boundary values of {@code val}.
*/
private static int searchIndex(double c, double[] val) {
final int r = Arrays.binarySearch(val, c);
if (r == -1 ||
r == -val.length - 1) {
throw new OutOfRangeException(c, val[0], val[val.length - 1]);
}
if (r < 0) {
// "c" in within an interpolation sub-interval: Return the
// index of the sample at the lower end of the sub-interval.
return -r - 2;
}
final int last = val.length - 1;
if (r == last) {
// "c" is the last sample of the range: Return the index
// of the sample at the lower end of the last sub-interval.
return last - 1;
}
// "c" is another sample point.
return r;
}
/**
* Compute the spline coefficients from the list of function values and
* function partial derivatives values at the four corners of a grid
* element. They must be specified in the following order:
* <ul>
* <li>f(0,0)</li>
* <li>f(1,0)</li>
* <li>f(0,1)</li>
* <li>f(1,1)</li>
* <li>f<sub>x</sub>(0,0)</li>
* <li>f<sub>x</sub>(1,0)</li>
* <li>f<sub>x</sub>(0,1)</li>
* <li>f<sub>x</sub>(1,1)</li>
* <li>f<sub>y</sub>(0,0)</li>
* <li>f<sub>y</sub>(1,0)</li>
* <li>f<sub>y</sub>(0,1)</li>
* <li>f<sub>y</sub>(1,1)</li>
* <li>f<sub>xy</sub>(0,0)</li>
* <li>f<sub>xy</sub>(1,0)</li>
* <li>f<sub>xy</sub>(0,1)</li>
* <li>f<sub>xy</sub>(1,1)</li>
* </ul>
* where the subscripts indicate the partial derivative with respect to
* the corresponding variable(s).
*
* @param beta List of function values and function partial derivatives
* values.
* @return the spline coefficients.
*/
private static double[] computeSplineCoefficients(double[] beta) {
final double[] a = new double[NUM_COEFF];
for (int i = 0; i < NUM_COEFF; i++) {
double result = 0;
final double[] row = AINV[i];
for (int j = 0; j < NUM_COEFF; j++) {
result += row[j] * beta[j];
}
a[i] = result;
}
return a;
}
/**
* Bicubic function.
*/
static final class BicubicFunction implements BivariateFunction {
/** Number of points. */
private static final short N = 4;
/** Coefficients. */
private final double[][] a;
/** First partial derivative along x. */
private final BivariateFunction partialDerivativeX;
/** First partial derivative along y. */
private final BivariateFunction partialDerivativeY;
/** Second partial derivative along x. */
private final BivariateFunction partialDerivativeXX;
/** Second partial derivative along y. */
private final BivariateFunction partialDerivativeYY;
/** Second crossed partial derivative. */
private final BivariateFunction partialDerivativeXY;
/**
* Simple constructor.
*
* @param coeff Spline coefficients.
* @param xR x spacing.
* @param yR y spacing.
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* this function.
*/
BicubicFunction(double[] coeff,
double xR,
double yR,
boolean initializeDerivatives) {
a = new double[N][N];
for (int j = 0; j < N; j++) {
final double[] aJ = a[j];
for (int i = 0; i < N; i++) {
aJ[i] = coeff[i * N + j];
}
}
if (initializeDerivatives) {
// Compute all partial derivatives functions.
final double[][] aX = new double[N][N];
final double[][] aY = new double[N][N];
final double[][] aXX = new double[N][N];
final double[][] aYY = new double[N][N];
final double[][] aXY = new double[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
final double c = a[i][j];
aX[i][j] = i * c;
aY[i][j] = j * c;
aXX[i][j] = (i - 1) * aX[i][j];
aYY[i][j] = (j - 1) * aY[i][j];
aXY[i][j] = j * aX[i][j];
}
}
partialDerivativeX = (double x, double y) -> {
final double x2 = x * x;
final double[] pX = {0, 1, x, x2};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, 1, pY, 0, aX) / xR;
};
partialDerivativeY = (double x, double y) -> {
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double y2 = y * y;
final double[] pY = {0, 1, y, y2};
return apply(pX, 0, pY, 1, aY) / yR;
};
partialDerivativeXX = (double x, double y) -> {
final double[] pX = {0, 0, 1, x};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, 2, pY, 0, aXX) / (xR * xR);
};
partialDerivativeYY = (double x, double y) -> {
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double[] pY = {0, 0, 1, y};
return apply(pX, 0, pY, 2, aYY) / (yR * yR);
};
partialDerivativeXY = (double x, double y) -> {
final double x2 = x * x;
final double[] pX = {0, 1, x, x2};
final double y2 = y * y;
final double[] pY = {0, 1, y, y2};
return apply(pX, 1, pY, 1, aXY) / (xR * yR);
};
} else {
partialDerivativeX = null;
partialDerivativeY = null;
partialDerivativeXX = null;
partialDerivativeYY = null;
partialDerivativeXY = null;
}
}
/**
* {@inheritDoc}
*/
@Override
public double value(double x, double y) {
if (x < 0 || x > 1) {
throw new OutOfRangeException(x, 0, 1);
}
if (y < 0 || y > 1) {
throw new OutOfRangeException(y, 0, 1);
}
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, 0, pY, 0, a);
}
/**
* Compute the value of the bicubic polynomial.
*
* <p>Assumes the powers are zero below the provided index, and 1 at the provided
* index. This allows skipping some zero products and optimising multiplication
* by one.
*
* @param pX Powers of the x-coordinate.
* @param i Index of pX[i] == 1
* @param pY Powers of the y-coordinate.
* @param j Index of pX[j] == 1
* @param coeff Spline coefficients.
* @return the interpolated value.
*/
private static double apply(double[] pX, int i, double[] pY, int j, double[][] coeff) {
// assert pX[i] == 1
double result = sumOfProducts(coeff[i], pY, j);
while (++i < N) {
final double r = sumOfProducts(coeff[i], pY, j);
result += r * pX[i];
}
return result;
}
/**
* Compute the sum of products starting from the provided index.
* Assumes that factor {@code b[j] == 1}.
*
* @param a Factors.
* @param b Factors.
* @param j Index to initialise the sum.
* @return the double
*/
private static double sumOfProducts(double[] a, double[] b, int j) {
// assert b[j] == 1
final Sum sum = Sum.of(a[j]);
while (++j < N) {
sum.addProduct(a[j], b[j]);
}
return sum.getAsDouble();
}
/**
* @return the partial derivative wrt {@code x}.
*/
BivariateFunction partialDerivativeX() {
return partialDerivativeX;
}
/**
* @return the partial derivative wrt {@code y}.
*/
BivariateFunction partialDerivativeY() {
return partialDerivativeY;
}
/**
* @return the second partial derivative wrt {@code x}.
*/
BivariateFunction partialDerivativeXX() {
return partialDerivativeXX;
}
/**
* @return the second partial derivative wrt {@code y}.
*/
BivariateFunction partialDerivativeYY() {
return partialDerivativeYY;
}
/**
* @return the second partial cross-derivative.
*/
BivariateFunction partialDerivativeXY() {
return partialDerivativeXY;
}
}
}