BicubicSplineInterpolator implements BivariateGridInterpolator { public BicubicSplineInterpolatingFunction interpolate(final double[] xval, final double[] yval, final double[][] fval) throws NoDataException, DimensionMismatchException, NonMonotonicSequenceException, NumberIsTooSmallException { if (xval.length == 0 || yval.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); final int xLen = xval.length; final int yLen = yval.length; final double[][] fX = new double[yLen][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++) { fX[j][i] = fval[i][j]; } } final SplineInterpolator spInterpolator = new SplineInterpolator(); final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; for (int j = 0; j < yLen; j++) { ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); } final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; for (int i = 0; i < xLen; i++) { xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); } final double[][] dFdX = new double[xLen][yLen]; for (int j = 0; j < yLen; j++) { final UnivariateFunction f = ySplineX[j].derivative(); for (int i = 0; i < xLen; i++) { dFdX[i][j] = f.value(xval[i]); } } final double[][] dFdY = new double[xLen][yLen]; for (int i = 0; i < xLen; i++) { final UnivariateFunction f = xSplineY[i].derivative(); for (int j = 0; j < yLen; j++) { dFdY[i][j] = f.value(yval[j]); } } final double[][] d2FdXdY = new double[xLen][yLen]; 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); d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); } } return new BicubicSplineInterpolatingFunction(xval, yval, fval, dFdX, dFdY, d2FdXdY, initializeDerivatives); } BicubicSplineInterpolator();  BicubicSplineInterpolator(boolean initializeDerivatives); BicubicSplineInterpolatingFunction interpolate(final double[] xval,
                                                          final double[] yval,
                                                          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[xval.length][yval.length]; BivariateGridInterpolator interpolator = new BicubicSplineInterpolator(); @SuppressWarnings("unused") BivariateFunction p = interpolator.interpolate(xval, yval, zval); double[] wxval = new double[] {3, 2, 5, 6.5}; try { p = interpolator.interpolate(wxval, yval, zval); 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); Assert.fail("an exception should have been thrown"); } catch (MathIllegalArgumentException e) { } double[][] wzval = new double[xval.length][yval.length + 1]; try { p = interpolator.interpolate(xval, yval, wzval); Assert.fail("an exception should have been thrown"); } catch (DimensionMismatchException e) { } wzval = new double[xval.length - 1][yval.length]; try { p = interpolator.interpolate(xval, yval, wzval); Assert.fail("an exception should have been thrown"); } catch (DimensionMismatchException e) { } }
@Test public void testInterpolation1() { final int sz = 21; double[] xval = new double[sz]; double[] yval = new double[sz]; final double delta = 1d / (sz - 1); for (int i = 0; i < sz; i++) { xval[i] = -1 + 15 * i * delta; yval[i] = -20 + 30 * i * delta; } BivariateFunction f = new BivariateFunction() { public double value(double x, double y) { return 2 * x - 3 * y + 5; } }; double[][] zval = new double[xval.length][yval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { zval[i][j] = f.value(xval[i], yval[j]); } } BivariateGridInterpolator interpolator = new BicubicSplineInterpolator(); BivariateFunction p = interpolator.interpolate(xval, yval, zval); double x, y; final RandomGenerator rng = new Well19937c(1234567L); final UniformRealDistribution distX = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); final UniformRealDistribution distY = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); final int numSamples = 50; final double tol = 6; for (int i = 0; i < numSamples; i++) { x = distX.sample(); for (int j = 0; j < numSamples; j++) { y = distY.sample(); Assert.assertEquals(f.value(x, y), p.value(x, y), tol); } } }
@Test public void testInterpolation2() { final int sz = 21; double[] xval = new double[sz]; double[] yval = new double[sz]; final double delta = 1d / (sz - 1); for (int i = 0; i < sz; i++) { xval[i] = -1 + 15 * i * delta; yval[i] = -20 + 30 * i * delta; } BivariateFunction f = new BivariateFunction() { public double value(double x, double y) { return 2 * x * x - 3 * y * y + 4 * x * y - 5; } }; double[][] zval = new double[xval.length][yval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { zval[i][j] = f.value(xval[i], yval[j]); } } BivariateGridInterpolator interpolator = new BicubicSplineInterpolator(); BivariateFunction p = interpolator.interpolate(xval, yval, zval); double x, y; final RandomGenerator rng = new Well19937c(1234567L); final UniformRealDistribution distX = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); final UniformRealDistribution distY = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); final int numSamples = 50; final double tol = 251; for (int i = 0; i < numSamples; i++) { x = distX.sample(); for (int j = 0; j < numSamples; j++) { y = distY.sample(); Assert.assertEquals(f.value(x, y), p.value(x, y), tol); } } }