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

  19. import java.util.Arrays;
  20. import org.apache.commons.math3.analysis.BivariateFunction;
  21. import org.apache.commons.math3.exception.DimensionMismatchException;
  22. import org.apache.commons.math3.exception.NoDataException;
  23. import org.apache.commons.math3.exception.OutOfRangeException;
  24. import org.apache.commons.math3.exception.NonMonotonicSequenceException;
  25. import org.apache.commons.math3.util.MathArrays;

  27. /**
  28.  * Function that implements the
  29.  * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
  30.  * bicubic spline interpolation</a>. Due to numerical accuracy issues this should not
  31.  * be used.
  32.  *
  33.  * @since 2.1
  34.  * @deprecated as of 3.4 replaced by
  35.  * {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction}
  36.  */
  37. @Deprecated
  38. public class BicubicSplineInterpolatingFunction
  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 BicubicSplineFunction[][] splines;
  71.     /**
  72.      * Partial derivatives.
  73.      * The value of the first index determines the kind of derivatives:
  74.      * 0 = first partial derivatives wrt x
  75.      * 1 = first partial derivatives wrt y
  76.      * 2 = second partial derivatives wrt x
  77.      * 3 = second partial derivatives wrt y
  78.      * 4 = cross partial derivatives
  79.      */
  80.     private final BivariateFunction[][][] partialDerivatives;

  82.     /**
  83.      * @param x Sample values of the x-coordinate, in increasing order.
  84.      * @param y Sample values of the y-coordinate, in increasing order.
  85.      * @param f Values of the function on every grid point.
  86.      * @param dFdX Values of the partial derivative of function with respect
  87.      * to x on every grid point.
  88.      * @param dFdY Values of the partial derivative of function with respect
  89.      * to y on every grid point.
  90.      * @param d2FdXdY Values of the cross partial derivative of function on
  91.      * every grid point.
  92.      * @throws DimensionMismatchException if the various arrays do not contain
  93.      * the expected number of elements.
  94.      * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
  95.      * not strictly increasing.
  96.      * @throws NoDataException if any of the arrays has zero length.
  97.      */
  98.     public BicubicSplineInterpolatingFunction(double[] x,
  99.                                               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.     }

  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;

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

  164.         MathArrays.checkOrder(x);
  165.         MathArrays.checkOrder(y);

  167.         xval = x.clone();
  168.         yval = y.clone();

  170.         final int lastI = xLen - 1;
  171.         final int lastJ = yLen - 1;
  172.         splines = new BicubicSplineFunction[lastI][lastJ];

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

  197.                 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
  198.                                                           initializeDerivatives);
  199.             }
  200.         }

  202.         if (initializeDerivatives) {
  203.             // Compute all partial derivatives.
  204.             partialDerivatives = new BivariateFunction[5][lastI][lastJ];

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

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

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

  233.         return splines[i][j].value(xN, yN);
  234.     }

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

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

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

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

  355.         return partialDerivatives[which][i][j].value(xN, yN);
  356.     }

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

  369.         if (r == -1 ||
  370.             r == -val.length - 1) {
  371.             throw new OutOfRangeException(c, val[0], val[val.length - 1]);
  372.         }

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

  386.         // "c" is another sample point.
  387.         return r;
  388.     }

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

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

  431.         return a;
  432.     }
  433. }

  435. /**
  436.  * 2D-spline function.
  437.  *
  438.  */
  439. class 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;

  455.     /**
  456.      * Simple constructor.
  457.      *
  458.      * @param coeff Spline coefficients.
  459.      */
  460.     public BicubicSplineFunction(double[] coeff) {
  461.         this(coeff, false);
  462.     }

  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.     public BicubicSplineFunction(double[] coeff,
  473.                                  boolean initializeDerivatives) {
  474.         a = new double[N][N];
  475.         for (int i = 0; i < N; i++) {
  476.             for (int j = 0; j < N; j++) {
  477.                 a[i][j] = coeff[i * N + j];
  478.             }
  479.         }

  481.         if (initializeDerivatives) {
  482.             // Compute all partial derivatives functions.
  483.             final double[][] aX = new double[N][N];
  484.             final double[][] aY = new double[N][N];
  485.             final double[][] aXX = new double[N][N];
  486.             final double[][] aYY = new double[N][N];
  487.             final double[][] aXY = new double[N][N];

  489.             for (int i = 0; i < N; i++) {
  490.                 for (int j = 0; j < N; j++) {
  491.                     final double c = a[i][j];
  492.                     aX[i][j] = i * c;
  493.                     aY[i][j] = j * c;
  494.                     aXX[i][j] = (i - 1) * aX[i][j];
  495.                     aYY[i][j] = (j - 1) * aY[i][j];
  496.                     aXY[i][j] = j * aX[i][j];
  497.                 }
  498.             }

  500.             partialDerivativeX = new BivariateFunction() {
  501.                     public double value(double x, double y)  {
  502.                         final double x2 = x * x;
  503.                         final double[] pX = {0, 1, x, x2};

  505.                         final double y2 = y * y;
  506.                         final double y3 = y2 * y;
  507.                         final double[] pY = {1, y, y2, y3};

  509.                         return apply(pX, pY, aX);
  510.                     }
  511.                 };
  512.             partialDerivativeY = new BivariateFunction() {
  513.                     public double value(double x, double y)  {
  514.                         final double x2 = x * x;
  515.                         final double x3 = x2 * x;
  516.                         final double[] pX = {1, x, x2, x3};

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

  521.                         return apply(pX, pY, aY);
  522.                     }
  523.                 };
  524.             partialDerivativeXX = new BivariateFunction() {
  525.                     public double value(double x, double y)  {
  526.                         final double[] pX = {0, 0, 1, x};

  528.                         final double y2 = y * y;
  529.                         final double y3 = y2 * y;
  530.                         final double[] pY = {1, y, y2, y3};

  532.                         return apply(pX, pY, aXX);
  533.                     }
  534.                 };
  535.             partialDerivativeYY = new BivariateFunction() {
  536.                     public double value(double x, double y)  {
  537.                         final double x2 = x * x;
  538.                         final double x3 = x2 * x;
  539.                         final double[] pX = {1, x, x2, x3};

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

  543.                         return apply(pX, pY, aYY);
  544.                     }
  545.                 };
  546.             partialDerivativeXY = new BivariateFunction() {
  547.                     public double value(double x, double y)  {
  548.                         final double x2 = x * x;
  549.                         final double[] pX = {0, 1, x, x2};

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

  554.                         return apply(pX, pY, aXY);
  555.                     }
  556.                 };
  557.         } else {
  558.             partialDerivativeX = null;
  559.             partialDerivativeY = null;
  560.             partialDerivativeXX = null;
  561.             partialDerivativeYY = null;
  562.             partialDerivativeXY = null;
  563.         }
  564.     }

  566.     /**
  567.      * {@inheritDoc}
  568.      */
  569.     public double value(double x, double y) {
  570.         if (x < 0 || x > 1) {
  571.             throw new OutOfRangeException(x, 0, 1);
  572.         }
  573.         if (y < 0 || y > 1) {
  574.             throw new OutOfRangeException(y, 0, 1);
  575.         }

  577.         final double x2 = x * x;
  578.         final double x3 = x2 * x;
  579.         final double[] pX = {1, x, x2, x3};

  581.         final double y2 = y * y;
  582.         final double y3 = y2 * y;
  583.         final double[] pY = {1, y, y2, y3};

  585.         return apply(pX, pY, a);
  586.     }

  588.     /**
  589.      * Compute the value of the bicubic polynomial.
  590.      *
  591.      * @param pX Powers of the x-coordinate.
  592.      * @param pY Powers of the y-coordinate.
  593.      * @param coeff Spline coefficients.
  594.      * @return the interpolated value.
  595.      */
  596.     private double apply(double[] pX, double[] pY, double[][] coeff) {
  597.         double result = 0;
  598.         for (int i = 0; i < N; i++) {
  599.             for (int j = 0; j < N; j++) {
  600.                 result += coeff[i][j] * pX[i] * pY[j];
  601.             }
  602.         }

  604.         return result;
  605.     }

  607.     /**
  608.      * @return the partial derivative wrt {@code x}.
  609.      */
  610.     public BivariateFunction partialDerivativeX() {
  611.         return partialDerivativeX;
  612.     }
  613.     /**
  614.      * @return the partial derivative wrt {@code y}.
  615.      */
  616.     public BivariateFunction partialDerivativeY() {
  617.         return partialDerivativeY;
  618.     }
  619.     /**
  620.      * @return the second partial derivative wrt {@code x}.
  621.      */
  622.     public BivariateFunction partialDerivativeXX() {
  623.         return partialDerivativeXX;
  624.     }
  625.     /**
  626.      * @return the second partial derivative wrt {@code y}.
  627.      */
  628.     public BivariateFunction partialDerivativeYY() {
  629.         return partialDerivativeYY;
  630.     }
  631.     /**
  632.      * @return the second partial cross-derivative.
  633.      */
  634.     public BivariateFunction partialDerivativeXY() {
  635.         return partialDerivativeXY;
  636.     }
  637. }
Created with JaCoCo 0.7.2.201409121644