Repository: commons-math Updated Branches: refs/heads/master 0802df0f7 -> e11c00085
MATH-1166: Bicubic interpolation New classes to replace "BicubicSplineInterpolatingFunction" and "BicubicSplineInterpolator". Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/136bf342 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/136bf342 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/136bf342 Branch: refs/heads/master Commit: 136bf34297a9ac8049e451800df819ebf092cdf2 Parents: 8aec465 Author: Gilles <er...@apache.org> Authored: Fri Nov 28 15:16:51 2014 +0100 Committer: Gilles <er...@apache.org> Committed: Fri Nov 28 15:16:51 2014 +0100 ---------------------------------------------------------------------- .../BicubicInterpolatingFunction.java | 325 ++++++++++++++++ .../interpolation/BicubicInterpolator.java | 112 ++++++ .../BicubicInterpolatingFunctionTest.java | 383 +++++++++++++++++++ .../interpolation/BicubicInterpolatorTest.java | 175 +++++++++ 4 files changed, 995 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/136bf342/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java new file mode 100644 index 0000000..7fe947f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java @@ -0,0 +1,325 @@ +/* + * 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 java.util.Arrays; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Function that implements the + * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation";> + * bicubic spline interpolation</a>. + * + * @since 3.4 + */ +public class BicubicInterpolatingFunction + implements BivariateFunction { + /** Number of coefficients. */ + private static final int NUM_COEFF = 16; + /** + * Matrix to compute the spline coefficients from the function values + * and function derivatives values + */ + private static final double[][] AINV = { + { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, + { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, + { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, + { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, + { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, + { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, + { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, + { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, + { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, + { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, + { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } + }; + + /** Samples x-coordinates */ + private final double[] xval; + /** Samples y-coordinates */ + private final double[] yval; + /** Set of cubic splines patching the whole data grid */ + private final BicubicFunction[][] splines; + + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect + * to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect + * to y on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on + * every grid point. + * @throws DimensionMismatchException if the various arrays do not contain + * the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x} or {@code y} are + * not strictly increasing. + * @throws NoDataException if any of the arrays has zero length. + */ + public BicubicInterpolatingFunction(double[] x, + double[] y, + double[][] f, + double[][] dFdX, + double[][] dFdY, + double[][] d2FdXdY) + throws DimensionMismatchException, + NoDataException, + NonMonotonicSequenceException { + final int xLen = x.length; + final int yLen = y.length; + + if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { + throw new NoDataException(); + } + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + if (xLen != dFdX.length) { + throw new DimensionMismatchException(xLen, dFdX.length); + } + if (xLen != dFdY.length) { + throw new DimensionMismatchException(xLen, dFdY.length); + } + if (xLen != d2FdXdY.length) { + throw new DimensionMismatchException(xLen, d2FdXdY.length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + + xval = x.clone(); + yval = y.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + splines = new BicubicFunction[lastI][lastJ]; + + for (int i = 0; i < lastI; i++) { + if (f[i].length != yLen) { + throw new DimensionMismatchException(f[i].length, yLen); + } + if (dFdX[i].length != yLen) { + throw new DimensionMismatchException(dFdX[i].length, yLen); + } + if (dFdY[i].length != yLen) { + throw new DimensionMismatchException(dFdY[i].length, yLen); + } + if (d2FdXdY[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdY[i].length, yLen); + } + final int ip1 = i + 1; + final double xR = xval[ip1] - xval[i]; + for (int j = 0; j < lastJ; j++) { + final int jp1 = j + 1; + final double yR = yval[jp1] - yval[j]; + final double xRyR = xR * yR; + final double[] beta = new double[] { + f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], + dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR, + dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR, + d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR + }; + + splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta)); + } + } + } + + /** + * {@inheritDoc} + */ + public double value(double x, double y) + throws OutOfRangeException { + final int i = searchIndex(x, xval); + final int j = searchIndex(y, yval); + + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + + return splines[i][j].value(xN, yN); + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @return {@code true} if (x, y) is a valid point. + */ + public boolean isValidPoint(double x, double y) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @return the index in {@code val} corresponding to the interval + * containing {@code c}. + * @throws OutOfRangeException if {@code c} is out of the + * range defined by the boundary values of {@code val}. + */ + private int searchIndex(double c, double[] val) { + final int r = Arrays.binarySearch(val, c); + + if (r == -1 || + r == -val.length - 1) { + throw new OutOfRangeException(c, val[0], val[val.length - 1]); + } + + if (r < 0) { + // "c" in within an interpolation sub-interval: Return the + // index of the sample at the lower end of the sub-interval. + return -r - 2; + } + final int last = val.length - 1; + if (r == last) { + // "c" is the last sample of the range: Return the index + // of the sample at the lower end of the last sub-interval. + return last - 1; + } + + // "c" is another sample point. + return r; + } + + /** + * Compute the spline coefficients from the list of function values and + * function partial derivatives values at the four corners of a grid + * element. They must be specified in the following order: + * <ul> + * <li>f(0,0)</li> + * <li>f(1,0)</li> + * <li>f(0,1)</li> + * <li>f(1,1)</li> + * <li>f<sub>x</sub>(0,0)</li> + * <li>f<sub>x</sub>(1,0)</li> + * <li>f<sub>x</sub>(0,1)</li> + * <li>f<sub>x</sub>(1,1)</li> + * <li>f<sub>y</sub>(0,0)</li> + * <li>f<sub>y</sub>(1,0)</li> + * <li>f<sub>y</sub>(0,1)</li> + * <li>f<sub>y</sub>(1,1)</li> + * <li>f<sub>xy</sub>(0,0)</li> + * <li>f<sub>xy</sub>(1,0)</li> + * <li>f<sub>xy</sub>(0,1)</li> + * <li>f<sub>xy</sub>(1,1)</li> + * </ul> + * where the subscripts indicate the partial derivative with respect to + * the corresponding variable(s). + * + * @param beta List of function values and function partial derivatives + * values. + * @return the spline coefficients. + */ + private double[] computeSplineCoefficients(double[] beta) { + final double[] a = new double[NUM_COEFF]; + + for (int i = 0; i < NUM_COEFF; i++) { + double result = 0; + final double[] row = AINV[i]; + for (int j = 0; j < NUM_COEFF; j++) { + result += row[j] * beta[j]; + } + a[i] = result; + } + + return a; + } +} + +/** + * Bicubic function. + */ +class BicubicFunction implements BivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][] a; + + /** + * Simple constructor. + * + * @param coeff Spline coefficients. + */ + public BicubicFunction(double[] coeff) { + a = new double[N][N]; + for (int j = 0; j < N; j++) { + final double[] aJ = a[j]; + for (int i = 0; i < N; i++) { + aJ[i] = coeff[i * N + j]; + } + } + } + + /** + * {@inheritDoc} + */ + public double value(double x, double y) { + if (x < 0 || x > 1) { + throw new OutOfRangeException(x, 0, 1); + } + if (y < 0 || y > 1) { + throw new OutOfRangeException(y, 0, 1); + } + + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = {1, x, x2, x3}; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = {1, y, y2, y3}; + + return apply(pX, pY, a); + } + + /** + * Compute the value of the bicubic polynomial. + * + * @param pX Powers of the x-coordinate. + * @param pY Powers of the y-coordinate. + * @param coeff Spline coefficients. + * @return the interpolated value. + */ + private double apply(double[] pX, double[] pY, double[][] coeff) { + double result = 0; + for (int i = 0; i < N; i++) { + final double r = MathArrays.linearCombination(coeff[i], pY); + result += r * pX[i]; + } + + return result; + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/136bf342/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java new file mode 100644 index 0000000..7a3f948 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java @@ -0,0 +1,112 @@ +/* + * 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 {@link BicubicInterpolatingFunction bicubic interpolating + * function}. + * <p> + * Caveat: Because the interpolation scheme requires that derivatives be + * specified at the sample points, those are approximated with finite + * differences (using the 2-points symmetric formulae). + * Since their values are undefined at the borders of the provided + * interpolation ranges, the interpolated values will be wrong at the + * edges of the patch. + * The {@code interpolate} method will return a function that overrides + * {@link isValidPoint(double,double)} to indicate points where the + * interpolation will be inaccurate. + * </p> + * + * @since 3.4 + */ +public class BicubicInterpolator + implements BivariateGridInterpolator { + /** + * {@inheritDoc} + */ + public BicubicInterpolatingFunction 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; + + // Approximation to the partial derivatives using finite differences. + final double[][] dFdX = new double[xLen][yLen]; + final double[][] dFdY = new double[xLen][yLen]; + final double[][] d2FdXdY = new double[xLen][yLen]; + for (int i = 1; i < xLen - 1; i++) { + final int nI = i + 1; + final int pI = i - 1; + + final double nX = xval[nI]; + final double pX = xval[pI]; + + final double deltaX = nX - pX; + + for (int j = 1; j < yLen - 1; j++) { + final int nJ = j + 1; + final int pJ = j - 1; + + final double nY = yval[nJ]; + final double pY = yval[pJ]; + + final double deltaY = nY - pY; + + dFdX[i][j] = (fval[nI][j] - fval[pI][j]) / deltaX; + dFdY[i][j] = (fval[i][nJ] - fval[i][pJ]) / deltaY; + + final double deltaXY = deltaX * deltaY; + + d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / deltaXY; + } + } + + // Create the interpolating function. + return new BicubicInterpolatingFunction(xval, yval, fval, + dFdX, dFdY, d2FdXdY) { + @Override + public boolean isValidPoint(double x, double y) { + if (x < xval[1] || + x > xval[xval.length - 2] || + y < yval[1] || + y > yval[yval.length - 2]) { + return false; + } else { + return true; + } + } + }; + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/136bf342/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunctionTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunctionTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunctionTest.java new file mode 100644 index 0000000..6d56ed8 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunctionTest.java @@ -0,0 +1,383 @@ +/* + * 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.MathIllegalArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; +import org.junit.Assert; +import org.junit.Test; +import org.junit.Ignore; + +/** + * Test case for the bicubic function. + */ +public final class BicubicInterpolatingFunctionTest { + /** + * Test preconditions. + */ + @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]; + + @SuppressWarnings("unused") + BivariateFunction bcf = new BicubicInterpolatingFunction(xval, yval, zval, + zval, zval, zval); + + double[] wxval = new double[] {3, 2, 5, 6.5}; + try { + bcf = new BicubicInterpolatingFunction(wxval, yval, zval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[] wyval = new double[] {-4, -1, -1, 2.5}; + try { + bcf = new BicubicInterpolatingFunction(xval, wyval, zval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[][] wzval = new double[xval.length][yval.length - 1]; + try { + bcf = new BicubicInterpolatingFunction(xval, yval, wzval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, wzval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, wzval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, zval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + + wzval = new double[xval.length - 1][yval.length]; + try { + bcf = new BicubicInterpolatingFunction(xval, yval, wzval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, wzval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, wzval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, zval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + } + + @Test + public void testIsValidPoint() { + final double xMin = -12; + final double xMax = 34; + final double yMin = 5; + final double yMax = 67; + final double[] xval = new double[] { xMin, xMax }; + final double[] yval = new double[] { yMin, yMax }; + final double[][] f = new double[][] { { 1, 2 }, + { 3, 4 } }; + final double[][] dFdX = f; + final double[][] dFdY = f; + final double[][] dFdXdY = f; + + final BicubicInterpolatingFunction bcf + = new BicubicInterpolatingFunction(xval, yval, f, + dFdX, dFdY, dFdXdY); + + double x, y; + + x = xMin; + y = yMin; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + x = xMax; + y = yMax; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + final double xRange = xMax - xMin; + final double yRange = yMax - yMin; + x = xMin + xRange / 3.4; + y = yMin + yRange / 1.2; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + final double small = 1e-8; + x = xMin - small; + y = yMax; + Assert.assertFalse(bcf.isValidPoint(x, y)); + // Ensure that an exception would have been thrown. + try { + bcf.value(x, y); + Assert.fail("OutOfRangeException expected"); + } catch (OutOfRangeException expected) {} + + x = xMin; + y = yMax + small; + Assert.assertFalse(bcf.isValidPoint(x, y)); + // Ensure that an exception would have been thrown. + try { + bcf.value(x, y); + Assert.fail("OutOfRangeException expected"); + } catch (OutOfRangeException expected) {} + } + + /** + * Interpolating a plane. + * <p> + * z = 2 x - 3 y + 5 + */ + @Test + public void testPlane() { + final int numberOfElements = 10; + final double minimumX = -10; + final double maximumX = 10; + final double minimumY = -10; + final double maximumY = 10; + final int numberOfSamples = 1000; + + final double interpolationTolerance = 1e-15; + final double maxTolerance = 1e-14; + + // Function values + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x - 3 * y + 5; + } + }; + BivariateFunction dfdx = new BivariateFunction() { + public double value(double x, double y) { + return 2; + } + }; + BivariateFunction dfdy = new BivariateFunction() { + public double value(double x, double y) { + return -3; + } + }; + BivariateFunction d2fdxdy = new BivariateFunction() { + public double value(double x, double y) { + return 0; + } + }; + + testInterpolation(minimumX, + maximumX, + minimumY, + maximumY, + numberOfElements, + numberOfSamples, + f, + dfdx, + dfdy, + d2fdxdy, + interpolationTolerance, + maxTolerance, + false); + } + + /** + * Interpolating a paraboloid. + * <p> + * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5 + */ + @Test + public void testParaboloid() { + final int numberOfElements = 10; + final double minimumX = -10; + final double maximumX = 10; + final double minimumY = -10; + final double maximumY = 10; + final int numberOfSamples = 1000; + + final double interpolationTolerance = 2e-14; + final double maxTolerance = 1e-12; + + // Function values + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x * x - 3 * y * y + 4 * x * y - 5; + } + }; + BivariateFunction dfdx = new BivariateFunction() { + public double value(double x, double y) { + return 4 * (x + y); + } + }; + BivariateFunction dfdy = new BivariateFunction() { + public double value(double x, double y) { + return 4 * x - 6 * y; + } + }; + BivariateFunction d2fdxdy = new BivariateFunction() { + public double value(double x, double y) { + return 4; + } + }; + + testInterpolation(minimumX, + maximumX, + minimumY, + maximumY, + numberOfElements, + numberOfSamples, + f, + dfdx, + dfdy, + d2fdxdy, + interpolationTolerance, + maxTolerance, + false); + } + + /** + * @param minimumX Lower bound of interpolation range along the x-coordinate. + * @param maximumX Higher bound of interpolation range along the x-coordinate. + * @param minimumY Lower bound of interpolation range along the y-coordinate. + * @param maximumY Higher bound of interpolation range along the y-coordinate. + * @param numberOfElements Number of data points (along each dimension). + * @param numberOfSamples Number of test points. + * @param f Function to test. + * @param dfdx Partial derivative w.r.t. x of the function to test. + * @param dfdy Partial derivative w.r.t. y of the function to test. + * @param d2fdxdy Second partial cross-derivative of the function to test. + * @param meanTolerance Allowed average error (mean error on all interpolated values). + * @param maxTolerance Allowed error on each interpolated value. + */ + private void testInterpolation(double minimumX, + double maximumX, + double minimumY, + double maximumY, + int numberOfElements, + int numberOfSamples, + BivariateFunction f, + BivariateFunction dfdx, + BivariateFunction dfdy, + BivariateFunction d2fdxdy, + double meanTolerance, + double maxTolerance, + boolean print) { + double expected; + double actual; + double currentX; + double currentY; + final double deltaX = (maximumX - minimumX) / numberOfElements; + final double deltaY = (maximumY - minimumY) / numberOfElements; + final double[] xValues = new double[numberOfElements]; + final double[] yValues = new double[numberOfElements]; + final double[][] zValues = new double[numberOfElements][numberOfElements]; + final double[][] dzdx = new double[numberOfElements][numberOfElements]; + final double[][] dzdy = new double[numberOfElements][numberOfElements]; + final double[][] d2zdxdy = new double[numberOfElements][numberOfElements]; + + for (int i = 0; i < numberOfElements; i++) { + xValues[i] = minimumX + deltaX * i; + final double x = xValues[i]; + for (int j = 0; j < numberOfElements; j++) { + yValues[j] = minimumY + deltaY * j; + final double y = yValues[j]; + zValues[i][j] = f.value(x, y); + dzdx[i][j] = dfdx.value(x, y); + dzdy[i][j] = dfdy.value(x, y); + d2zdxdy[i][j] = d2fdxdy.value(x, y); + } + } + + final BivariateFunction interpolation + = new BicubicInterpolatingFunction(xValues, + yValues, + zValues, + dzdx, + dzdy, + d2zdxdy); + + for (int i = 0; i < numberOfElements; i++) { + currentX = xValues[i]; + for (int j = 0; j < numberOfElements; j++) { + currentY = yValues[j]; + expected = f.value(currentX, currentY); + actual = interpolation.value(currentX, currentY); + Assert.assertTrue("On data point: " + expected + " != " + actual, + Precision.equals(expected, actual)); + } + } + + final RandomGenerator rng = new Well19937c(1234567L); + final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]); + final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]); + + double sumError = 0; + for (int i = 0; i < numberOfSamples; i++) { + currentX = distX.sample(); + currentY = distY.sample(); + expected = f.value(currentX, currentY); + + if (print) { + System.out.println(currentX + " " + currentY + " -> "); + } + + actual = interpolation.value(currentX, currentY); + sumError += FastMath.abs(actual - expected); + + if (print) { + System.out.println(actual + " (diff=" + (expected - actual) + ")"); + } + + Assert.assertEquals(expected, actual, maxTolerance); + } + + final double meanError = sumError / numberOfSamples; + Assert.assertEquals(0, meanError, meanTolerance); + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/136bf342/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatorTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatorTest.java new file mode 100644 index 0000000..e2a3550 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatorTest.java @@ -0,0 +1,175 @@ +/* + * 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.MathIllegalArgumentException; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test case for the bicubic interpolator. + */ +public final class BicubicInterpolatorTest { + /** + * Test preconditions. + */ + @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 BicubicInterpolator(); + + @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) { + // Expected + } + + 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) { + // Expected + } + + 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) { + // Expected + } + 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) { + // Expected + } + } + + /** + * Interpolating a plane. + * <p> + * z = 2 x - 3 y + 5 + */ + @Test + public void testPlane() { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x - 3 * y + 5; + } + }; + + testInterpolation(3000, + 1e-13, + f, + false); + } + + /** + * Interpolating a paraboloid. + * <p> + * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5 + */ + @Test + public void testParaboloid() { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x * x - 3 * y * y + 4 * x * y - 5; + } + }; + + testInterpolation(3000, + 1e-12, + f, + false); + } + + /** + * @param numSamples Number of test samples. + * @param tolerance Allowed tolerance on the interpolated value. + * @param f Test function. + * @param print Whether to print debugging output to the console. + */ + private void testInterpolation(int numSamples, + double tolerance, + BivariateFunction f, + boolean print) { + final int sz = 21; + final double[] xval = new double[sz]; + final double[] yval = new double[sz]; + // Coordinate values + 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; + } + + final 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]); + } + } + + final BicubicInterpolator interpolator = new BicubicInterpolator(); + final BicubicInterpolatingFunction p = interpolator.interpolate(xval, yval, zval); + double x, y; + + final RandomGenerator rng = new Well19937c(); + final UniformRealDistribution distX = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); + final UniformRealDistribution distY = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); + + int count = 0; + while (true) { + x = distX.sample(); + y = distY.sample(); + if (!p.isValidPoint(x, y)) { + if (print) { + System.out.println("# " + x + " " + y); + } + continue; + } + + if (count++ > numSamples) { + break; + } + final double expected = f.value(x, y); + final double actual = p.value(x, y); + + if (print) { + System.out.println(x + " " + y + " " + expected + " " + actual); + } + + Assert.assertEquals(expected, actual, tolerance); + } + } +}
