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

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

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

  58.     /** Samples x-coordinates */
  59.     private final double[] xval;
  60.     /** Samples y-coordinates */
  61.     private final double[] yval;
  62.     /** Set of cubic splines patching the whole data grid */
  63.     private final BicubicFunction[][] splines;

  64.     /**
  65.      * @param x Sample values of the x-coordinate, in increasing order.
  66.      * @param y Sample values of the y-coordinate, in increasing order.
  67.      * @param f Values of the function on every grid point.
  68.      * @param dFdX Values of the partial derivative of function with respect
  69.      * to x on every grid point.
  70.      * @param dFdY Values of the partial derivative of function with respect
  71.      * to y on every grid point.
  72.      * @param d2FdXdY Values of the cross partial derivative of function on
  73.      * every grid point.
  74.      * @throws DimensionMismatchException if the various arrays do not contain
  75.      * the expected number of elements.
  76.      * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
  77.      * not strictly increasing.
  78.      * @throws NoDataException if any of the arrays has zero length.
  79.      */
  80.     public BicubicInterpolatingFunction(double[] x,
  81.                                         double[] y,
  82.                                         double[][] f,
  83.                                         double[][] dFdX,
  84.                                         double[][] dFdY,
  85.                                         double[][] d2FdXdY)
  86.         throws DimensionMismatchException,
  87.                NoDataException,
  88.                NonMonotonicSequenceException {
  89.         final int xLen = x.length;
  90.         final int yLen = y.length;

  91.         if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
  92.             throw new NoDataException();
  93.         }
  94.         if (xLen != f.length) {
  95.             throw new DimensionMismatchException(xLen, f.length);
  96.         }
  97.         if (xLen != dFdX.length) {
  98.             throw new DimensionMismatchException(xLen, dFdX.length);
  99.         }
  100.         if (xLen != dFdY.length) {
  101.             throw new DimensionMismatchException(xLen, dFdY.length);
  102.         }
  103.         if (xLen != d2FdXdY.length) {
  104.             throw new DimensionMismatchException(xLen, d2FdXdY.length);
  105.         }

  106.         MathArrays.checkOrder(x);
  107.         MathArrays.checkOrder(y);

  108.         xval = x.clone();
  109.         yval = y.clone();

  110.         final int lastI = xLen - 1;
  111.         final int lastJ = yLen - 1;
  112.         splines = new BicubicFunction[lastI][lastJ];

  113.         for (int i = 0; i < lastI; i++) {
  114.             if (f[i].length != yLen) {
  115.                 throw new DimensionMismatchException(f[i].length, yLen);
  116.             }
  117.             if (dFdX[i].length != yLen) {
  118.                 throw new DimensionMismatchException(dFdX[i].length, yLen);
  119.             }
  120.             if (dFdY[i].length != yLen) {
  121.                 throw new DimensionMismatchException(dFdY[i].length, yLen);
  122.             }
  123.             if (d2FdXdY[i].length != yLen) {
  124.                 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
  125.             }
  126.             final int ip1 = i + 1;
  127.             final double xR = xval[ip1] - xval[i];
  128.             for (int j = 0; j < lastJ; j++) {
  129.                 final int jp1 = j + 1;
  130.                 final double yR = yval[jp1] - yval[j];
  131.                 final double xRyR = xR * yR;
  132.                 final double[] beta = new double[] {
  133.                     f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
  134.                     dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
  135.                     dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
  136.                     d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
  137.                 };

  138.                 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta));
  139.             }
  140.         }
  141.     }

  142.     /**
  143.      * {@inheritDoc}
  144.      */
  145.     public double value(double x, double y)
  146.         throws OutOfRangeException {
  147.         final int i = searchIndex(x, xval);
  148.         final int j = searchIndex(y, yval);

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

  151.         return splines[i][j].value(xN, yN);
  152.     }

  153.     /**
  154.      * Indicates whether a point is within the interpolation range.
  155.      *
  156.      * @param x First coordinate.
  157.      * @param y Second coordinate.
  158.      * @return {@code true} if (x, y) is a valid point.
  159.      */
  160.     public boolean isValidPoint(double x, double y) {
  161.         if (x < xval[0] ||
  162.             x > xval[xval.length - 1] ||
  163.             y < yval[0] ||
  164.             y > yval[yval.length - 1]) {
  165.             return false;
  166.         } else {
  167.             return true;
  168.         }
  169.     }

  170.     /**
  171.      * @param c Coordinate.
  172.      * @param val Coordinate samples.
  173.      * @return the index in {@code val} corresponding to the interval
  174.      * containing {@code c}.
  175.      * @throws OutOfRangeException if {@code c} is out of the
  176.      * range defined by the boundary values of {@code val}.
  177.      */
  178.     private int searchIndex(double c, double[] val) {
  179.         final int r = Arrays.binarySearch(val, c);

  180.         if (r == -1 ||
  181.             r == -val.length - 1) {
  182.             throw new OutOfRangeException(c, val[0], val[val.length - 1]);
  183.         }

  184.         if (r < 0) {
  185.             // "c" in within an interpolation sub-interval: Return the
  186.             // index of the sample at the lower end of the sub-interval.
  187.             return -r - 2;
  188.         }
  189.         final int last = val.length - 1;
  190.         if (r == last) {
  191.             // "c" is the last sample of the range: Return the index
  192.             // of the sample at the lower end of the last sub-interval.
  193.             return last - 1;
  194.         }

  195.         // "c" is another sample point.
  196.         return r;
  197.     }

  198.     /**
  199.      * Compute the spline coefficients from the list of function values and
  200.      * function partial derivatives values at the four corners of a grid
  201.      * element. They must be specified in the following order:
  202.      * <ul>
  203.      *  <li>f(0,0)</li>
  204.      *  <li>f(1,0)</li>
  205.      *  <li>f(0,1)</li>
  206.      *  <li>f(1,1)</li>
  207.      *  <li>f<sub>x</sub>(0,0)</li>
  208.      *  <li>f<sub>x</sub>(1,0)</li>
  209.      *  <li>f<sub>x</sub>(0,1)</li>
  210.      *  <li>f<sub>x</sub>(1,1)</li>
  211.      *  <li>f<sub>y</sub>(0,0)</li>
  212.      *  <li>f<sub>y</sub>(1,0)</li>
  213.      *  <li>f<sub>y</sub>(0,1)</li>
  214.      *  <li>f<sub>y</sub>(1,1)</li>
  215.      *  <li>f<sub>xy</sub>(0,0)</li>
  216.      *  <li>f<sub>xy</sub>(1,0)</li>
  217.      *  <li>f<sub>xy</sub>(0,1)</li>
  218.      *  <li>f<sub>xy</sub>(1,1)</li>
  219.      * </ul>
  220.      * where the subscripts indicate the partial derivative with respect to
  221.      * the corresponding variable(s).
  222.      *
  223.      * @param beta List of function values and function partial derivatives
  224.      * values.
  225.      * @return the spline coefficients.
  226.      */
  227.     private double[] computeSplineCoefficients(double[] beta) {
  228.         final double[] a = new double[NUM_COEFF];

  229.         for (int i = 0; i < NUM_COEFF; i++) {
  230.             double result = 0;
  231.             final double[] row = AINV[i];
  232.             for (int j = 0; j < NUM_COEFF; j++) {
  233.                 result += row[j] * beta[j];
  234.             }
  235.             a[i] = result;
  236.         }

  237.         return a;
  238.     }
  239. }

  240. /**
  241.  * Bicubic function.
  242.  */
  243. class BicubicFunction implements BivariateFunction {
  244.     /** Number of points. */
  245.     private static final short N = 4;
  246.     /** Coefficients */
  247.     private final double[][] a;

  248.     /**
  249.      * Simple constructor.
  250.      *
  251.      * @param coeff Spline coefficients.
  252.      */
  253.     public BicubicFunction(double[] coeff) {
  254.         a = new double[N][N];
  255.         for (int j = 0; j < N; j++) {
  256.             final double[] aJ = a[j];
  257.             for (int i = 0; i < N; i++) {
  258.                 aJ[i] = coeff[i * N + j];
  259.             }
  260.         }
  261.     }

  262.     /**
  263.      * {@inheritDoc}
  264.      */
  265.     public double value(double x, double y) {
  266.         if (x < 0 || x > 1) {
  267.             throw new OutOfRangeException(x, 0, 1);
  268.         }
  269.         if (y < 0 || y > 1) {
  270.             throw new OutOfRangeException(y, 0, 1);
  271.         }

  272.         final double x2 = x * x;
  273.         final double x3 = x2 * x;
  274.         final double[] pX = {1, x, x2, x3};

  275.         final double y2 = y * y;
  276.         final double y3 = y2 * y;
  277.         final double[] pY = {1, y, y2, y3};

  278.         return apply(pX, pY, a);
  279.     }

  280.     /**
  281.      * Compute the value of the bicubic polynomial.
  282.      *
  283.      * @param pX Powers of the x-coordinate.
  284.      * @param pY Powers of the y-coordinate.
  285.      * @param coeff Spline coefficients.
  286.      * @return the interpolated value.
  287.      */
  288.     private double apply(double[] pX, double[] pY, double[][] coeff) {
  289.         double result = 0;
  290.         for (int i = 0; i < N; i++) {
  291.             final double r = MathArrays.linearCombination(coeff[i], pY);
  292.             result += r * pX[i];
  293.         }

  294.         return result;
  295.     }
  296. }