/*
 * 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.math3.analysis.interpolation;

import java.util.Arrays;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.util.MathArrays;

Function that implements the bicubic spline interpolation.
Since:3.4
/** * 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.
/** Number of coefficients. */
private static final int NUM_COEFF = 16;
Matrix to compute the spline coefficients from the function values and function derivatives values
/** * 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
/** Samples x-coordinates */
private final double[] xval;
Samples y-coordinates
/** Samples y-coordinates */
private final double[] yval;
Set of cubic splines patching the whole data grid
/** Set of cubic splines patching the whole data grid */
private final BicubicFunction[][] splines;
Params:
  • x – Sample values of the x-coordinate, in increasing order.
  • y – Sample values of the y-coordinate, in increasing order.
  • f – Values of the function on every grid point.
  • dFdX – Values of the partial derivative of function with respect to x on every grid point.
  • dFdY – Values of the partial derivative of function with respect to y on every grid point.
  • d2FdXdY – Values of the cross partial derivative of function on every grid point.
Throws:
/** * @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 { 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)); } } }
{@inheritDoc}
/** * {@inheritDoc} */
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.
Params:
  • x – First coordinate.
  • y – Second coordinate.
Returns:true if (x, y) is a valid point.
/** * 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) { if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] || y > yval[yval.length - 1]) { return false; } else { return true; } }
Params:
  • c – Coordinate.
  • val – Coordinate samples.
Throws:
Returns:the index in val corresponding to the interval containing c.
/** * @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 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:
  • f(0,0)
  • f(1,0)
  • f(0,1)
  • f(1,1)
  • fx(0,0)
  • fx(1,0)
  • fx(0,1)
  • fx(1,1)
  • fy(0,0)
  • fy(1,0)
  • fy(0,1)
  • fy(1,1)
  • fxy(0,0)
  • fxy(1,0)
  • fxy(0,1)
  • fxy(1,1)
where the subscripts indicate the partial derivative with respect to the corresponding variable(s).
Params:
  • beta – List of function values and function partial derivatives values.
Returns:the spline coefficients.
/** * 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 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.
/** * Bicubic function. */
class BicubicFunction implements BivariateFunction {
Number of points.
/** Number of points. */
private static final short N = 4;
Coefficients
/** Coefficients */
private final double[][] a;
Simple constructor.
Params:
  • coeff – Spline coefficients.
/** * Simple constructor. * * @param coeff Spline coefficients. */
BicubicFunction(double[] coeff) { 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]; } } }
{@inheritDoc}
/** * {@inheritDoc} */
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, pY, a); }
Compute the value of the bicubic polynomial.
Params:
  • pX – Powers of the x-coordinate.
  • pY – Powers of the y-coordinate.
  • coeff – Spline coefficients.
Returns:the interpolated value.
/** * Compute the value of the bicubic polynomial. * * @param pX Powers of the x-coordinate. * @param pY Powers of the y-coordinate. * @param coeff Spline coefficients. * @return the interpolated value. */
private double apply(double[] pX, double[] pY, double[][] coeff) { double result = 0; for (int i = 0; i < N; i++) { final double r = MathArrays.linearCombination(coeff[i], pY); result += r * pX[i]; } return result; } }