TricubicSplineInterpolator implements TrivariateGridInterpolator { 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; 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); final BicubicSplineInterpolatingFunction[] xSplineYZ = new BicubicSplineInterpolatingFunction[xLen]; for (int i = 0; i < xLen; i++) { xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); } final BicubicSplineInterpolatingFunction[] ySplineZX = new BicubicSplineInterpolatingFunction[yLen]; for (int j = 0; j < yLen; j++) { ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); } final BicubicSplineInterpolatingFunction[] zSplineXY = new BicubicSplineInterpolatingFunction[zLen]; for (int k = 0; k < zLen; k++) { zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); } 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); } } } 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); } } } 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); } } } 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); 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])) ; } } } return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, dFdX, dFdY, dFdZ, d2FdXdY, d2FdZdX, d2FdYdZ, d3FdXdYdZ); }  TricubicSplineInterpolatingFunction interpolate(final double[] xval,
                                                           final double[] yval,
                                                           final double[] zval,
                                                           final double[][][] fval);  }
@Test public void testPreconditions() { double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2.5}; double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5}; double[][][] fval = new double[xval.length][yval.length][zval.length]; TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator(); @SuppressWarnings("unused") TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval); double[] wxval = new double[] {3, 2, 5, 6.5}; try { p = interpolator.interpolate(wxval, yval, zval, fval); Assert.fail("an exception should have been thrown"); } catch (MathIllegalArgumentException e) { } double[] wyval = new double[] {-4, -3, -1, -1}; try { p = interpolator.interpolate(xval, wyval, zval, fval); Assert.fail("an exception should have been thrown"); } catch (MathIllegalArgumentException e) { } double[] wzval = new double[] {-12, -8, -5.5, -3, -4, 2.5}; try { p = interpolator.interpolate(xval, yval, wzval, fval); Assert.fail("an exception should have been thrown"); } catch (MathIllegalArgumentException e) { } double[][][] wfval = new double[xval.length][yval.length + 1][zval.length]; try { p = interpolator.interpolate(xval, yval, zval, wfval); Assert.fail("an exception should have been thrown"); } catch (DimensionMismatchException e) { } wfval = new double[xval.length - 1][yval.length][zval.length]; try { p = interpolator.interpolate(xval, yval, zval, wfval); Assert.fail("an exception should have been thrown"); } catch (DimensionMismatchException e) { } wfval = new double[xval.length][yval.length][zval.length - 1]; try { p = interpolator.interpolate(xval, yval, zval, wfval); Assert.fail("an exception should have been thrown"); } catch (DimensionMismatchException e) { } }
@Ignore@Test public void testPlane() { TrivariateFunction f = new TrivariateFunction() { public double value(double x, double y, double z) { return 2 * x - 3 * y - z + 5; } }; TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator(); double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2, 2.5}; double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5}; double[][][] fval = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { fval[i][j][k] = f.value(xval[i], yval[j], zval[k]); } } } TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval); double x, y, z; double expected, result; x = 4; y = -3; z = 0; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("On sample point", expected, result, 1e-15); x = 4.5; y = -1.5; z = -4.25; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 0.3); x = 3.5; y = -3.5; z = -10; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 0.3); }
@Ignore@Test public void testWave() { double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2, 2.5}; double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4}; final double a = 0.2; final double omega = 0.5; final double kx = 2; final double ky = 1; TrivariateFunction f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.cos(omega * z - kx * x - ky * y); } }; double[][][] fval = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { fval[i][j][k] = f.value(xval[i], yval[j], zval[k]); } } } TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator(); TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval); double x, y, z; double expected, result; x = 4; y = -3; z = 0; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("On sample point", expected, result, 1e-12); x = 4.5; y = -1.5; z = -4.25; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("Half-way between sample points (middle of the patch)", expected, result, 0.1); x = 3.5; y = -3.5; z = -10; expected = f.value(x, y, z); result = p.value(x, y, z); Assert.assertEquals("Half-way between sample points (border of the patch)", expected, result, 0.1); }