BicubicInterpolatingFunction.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.analysis.interpolation;

  19. import java.util.Arrays;
  20. import java.util.function.DoubleBinaryOperator;
  21. import java.util.function.Function;

  23. import org.apache.commons.numbers.core.Sum;
  24. import org.apache.commons.math4.legacy.analysis.BivariateFunction;
  25. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  26. import org.apache.commons.math4.legacy.exception.NoDataException;
  27. import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
  28. import org.apache.commons.math4.legacy.exception.OutOfRangeException;
  29. import org.apache.commons.math4.legacy.core.MathArrays;

  31. /**
  32.  * Function that implements the
  33.  * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
  34.  * bicubic spline interpolation</a>.
  35.  *
  36.  * @since 3.4
  37.  */
  38. public class BicubicInterpolatingFunction
  39.     implements BivariateFunction {
  40.     /** Number of coefficients. */
  41.     private static final int NUM_COEFF = 16;
  42.     /**
  43.      * Matrix to compute the spline coefficients from the function values
  44.      * and function derivatives values.
  45.      */
  46.     private static final double[][] AINV = {
  47.         { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  48.         { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
  49.         { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
  50.         { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
  51.         { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
  52.         { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
  53.         { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
  54.         { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
  55.         { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
  56.         { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
  57.         { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
  58.         { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
  59.         { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
  60.         { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
  61.         { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
  62.         { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
  63.     };

  65.     /** Samples x-coordinates. */
  66.     private final double[] xval;
  67.     /** Samples y-coordinates. */
  68.     private final double[] yval;
  69.     /** Set of cubic splines patching the whole data grid. */
  70.     private final BicubicFunction[][] splines;

  72.     /**
  73.      * @param x Sample values of the x-coordinate, in increasing order.
  74.      * @param y Sample values of the y-coordinate, in increasing order.
  75.      * @param f Values of the function on every grid point.
  76.      * @param dFdX Values of the partial derivative of function with respect
  77.      * to x on every grid point.
  78.      * @param dFdY Values of the partial derivative of function with respect
  79.      * to y on every grid point.
  80.      * @param d2FdXdY Values of the cross partial derivative of function on
  81.      * every grid point.
  82.      * @throws DimensionMismatchException if the various arrays do not contain
  83.      * the expected number of elements.
  84.      * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
  85.      * not strictly increasing.
  86.      * @throws NoDataException if any of the arrays has zero length.
  87.      */
  88.     public BicubicInterpolatingFunction(double[] x,
  89.                                         double[] y,
  90.                                         double[][] f,
  91.                                         double[][] dFdX,
  92.                                         double[][] dFdY,
  93.                                         double[][] d2FdXdY)
  94.         throws DimensionMismatchException,
  95.                NoDataException,
  96.                NonMonotonicSequenceException {
  97.         this(x, y, f, dFdX, dFdY, d2FdXdY, false);
  98.     }

  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;

  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.         }

  148.         MathArrays.checkOrder(x);
  149.         MathArrays.checkOrder(y);

  151.         xval = x.clone();
  152.         yval = y.clone();

  154.         final int lastI = xLen - 1;
  155.         final int lastJ = yLen - 1;
  156.         splines = new BicubicFunction[lastI][lastJ];

  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.                 };

  184.                 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta),
  185.                                                     xR,
  186.                                                     yR,
  187.                                                     initializeDerivatives);
  188.             }
  189.         }
  190.     }

  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);

  201.         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
  202.         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);

  204.         return splines[i][j].value(xN, yN);
  205.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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);

  283.             final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
  284.             final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);

  286.             return which.apply(splines[i][j]).value(xN, yN);
  287.         };
  288.     }

  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);

  301.         if (r == -1 ||
  302.             r == -val.length - 1) {
  303.             throw new OutOfRangeException(c, val[0], val[val.length - 1]);
  304.         }

  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.         }

  318.         // "c" is another sample point.
  319.         return r;
  320.     }

  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];

  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.         }

  363.         return a;
  364.     }

  366.     /**
  367.      * Bicubic function.
  368.      */
  369.     static final class BicubicFunction implements BivariateFunction {
  370.         /** Number of points. */
  371.         private static final short N = 4;
  372.         /** Coefficients. */
  373.         private final double[][] a;
  374.         /** First partial derivative along x. */
  375.         private final BivariateFunction partialDerivativeX;
  376.         /** First partial derivative along y. */
  377.         private final BivariateFunction partialDerivativeY;
  378.         /** Second partial derivative along x. */
  379.         private final BivariateFunction partialDerivativeXX;
  380.         /** Second partial derivative along y. */
  381.         private final BivariateFunction partialDerivativeYY;
  382.         /** Second crossed partial derivative. */
  383.         private final BivariateFunction partialDerivativeXY;

  385.         /**
  386.          * Simple constructor.
  387.          *
  388.          * @param coeff Spline coefficients.
  389.          * @param xR x spacing.
  390.          * @param yR y spacing.
  391.          * @param initializeDerivatives Whether to initialize the internal data
  392.          * needed for calling any of the methods that compute the partial derivatives
  393.          * this function.
  394.          */
  395.         BicubicFunction(double[] coeff,
  396.                         double xR,
  397.                         double yR,
  398.                         boolean initializeDerivatives) {
  399.             a = new double[N][N];
  400.             for (int j = 0; j < N; j++) {
  401.                 final double[] aJ = a[j];
  402.                 for (int i = 0; i < N; i++) {
  403.                     aJ[i] = coeff[i * N + j];
  404.                 }
  405.             }

  407.             if (initializeDerivatives) {
  408.                 // Compute all partial derivatives functions.
  409.                 final double[][] aX = new double[N][N];
  410.                 final double[][] aY = new double[N][N];
  411.                 final double[][] aXX = new double[N][N];
  412.                 final double[][] aYY = new double[N][N];
  413.                 final double[][] aXY = new double[N][N];

  415.                 for (int i = 0; i < N; i++) {
  416.                     for (int j = 0; j < N; j++) {
  417.                         final double c = a[i][j];
  418.                         aX[i][j] = i * c;
  419.                         aY[i][j] = j * c;
  420.                         aXX[i][j] = (i - 1) * aX[i][j];
  421.                         aYY[i][j] = (j - 1) * aY[i][j];
  422.                         aXY[i][j] = j * aX[i][j];
  423.                     }
  424.                 }

  426.                 partialDerivativeX = (double x, double y) -> {
  427.                     final double x2 = x * x;
  428.                     final double[] pX = {0, 1, x, x2};

  430.                     final double y2 = y * y;
  431.                     final double y3 = y2 * y;
  432.                     final double[] pY = {1, y, y2, y3};

  434.                     return apply(pX, 1, pY, 0, aX) / xR;
  435.                 };
  436.                 partialDerivativeY = (double x, double y) -> {
  437.                     final double x2 = x * x;
  438.                     final double x3 = x2 * x;
  439.                     final double[] pX = {1, x, x2, x3};

  441.                     final double y2 = y * y;
  442.                     final double[] pY = {0, 1, y, y2};

  444.                     return apply(pX, 0, pY, 1, aY) / yR;
  445.                 };
  446.                 partialDerivativeXX = (double x, double y) -> {
  447.                     final double[] pX = {0, 0, 1, x};

  449.                     final double y2 = y * y;
  450.                     final double y3 = y2 * y;
  451.                     final double[] pY = {1, y, y2, y3};

  453.                     return apply(pX, 2, pY, 0, aXX) / (xR * xR);
  454.                 };
  455.                 partialDerivativeYY = (double x, double y) -> {
  456.                     final double x2 = x * x;
  457.                     final double x3 = x2 * x;
  458.                     final double[] pX = {1, x, x2, x3};

  460.                     final double[] pY = {0, 0, 1, y};

  462.                     return apply(pX, 0, pY, 2, aYY) / (yR * yR);
  463.                 };
  464.                 partialDerivativeXY = (double x, double y) -> {
  465.                     final double x2 = x * x;
  466.                     final double[] pX = {0, 1, x, x2};

  468.                     final double y2 = y * y;
  469.                     final double[] pY = {0, 1, y, y2};

  471.                     return apply(pX, 1, pY, 1, aXY) / (xR * yR);
  472.                 };
  473.             } else {
  474.                 partialDerivativeX = null;
  475.                 partialDerivativeY = null;
  476.                 partialDerivativeXX = null;
  477.                 partialDerivativeYY = null;
  478.                 partialDerivativeXY = null;
  479.             }
  480.         }

  482.         /**
  483.          * {@inheritDoc}
  484.          */
  485.         @Override
  486.         public double value(double x, double y) {
  487.             if (x < 0 || x > 1) {
  488.                 throw new OutOfRangeException(x, 0, 1);
  489.             }
  490.             if (y < 0 || y > 1) {
  491.                 throw new OutOfRangeException(y, 0, 1);
  492.             }

  494.             final double x2 = x * x;
  495.             final double x3 = x2 * x;
  496.             final double[] pX = {1, x, x2, x3};

  498.             final double y2 = y * y;
  499.             final double y3 = y2 * y;
  500.             final double[] pY = {1, y, y2, y3};

  502.             return apply(pX, 0, pY, 0, a);
  503.         }

  505.         /**
  506.          * Compute the value of the bicubic polynomial.
  507.          *
  508.          * <p>Assumes the powers are zero below the provided index, and 1 at the provided
  509.          * index. This allows skipping some zero products and optimising multiplication
  510.          * by one.
  511.          *
  512.          * @param pX Powers of the x-coordinate.
  513.          * @param i Index of pX[i] == 1
  514.          * @param pY Powers of the y-coordinate.
  515.          * @param j Index of pX[j] == 1
  516.          * @param coeff Spline coefficients.
  517.          * @return the interpolated value.
  518.          */
  519.         private static double apply(double[] pX, int i, double[] pY, int j, double[][] coeff) {
  520.             // assert pX[i] == 1
  521.             double result = sumOfProducts(coeff[i], pY, j);
  522.             while (++i < N) {
  523.                 final double r = sumOfProducts(coeff[i], pY, j);
  524.                 result += r * pX[i];
  525.             }
  526.             return result;
  527.         }

  529.         /**
  530.          * Compute the sum of products starting from the provided index.
  531.          * Assumes that factor {@code b[j] == 1}.
  532.          *
  533.          * @param a Factors.
  534.          * @param b Factors.
  535.          * @param j Index to initialise the sum.
  536.          * @return the double
  537.          */
  538.         private static double sumOfProducts(double[] a, double[] b, int j) {
  539.             // assert b[j] == 1
  540.             final Sum sum = Sum.of(a[j]);
  541.             while (++j < N) {
  542.                 sum.addProduct(a[j], b[j]);
  543.             }
  544.             return sum.getAsDouble();
  545.         }

  547.         /**
  548.          * @return the partial derivative wrt {@code x}.
  549.          */
  550.         BivariateFunction partialDerivativeX() {
  551.             return partialDerivativeX;
  552.         }

  554.         /**
  555.          * @return the partial derivative wrt {@code y}.
  556.          */
  557.         BivariateFunction partialDerivativeY() {
  558.             return partialDerivativeY;
  559.         }

  561.         /**
  562.          * @return the second partial derivative wrt {@code x}.
  563.          */
  564.         BivariateFunction partialDerivativeXX() {
  565.             return partialDerivativeXX;
  566.         }

  568.         /**
  569.          * @return the second partial derivative wrt {@code y}.
  570.          */
  571.         BivariateFunction partialDerivativeYY() {
  572.             return partialDerivativeYY;
  573.         }

  575.         /**
  576.          * @return the second partial cross-derivative.
  577.          */
  578.         BivariateFunction partialDerivativeXY() {
  579.             return partialDerivativeXY;
  580.         }
  581.     }
  582. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT