diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java | 201 |
1 files changed, 201 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java new file mode 100644 index 0000000..7f43e6f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java @@ -0,0 +1,201 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a tricubic interpolating function. + * + * @since 2.2 + * @deprecated To be removed in 4.0 (see MATH-1166). + */ +@Deprecated +public class TricubicSplineInterpolator + implements TrivariateGridInterpolator { + /** + * {@inheritDoc} + */ + public TricubicSplineInterpolatingFunction interpolate(final double[] xval, + final double[] yval, + final double[] zval, + final double[][][] fval) + throws NoDataException, NumberIsTooSmallException, + DimensionMismatchException, NonMonotonicSequenceException { + if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { + throw new NoDataException(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + MathArrays.checkOrder(zval); + + final int xLen = xval.length; + final int yLen = yval.length; + final int zLen = zval.length; + + // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets + // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k]) + // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k]) + final double[][][] fvalXY = new double[zLen][xLen][yLen]; + final double[][][] fvalZX = new double[yLen][zLen][xLen]; + for (int i = 0; i < xLen; i++) { + if (fval[i].length != yLen) { + throw new DimensionMismatchException(fval[i].length, yLen); + } + + for (int j = 0; j < yLen; j++) { + if (fval[i][j].length != zLen) { + throw new DimensionMismatchException(fval[i][j].length, zLen); + } + + for (int k = 0; k < zLen; k++) { + final double v = fval[i][j][k]; + fvalXY[k][i][j] = v; + fvalZX[j][k][i] = v; + } + } + } + + final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true); + + // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z + final BicubicSplineInterpolatingFunction[] xSplineYZ + = new BicubicSplineInterpolatingFunction[xLen]; + for (int i = 0; i < xLen; i++) { + xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); + } + + // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x + final BicubicSplineInterpolatingFunction[] ySplineZX + = new BicubicSplineInterpolatingFunction[yLen]; + for (int j = 0; j < yLen; j++) { + ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); + } + + // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y + final BicubicSplineInterpolatingFunction[] zSplineXY + = new BicubicSplineInterpolatingFunction[zLen]; + for (int k = 0; k < zLen; k++) { + zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); + } + + // Partial derivatives wrt x and wrt y + final double[][][] dFdX = new double[xLen][yLen][zLen]; + final double[][][] dFdY = new double[xLen][yLen][zLen]; + final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; + for (int k = 0; k < zLen; k++) { + final BicubicSplineInterpolatingFunction f = zSplineXY[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + dFdX[i][j][k] = f.partialDerivativeX(x, y); + dFdY[i][j][k] = f.partialDerivativeY(x, y); + d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); + } + } + } + + // Partial derivatives wrt y and wrt z + final double[][][] dFdZ = new double[xLen][yLen][zLen]; + final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; + for (int i = 0; i < xLen; i++) { + final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + dFdZ[i][j][k] = f.partialDerivativeY(y, z); + d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); + } + } + } + + // Partial derivatives wrt x and wrt z + final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; + for (int j = 0; j < yLen; j++) { + final BicubicSplineInterpolatingFunction f = ySplineZX[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); + } + } + } + + // Third partial cross-derivatives + final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; + for (int i = 0; i < xLen ; i++) { + final int nI = nextIndex(i, xLen); + final int pI = previousIndex(i); + for (int j = 0; j < yLen; j++) { + final int nJ = nextIndex(j, yLen); + final int pJ = previousIndex(j); + for (int k = 0; k < zLen; k++) { + final int nK = nextIndex(k, zLen); + final int pK = previousIndex(k); + + // XXX Not sure about this formula + d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - + fval[pI][nJ][nK] + fval[pI][pJ][nK] - + fval[nI][nJ][pK] + fval[nI][pJ][pK] + + fval[pI][nJ][pK] - fval[pI][pJ][pK]) / + ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ; + } + } + } + + // Create the interpolating splines + return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, + dFdX, dFdY, dFdZ, + d2FdXdY, d2FdZdX, d2FdYdZ, + d3FdXdYdZ); + } + + /** + * Compute the next index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i} is larger than or equal to 0. + * + * @param i Index + * @param max Upper limit of the array + * @return the next index + */ + private int nextIndex(int i, int max) { + final int index = i + 1; + return index < max ? index : index - 1; + } + /** + * Compute the previous index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i} is smaller than the size of the array. + * + * @param i Index + * @return the previous index + */ + private int previousIndex(int i) { + final int index = i - 1; + return index >= 0 ? index : 0; + } +} |