TricubicSplineInterpolator.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 org.apache.commons.math3.exception.DimensionMismatchException;
  19. import org.apache.commons.math3.exception.NoDataException;
  20. import org.apache.commons.math3.exception.NonMonotonicSequenceException;
  21. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math3.util.MathArrays;

  23. /**
  24.  * Generates a tricubic interpolating function.
  25.  *
  26.  * @since 2.2
  27.  * @deprecated To be removed in 4.0 (see MATH-1166).
  28.  */
  29. @Deprecated
  30. public class TricubicSplineInterpolator
  31.     implements TrivariateGridInterpolator {
  32.     /**
  33.      * {@inheritDoc}
  34.      */
  35.     public TricubicSplineInterpolatingFunction interpolate(final double[] xval,
  36.                                                            final double[] yval,
  37.                                                            final double[] zval,
  38.                                                            final double[][][] fval)
  39.         throws NoDataException, NumberIsTooSmallException,
  40.                DimensionMismatchException, NonMonotonicSequenceException {
  41.         if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
  42.             throw new NoDataException();
  43.         }
  44.         if (xval.length != fval.length) {
  45.             throw new DimensionMismatchException(xval.length, fval.length);
  46.         }

  47.         MathArrays.checkOrder(xval);
  48.         MathArrays.checkOrder(yval);
  49.         MathArrays.checkOrder(zval);

  50.         final int xLen = xval.length;
  51.         final int yLen = yval.length;
  52.         final int zLen = zval.length;

  53.         // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
  54.         // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
  55.         // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
  56.         final double[][][] fvalXY = new double[zLen][xLen][yLen];
  57.         final double[][][] fvalZX = new double[yLen][zLen][xLen];
  58.         for (int i = 0; i < xLen; i++) {
  59.             if (fval[i].length != yLen) {
  60.                 throw new DimensionMismatchException(fval[i].length, yLen);
  61.             }

  62.             for (int j = 0; j < yLen; j++) {
  63.                 if (fval[i][j].length != zLen) {
  64.                     throw new DimensionMismatchException(fval[i][j].length, zLen);
  65.                 }

  66.                 for (int k = 0; k < zLen; k++) {
  67.                     final double v = fval[i][j][k];
  68.                     fvalXY[k][i][j] = v;
  69.                     fvalZX[j][k][i] = v;
  70.                 }
  71.             }
  72.         }

  73.         final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true);

  74.         // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
  75.         final BicubicSplineInterpolatingFunction[] xSplineYZ
  76.             = new BicubicSplineInterpolatingFunction[xLen];
  77.         for (int i = 0; i < xLen; i++) {
  78.             xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]);
  79.         }

  80.         // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
  81.         final BicubicSplineInterpolatingFunction[] ySplineZX
  82.             = new BicubicSplineInterpolatingFunction[yLen];
  83.         for (int j = 0; j < yLen; j++) {
  84.             ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]);
  85.         }

  86.         // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
  87.         final BicubicSplineInterpolatingFunction[] zSplineXY
  88.             = new BicubicSplineInterpolatingFunction[zLen];
  89.         for (int k = 0; k < zLen; k++) {
  90.             zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
  91.         }

  92.         // Partial derivatives wrt x and wrt y
  93.         final double[][][] dFdX = new double[xLen][yLen][zLen];
  94.         final double[][][] dFdY = new double[xLen][yLen][zLen];
  95.         final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
  96.         for (int k = 0; k < zLen; k++) {
  97.             final BicubicSplineInterpolatingFunction f = zSplineXY[k];
  98.             for (int i = 0; i < xLen; i++) {
  99.                 final double x = xval[i];
  100.                 for (int j = 0; j < yLen; j++) {
  101.                     final double y = yval[j];
  102.                     dFdX[i][j][k] = f.partialDerivativeX(x, y);
  103.                     dFdY[i][j][k] = f.partialDerivativeY(x, y);
  104.                     d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
  105.                 }
  106.             }
  107.         }

  108.         // Partial derivatives wrt y and wrt z
  109.         final double[][][] dFdZ = new double[xLen][yLen][zLen];
  110.         final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
  111.         for (int i = 0; i < xLen; i++) {
  112.             final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
  113.             for (int j = 0; j < yLen; j++) {
  114.                 final double y = yval[j];
  115.                 for (int k = 0; k < zLen; k++) {
  116.                     final double z = zval[k];
  117.                     dFdZ[i][j][k] = f.partialDerivativeY(y, z);
  118.                     d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
  119.                 }
  120.             }
  121.         }

  122.         // Partial derivatives wrt x and wrt z
  123.         final double[][][] d2FdZdX = new double[xLen][yLen][zLen];
  124.         for (int j = 0; j < yLen; j++) {
  125.             final BicubicSplineInterpolatingFunction f = ySplineZX[j];
  126.             for (int k = 0; k < zLen; k++) {
  127.                 final double z = zval[k];
  128.                 for (int i = 0; i < xLen; i++) {
  129.                     final double x = xval[i];
  130.                     d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
  131.                 }
  132.             }
  133.         }

  134.         // Third partial cross-derivatives
  135.         final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
  136.         for (int i = 0; i < xLen ; i++) {
  137.             final int nI = nextIndex(i, xLen);
  138.             final int pI = previousIndex(i);
  139.             for (int j = 0; j < yLen; j++) {
  140.                 final int nJ = nextIndex(j, yLen);
  141.                 final int pJ = previousIndex(j);
  142.                 for (int k = 0; k < zLen; k++) {
  143.                     final int nK = nextIndex(k, zLen);
  144.                     final int pK = previousIndex(k);

  145.                     // XXX Not sure about this formula
  146.                     d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
  147.                                           fval[pI][nJ][nK] + fval[pI][pJ][nK] -
  148.                                           fval[nI][nJ][pK] + fval[nI][pJ][pK] +
  149.                                           fval[pI][nJ][pK] - fval[pI][pJ][pK]) /
  150.                         ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ;
  151.                 }
  152.             }
  153.         }

  154.         // Create the interpolating splines
  155.         return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval,
  156.                                                        dFdX, dFdY, dFdZ,
  157.                                                        d2FdXdY, d2FdZdX, d2FdYdZ,
  158.                                                        d3FdXdYdZ);
  159.     }

  160.     /**
  161.      * Compute the next index of an array, clipping if necessary.
  162.      * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
  163.      *
  164.      * @param i Index
  165.      * @param max Upper limit of the array
  166.      * @return the next index
  167.      */
  168.     private int nextIndex(int i, int max) {
  169.         final int index = i + 1;
  170.         return index < max ? index : index - 1;
  171.     }
  172.     /**
  173.      * Compute the previous index of an array, clipping if necessary.
  174.      * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
  175.      *
  176.      * @param i Index
  177.      * @return the previous index
  178.      */
  179.     private int previousIndex(int i) {
  180.         final int index = i - 1;
  181.         return index >= 0 ? index : 0;
  182.     }
  183. }