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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.analysis.interpolation;

  24. import java.util.Arrays;

  26. import org.hipparchus.analysis.BivariateFunction;
  27. import org.hipparchus.exception.LocalizedCoreFormats;
  28. import org.hipparchus.exception.MathIllegalArgumentException;
  29. import org.hipparchus.util.MathArrays;
  30. import org.hipparchus.util.MathUtils;

  32. /**
  33.  * Function that implements the
  34.  * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
  35.  * bicubic spline interpolation</a>.
  36.  *
  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.     /** Simple constructor.
  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 MathIllegalArgumentException if the various arrays do not contain
  83.      * the expected number of elements.
  84.      * @throws MathIllegalArgumentException if {@code x} or {@code y} are
  85.      * not strictly increasing.
  86.      * @throws MathIllegalArgumentException 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 MathIllegalArgumentException {
  95.         final int xLen = x.length;
  96.         final int yLen = y.length;

  98.         if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
  99.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  100.         }
  101.         MathUtils.checkDimension(xLen, f.length);
  102.         MathUtils.checkDimension(xLen, dFdX.length);
  103.         MathUtils.checkDimension(xLen, dFdY.length);
  104.         MathUtils.checkDimension(xLen, d2FdXdY.length);
  105.         MathArrays.checkOrder(x);
  106.         MathArrays.checkOrder(y);

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

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

  115.         for (int i = 0; i < lastI; i++) {
  116.             MathUtils.checkDimension(f[i].length, yLen);
  117.             MathUtils.checkDimension(dFdX[i].length, yLen);
  118.             MathUtils.checkDimension(dFdY[i].length, yLen);
  119.             MathUtils.checkDimension(d2FdXdY[i].length, yLen);

  121.             final int ip1 = i + 1;
  122.             final double xR = xval[ip1] - xval[i];
  123.             for (int j = 0; j < lastJ; j++) {
  124.                 final int jp1 = j + 1;
  125.                 final double yR = yval[jp1] - yval[j];
  126.                 final double xRyR = xR * yR;
  127.                 final double[] beta = {
  128.                     f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
  129.                     dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
  130.                     dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
  131.                     d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
  132.                 };

  134.                 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta));
  135.             }
  136.         }
  137.     }

  139.     /**
  140.      * {@inheritDoc}
  141.      */
  142.     @Override
  143.     public double value(double x, double y)
  144.         throws MathIllegalArgumentException {
  145.         final int i = searchIndex(x, xval);
  146.         final int j = searchIndex(y, yval);

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

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

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

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

  183.         if (r == -1 ||
  184.             r == -val.length - 1) {
  185.             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  186.                                                    c, val[0], val[val.length - 1]);
  187.         }

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

  201.         // "c" is another sample point.
  202.         return r;
  203.     }

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

  237.         for (int i = 0; i < NUM_COEFF; i++) {
  238.             double result = 0;
  239.             final double[] row = AINV[i];
  240.             for (int j = 0; j < NUM_COEFF; j++) {
  241.                 result += row[j] * beta[j];
  242.             }
  243.             a[i] = result;
  244.         }

  246.         return a;
  247.     }
  248. }

  250. /**
  251.  * Bicubic function.
  252.  */
  253. class BicubicFunction implements BivariateFunction {
  254.     /** Number of points. */
  255.     private static final short N = 4;
  256.     /** Coefficients */
  257.     private final double[][] a;

  259.     /**
  260.      * Simple constructor.
  261.      *
  262.      * @param coeff Spline coefficients.
  263.      */
  264.     BicubicFunction(double[] coeff) {
  265.         a = new double[N][N];
  266.         for (int j = 0; j < N; j++) {
  267.             final double[] aJ = a[j];
  268.             for (int i = 0; i < N; i++) {
  269.                 aJ[i] = coeff[i * N + j];
  270.             }
  271.         }
  272.     }

  274.     /**
  275.      * {@inheritDoc}
  276.      */
  277.     @Override
  278.     public double value(double x, double y) {
  279.         MathUtils.checkRangeInclusive(x, 0, 1);
  280.         MathUtils.checkRangeInclusive(y, 0, 1);

  282.         final double x2 = x * x;
  283.         final double x3 = x2 * x;
  284.         final double[] pX = {1, x, x2, x3};

  286.         final double y2 = y * y;
  287.         final double y3 = y2 * y;
  288.         final double[] pY = {1, y, y2, y3};

  290.         return apply(pX, pY, a);
  291.     }

  293.     /**
  294.      * Compute the value of the bicubic polynomial.
  295.      *
  296.      * @param pX Powers of the x-coordinate.
  297.      * @param pY Powers of the y-coordinate.
  298.      * @param coeff Spline coefficients.
  299.      * @return the interpolated value.
  300.      */
  301.     private double apply(double[] pX, double[] pY, double[][] coeff) {
  302.         double result = 0;
  303.         for (int i = 0; i < N; i++) {
  304.             final double r = MathArrays.linearCombination(coeff[i], pY);
  305.             result += r * pX[i];
  306.         }

  308.         return result;
  309.     }
  310. }
Created with JaCoCo 0.8.12.202403310830