You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2014/12/14 18:27:50 UTC
[math] MATH-1173: Tricubic interpolation
Repository: commons-math
Updated Branches:
refs/heads/master 09129d536 -> 753f278d1
MATH-1173: Tricubic interpolation
New class "TricubicInterpolator" to replace "TricubicSplineInterpolator".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/753f278d
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/753f278d
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/753f278d
Branch: refs/heads/master
Commit: 753f278d10fc0c92e912965379401ac182644101
Parents: 09129d5
Author: Gilles <er...@apache.org>
Authored: Sun Dec 14 18:25:08 2014 +0100
Committer: Gilles <er...@apache.org>
Committed: Sun Dec 14 18:25:08 2014 +0100
----------------------------------------------------------------------
src/changes/changes.xml | 11 +-
.../interpolation/TricubicInterpolator.java | 142 ++++++++++++
.../interpolation/TricubicInterpolatorTest.java | 225 +++++++++++++++++++
3 files changed, 375 insertions(+), 3 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/753f278d/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 22a5143..e54d71e 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -59,9 +59,10 @@ Most notable among the new features are:
new distributions (Gumbel, Laplace, Logistic, Nakagami), and
improvements on percentiles algorithms (better handling for NaNs
in the regular algorithm, plus a new storeless implementation).
- Bicubic splines interpolators have been fixed and new implementations
- added. There have been numerous bug fixes and several improvements
- on performances or robustness. See below for a full list)
+ Bicubic and tricubic interpolators have been fixed and new
+ implementations added. There have been numerous bug fixes and
+ several improvements on performances or robustness. See below
+ for a full list.
The minimum version of the Java platform required to compile and use
Apache Commons Math is Java 5.
@@ -78,6 +79,10 @@ Users are encouraged to upgrade to this version as this release not
2. A few methods in the FastMath class are in fact slower that their
counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
">
+ <action dev="erans" type="add" issue="MATH-1173">
+ New classes "TricubicInterpolatingFunction" and "TricubicInterpolator" to
+ replace "TricubicSplineInterpolatingFunction" and "TricubicSplineInterpolator".
+ </action>
<action dev="erans" type="fix" issue="MATH-1178" due-to="Dmitriy">
Fixed example in userguide ("stat" section).
</action>
http://git-wip-us.apache.org/repos/asf/commons-math/blob/753f278d/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java
new file mode 100644
index 0000000..531e736
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java
@@ -0,0 +1,142 @@
+/*
+ * 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 3.4
+ */
+public class TricubicInterpolator
+ implements TrivariateGridInterpolator {
+ /**
+ * {@inheritDoc}
+ */
+ public TricubicInterpolatingFunction 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;
+
+ // Approximation to the partial derivatives using finite differences.
+ final double[][][] dFdX = new double[xLen][yLen][zLen];
+ final double[][][] dFdY = new double[xLen][yLen][zLen];
+ final double[][][] dFdZ = new double[xLen][yLen][zLen];
+ final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
+ final double[][][] d2FdXdZ = new double[xLen][yLen][zLen];
+ final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
+ final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
+
+ for (int i = 1; i < xLen - 1; i++) {
+ if (yval.length != fval[i].length) {
+ throw new DimensionMismatchException(yval.length, fval[i].length);
+ }
+
+ 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++) {
+ if (zval.length != fval[i][j].length) {
+ throw new DimensionMismatchException(zval.length, fval[i][j].length);
+ }
+
+ 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;
+ final double deltaXY = deltaX * deltaY;
+
+ for (int k = 1; k < zLen - 1; k++) {
+ final int nK = k + 1;
+ final int pK = k - 1;
+
+ final double nZ = zval[nK];
+ final double pZ = zval[pK];
+
+ final double deltaZ = nZ - pZ;
+
+ dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX;
+ dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY;
+ dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ;
+
+ final double deltaXZ = deltaX * deltaZ;
+ final double deltaYZ = deltaY * deltaZ;
+
+ d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY;
+ d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ;
+ d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ;
+
+ final double deltaXYZ = deltaXY * deltaZ;
+
+ 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]) / deltaXYZ;
+ }
+ }
+ }
+
+ // Create the interpolating function.
+ return new TricubicInterpolatingFunction(xval, yval, zval, fval,
+ dFdX, dFdY, dFdZ,
+ d2FdXdY, d2FdXdZ, d2FdYdZ,
+ d3FdXdYdZ) {
+ @Override
+ public boolean isValidPoint(double x, double y, double z) {
+ if (x < xval[1] ||
+ x > xval[xval.length - 2] ||
+ y < yval[1] ||
+ y > yval[yval.length - 2] ||
+ z < zval[1] ||
+ z > zval[zval.length - 2]) {
+ return false;
+ } else {
+ return true;
+ }
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/753f278d/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatorTest.java
new file mode 100644
index 0000000..d66e2d6
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatorTest.java
@@ -0,0 +1,225 @@
+/*
+ * 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.util.FastMath;
+import org.apache.commons.math3.analysis.TrivariateFunction;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test case for the {@link TricubicInterpolator tricubic interpolator}.
+ */
+public final class TricubicInterpolatorTest {
+ /**
+ * 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[] {-12, -8, -5.5, -3, 0, 2.5};
+ double[][][] fval = new double[xval.length][yval.length][zval.length];
+
+ @SuppressWarnings("unused")
+ TrivariateFunction tcf = new TricubicInterpolator().interpolate(xval, yval, zval, fval);
+
+ double[] wxval = new double[] {3, 2, 5, 6.5};
+ try {
+ tcf = new TricubicInterpolator().interpolate(wxval, yval, zval, fval);
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException e) {
+ // Expected
+ }
+ double[] wyval = new double[] {-4, -1, -1, 2.5};
+ try {
+ tcf = new TricubicInterpolator().interpolate(xval, wyval, zval, fval);
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException e) {
+ // Expected
+ }
+ double[] wzval = new double[] {-12, -8, -9, -3, 0, 2.5};
+ try {
+ tcf = new TricubicInterpolator().interpolate(xval, yval, wzval, fval);
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException e) {
+ // Expected
+ }
+ double[][][] wfval = new double[xval.length - 1][yval.length][zval.length];
+ try {
+ tcf = new TricubicInterpolator().interpolate(xval, yval, zval, wfval);
+ Assert.fail("an exception should have been thrown");
+ } catch (DimensionMismatchException e) {
+ // Expected
+ }
+ wfval = new double[xval.length][yval.length - 1][zval.length];
+ try {
+ tcf = new TricubicInterpolator().interpolate(xval, yval, zval, wfval);
+ Assert.fail("an exception should have been thrown");
+ } catch (DimensionMismatchException e) {
+ // Expected
+ }
+ wfval = new double[xval.length][yval.length][zval.length - 1];
+ try {
+ tcf = new TricubicInterpolator().interpolate(xval, yval, zval, wfval);
+ Assert.fail("an exception should have been thrown");
+ } catch (DimensionMismatchException e) {
+ // Expected
+ }
+ }
+
+ public void testIsValid() {
+ 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];
+
+ @SuppressWarnings("unused")
+ TricubicInterpolatingFunction tcf = new TricubicInterpolator().interpolate(xval, yval, zval, fval);
+
+ // Valid.
+ Assert.assertTrue(tcf.isValidPoint(4, -3, -8));
+ Assert.assertTrue(tcf.isValidPoint(5, -3, -8));
+ Assert.assertTrue(tcf.isValidPoint(4, -1, -8));
+ Assert.assertTrue(tcf.isValidPoint(5, -1, -8));
+ Assert.assertTrue(tcf.isValidPoint(4, -3, 0));
+ Assert.assertTrue(tcf.isValidPoint(5, -3, 0));
+ Assert.assertTrue(tcf.isValidPoint(4, -1, 0));
+ Assert.assertTrue(tcf.isValidPoint(5, -1, 0));
+
+ // Invalid.
+ Assert.assertFalse(tcf.isValidPoint(3.5, -3, -8));
+ Assert.assertFalse(tcf.isValidPoint(4.5, -3.1, -8));
+ Assert.assertFalse(tcf.isValidPoint(4.5, -2, 0));
+ Assert.assertFalse(tcf.isValidPoint(4.5, 0, -3.5));
+ Assert.assertFalse(tcf.isValidPoint(-10, 4.1, -1));
+ }
+
+ /**
+ * Test for a plane.
+ * <p>
+ * f(x, y, z) = 2 x - 3 y - 4 z + 5
+ * </p>
+ */
+ @Test
+ public void testPlane() {
+ 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};
+
+ // Function values
+ TrivariateFunction f = new TrivariateFunction() {
+ public double value(double x, double y, double z) {
+ return 2 * x - 3 * y - 4 * z + 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 tcf = new TricubicInterpolator().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 = tcf.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 = tcf.value(x, y, z);
+ Assert.assertEquals("Half-way between sample points (middle of the patch)",
+ expected, result, 1e-14);
+ }
+
+ /**
+ * Sine wave.
+ * <p>
+ * f(x, y, z) = a cos [ω z - k<sub>y</sub> x - k<sub>y</sub> y]
+ * </p>
+ * with A = 0.2, ω = 0.5, k<sub>x</sub> = 2, k<sub>y</sub> = 1.
+ */
+ @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;
+
+ // Function values
+ 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]);
+ }
+ }
+ }
+
+ TrivariateFunction tcf = new TricubicInterpolator().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 = tcf.value(x, y, z);
+ Assert.assertEquals("On sample point",
+ expected, result, 1e-14);
+
+ x = 4.5;
+ y = -1.5;
+ z = -4.25;
+ expected = f.value(x, y, z);
+ result = tcf.value(x, y, z);
+ Assert.assertEquals("Half-way between sample points (middle of the patch)",
+ expected, result, 1e-1); // XXX Too high tolerance!
+ }
+}