commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [1/4] [math] MATH-1166: Bicubic interpolation
Date Tue, 02 Dec 2014 11:06:54 GMT
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 <erans@apache.org>
Authored: Fri Nov 28 15:16:51 2014 +0100
Committer: Gilles <erans@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);
+        }
+    }
+}


Mime
View raw message