(size);
+ for (FacetData fd : other.microsphereData) {
+ microsphereData.add(new FacetData(fd.illumination(), fd.sample()));
+ }
+ }
+
+ /**
+ * Perform a copy.
+ *
+ * @return a copy of this instance.
+ */
+ public InterpolatingMicrosphere copy() {
+ return new InterpolatingMicrosphere(this);
+ }
+
+ /**
+ * Get the space dimensionality.
+ *
+ * @return the number of space dimensions.
+ */
+ public int getDimension() {
+ return dimension;
+ }
+
+ /**
+ * Get the size of the sphere.
+ *
+ * @return the number of surface elements of the microspshere.
+ */
+ public int getSize() {
+ return size;
+ }
+
+ /**
+ * Estimate the value at the requested location.
+ * This microsphere is placed at the given {@code point}, contribution
+ * of the given {@code samplePoints} to each sphere facet is computed
+ * (illumination) and the interpolation is performed (integration of
+ * the illumination).
+ *
+ * @param point Interpolation point.
+ * @param samplePoints Sampling data points.
+ * @param sampleValues Sampling data values at the corresponding
+ * {@code samplePoints}.
+ * @param exponent Exponent used in the power law that computes
+ * the weights (distance dimming factor) of the sample data.
+ * @param noInterpolationTolerance When the distance between the
+ * {@code point} and one of the {@code samplePoints} is less than
+ * this value, no interpolation will be performed, and the value
+ * of the sample will just be returned.
+ * @return the estimated value at the given {@code point}.
+ * @throws NotPositiveException if {@code exponent < 0}.
+ */
+ public double value(double[] point,
+ double[][] samplePoints,
+ double[] sampleValues,
+ double exponent,
+ double noInterpolationTolerance) {
+ if (exponent < 0) {
+ throw new NotPositiveException(exponent);
+ }
+
+ clear();
+
+ // Contribution of each sample point to the illumination of the
+ // microsphere's facets.
+ final int numSamples = samplePoints.length;
+ for (int i = 0; i < numSamples; i++) {
+ // Vector between interpolation point and current sample point.
+ final double[] diff = MathArrays.ebeSubtract(samplePoints[i], point);
+ final double diffNorm = MathArrays.safeNorm(diff);
+
+ if (FastMath.abs(diffNorm) < noInterpolationTolerance) {
+ // No need to interpolate, as the interpolation point is
+ // actually (very close to) one of the sampled points.
+ return sampleValues[i];
+ }
+
+ final double weight = FastMath.pow(diffNorm, -exponent);
+ illuminate(diff, sampleValues[i], weight);
+ }
+
+ return interpolate();
+ }
+
+ /**
+ * Replace {@code i}-th facet of the microsphere.
+ * Method for initializing the microsphere facets.
+ *
+ * @param normal Facet's normal vector.
+ * @param copy Whether to copy the given array.
+ * @throws DimensionMismatchException if the length of {@code n}
+ * does not match the space dimension.
+ * @throws MaxCountExceededException if the method has been called
+ * more times than the size of the sphere.
+ */
+ protected void add(double[] normal,
+ boolean copy) {
+ if (microsphere.size() >= size) {
+ throw new MaxCountExceededException(size);
+ }
+ if (normal.length > dimension) {
+ throw new DimensionMismatchException(normal.length, dimension);
+ }
+
+ microsphere.add(new Facet(copy ? normal.clone() : normal));
+ microsphereData.add(new FacetData(0d, 0d));
+ }
+
+ /**
+ * Interpolation.
+ *
+ * @return the value estimated from the current illumination of the
+ * microsphere.
+ */
+ private double interpolate() {
+ // Number of non-illuminated facets.
+ int darkCount = 0;
+
+ double value = 0;
+ double totalWeight = 0;
+ for (FacetData fd : microsphereData) {
+ final double iV = fd.illumination();
+ if (iV != 0d) {
+ value += iV * fd.sample();
+ totalWeight += iV;
+ } else {
+ ++darkCount;
+ }
+ }
+
+ final double darkFraction = darkCount / (double) size;
+
+ return darkFraction <= maxDarkFraction ?
+ value / totalWeight :
+ background;
+ }
+
+ /**
+ * Illumination.
+ *
+ * @param sampleDirection Vector whose origin is at the interpolation
+ * point and tail is at the sample location.
+ * @param sampleValue Data value of the sample.
+ * @param weight Weight.
+ */
+ private void illuminate(double[] sampleDirection,
+ double sampleValue,
+ double weight) {
+ for (int i = 0; i < size; i++) {
+ final double[] n = microsphere.get(i).getNormal();
+ final double cos = MathArrays.cosAngle(n, sampleDirection);
+
+ if (cos > 0) {
+ final double illumination = cos * weight;
+
+ if (illumination > darkThreshold &&
+ illumination > microsphereData.get(i).illumination()) {
+ microsphereData.set(i, new FacetData(illumination, sampleValue));
+ }
+ }
+ }
+ }
+
+ /**
+ * Reset the all the {@link Facet facets} data to zero.
+ */
+ private void clear() {
+ for (int i = 0; i < size; i++) {
+ microsphereData.set(i, new FacetData(0d, 0d));
+ }
+ }
+
+ /**
+ * Microsphere "facet" (surface element).
+ */
+ private static class Facet {
+ /** Normal vector characterizing a surface element. */
+ private final double[] normal;
+
+ /**
+ * @param n Normal vector characterizing a surface element
+ * of the microsphere. No copy is made.
+ */
+ public Facet(double[] n) {
+ normal = n;
+ }
+
+ /**
+ * Return a reference to the vector normal to this facet.
+ *
+ * @return the normal vector.
+ */
+ public double[] getNormal() {
+ return normal;
+ }
+ }
+
+ /**
+ * Data associated with each {@link Facet}.
+ */
+ private static class FacetData {
+ /** Illumination received from the sample. */
+ private final double illumination;
+ /** Data value of the sample. */
+ private final double sample;
+
+ /**
+ * @param illumination Illumination.
+ * @param sample Data value.
+ */
+ public FacetData(double illumination,
+ double sample) {
+ this.illumination = illumination;
+ this.sample = sample;
+ }
+
+ /**
+ * Get the illumination.
+ * @return the illumination.
+ */
+ public double illumination() {
+ return illumination;
+ }
+
+ /**
+ * Get the data value.
+ * @return the data value.
+ */
+ public double sample() {
+ return sample;
+ }
+ }
+}
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java
new file mode 100644
index 000000000..fdc01b2f3
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java
@@ -0,0 +1,87 @@
+/*
+ * 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.util.FastMath;
+import org.apache.commons.math3.util.MathUtils;
+
+/**
+ * Utility class for the {@link MicrosphereProjectionInterpolator} algorithm.
+ * For 2D interpolation, this class constructs the microsphere as a series of
+ * evenly spaced facets (rather than generating random normals as in the
+ * base implementation).
+ *
+ * @since 3.6
+ */
+public class InterpolatingMicrosphere2D extends InterpolatingMicrosphere {
+ /** Space dimension. */
+ private static final int DIMENSION = 2;
+
+ /**
+ * Create a sphere from vectors regularly sampled around a circle.
+ *
+ * @param size Number of surface elements of the sphere.
+ * @param maxDarkFraction Maximum fraction of the facets that can be dark.
+ * If the fraction of "non-illuminated" facets is larger, no estimation
+ * of the value will be performed, and the {@code background} value will
+ * be returned instead.
+ * @param darkThreshold Value of the illumination below which a facet is
+ * considered dark.
+ * @param background Value returned when the {@code maxDarkFraction}
+ * threshold is exceeded.
+ * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
+ * if {@code size <= 0}.
+ * @throws org.apache.commons.math3.exception.NotPositiveException if
+ * {@code darkThreshold < 0}.
+ * @throws org.apache.commons.math3.exception.OutOfRangeException if
+ * {@code maxDarkFraction} does not belong to the interval {@code [0, 1]}.
+ */
+ public InterpolatingMicrosphere2D(int size,
+ double maxDarkFraction,
+ double darkThreshold,
+ double background) {
+ super(DIMENSION, size, maxDarkFraction, darkThreshold, background);
+
+ // Generate the microsphere normals.
+ for (int i = 0; i < size; i++) {
+ final double angle = i * MathUtils.TWO_PI / size;
+
+ add(new double[] { FastMath.cos(angle),
+ FastMath.sin(angle) },
+ false);
+ }
+ }
+
+ /**
+ * Copy constructor.
+ *
+ * @param other Instance to copy.
+ */
+ protected InterpolatingMicrosphere2D(InterpolatingMicrosphere2D other) {
+ super(other);
+ }
+
+ /**
+ * Perform a copy.
+ *
+ * @return a copy of this instance.
+ */
+ @Override
+ public InterpolatingMicrosphere2D copy() {
+ return new InterpolatingMicrosphere2D(this);
+ }
+}
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java
index b747841f9..58be772e8 100644
--- a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java
@@ -34,7 +34,10 @@ import org.apache.commons.math3.util.FastMath;
* Interpolating function that implements the
* Microsphere Projection.
*
+ * @deprecated Code will be removed in 4.0. Use {@link InterpolatingMicrosphere}
+ * and {@link MicrosphereProjectionInterpolator} instead.
*/
+@Deprecated
public class MicrosphereInterpolatingFunction
implements MultivariateFunction {
/**
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java
index c9881ce99..d9174bc5e 100644
--- a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java
@@ -30,7 +30,10 @@ import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
* MS thesis.
*
* @since 2.1
+ * @deprecated Code will be removed in 4.0. Use {@link InterpolatingMicrosphere}
+ * and {@link MicrosphereProjectionInterpolator} instead.
*/
+@Deprecated
public class MicrosphereInterpolator
implements MultivariateInterpolator {
/**
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java
new file mode 100644
index 000000000..700b2906f
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java
@@ -0,0 +1,166 @@
+/*
+ * 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.analysis.MultivariateFunction;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.NoDataException;
+import org.apache.commons.math3.exception.NotPositiveException;
+import org.apache.commons.math3.exception.NullArgumentException;
+import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
+
+/**
+ * Interpolator that implements the algorithm described in
+ * William Dudziak's
+ * MS thesis.
+ *
+ * @since 3.6
+ */
+public class MicrosphereProjectionInterpolator
+ implements MultivariateInterpolator {
+ /** Brightness exponent. */
+ private final double exponent;
+ /** Microsphere. */
+ private final InterpolatingMicrosphere microsphere;
+ /** Whether to share the sphere. */
+ private final boolean sharedSphere;
+ /** Tolerance value below which no interpolation is necessary. */
+ private final double noInterpolationTolerance;
+
+ /**
+ * Create a microsphere interpolator.
+ *
+ * @param dimension Space dimension.
+ * @param elements Number of surface elements of the microsphere.
+ * @param exponent Exponent used in the power law that computes the
+ * @param maxDarkFraction Maximum fraction of the facets that can be dark.
+ * If the fraction of "non-illuminated" facets is larger, no estimation
+ * of the value will be performed, and the {@code background} value will
+ * be returned instead.
+ * @param darkThreshold Value of the illumination below which a facet is
+ * considered dark.
+ * @param background Value returned when the {@code maxDarkFraction}
+ * threshold is exceeded.
+ * @param sharedSphere Whether the sphere can be shared among the
+ * interpolating function instances. If {@code true}, the instances
+ * will share the same data, and thus will not be thread-safe.
+ * @param noInterpolationTolerance When the distance between an
+ * interpolated point and one of the sample points is less than this
+ * value, no interpolation will be performed (the value of the sample
+ * will be returned).
+ * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
+ * if {@code dimension <= 0} or {@code elements <= 0}.
+ * @throws NotPositiveException if {@code exponent < 0}.
+ * @throws NotPositiveException if {@code darkThreshold < 0}.
+ * @throws org.apache.commons.math3.exception.OutOfRangeException if
+ * {@code maxDarkFraction} does not belong to the interval {@code [0, 1]}.
+ */
+ public MicrosphereProjectionInterpolator(int dimension,
+ int elements,
+ double maxDarkFraction,
+ double darkThreshold,
+ double background,
+ double exponent,
+ boolean sharedSphere,
+ double noInterpolationTolerance) {
+ this(new InterpolatingMicrosphere(dimension,
+ elements,
+ maxDarkFraction,
+ darkThreshold,
+ background,
+ new UnitSphereRandomVectorGenerator(dimension)),
+ exponent,
+ sharedSphere,
+ noInterpolationTolerance);
+ }
+
+ /**
+ * Create a microsphere interpolator.
+ *
+ * @param microsphere Microsphere.
+ * @param exponent Exponent used in the power law that computes the
+ * weights (distance dimming factor) of the sample data.
+ * @param sharedSphere Whether the sphere can be shared among the
+ * interpolating function instances. If {@code true}, the instances
+ * will share the same data, and thus will not be thread-safe.
+ * @param noInterpolationTolerance When the distance between an
+ * interpolated point and one of the sample points is less than this
+ * value, no interpolation will be performed (the value of the sample
+ * will be returned).
+ * @throws NotPositiveException if {@code exponent < 0}.
+ */
+ public MicrosphereProjectionInterpolator(InterpolatingMicrosphere microsphere,
+ double exponent,
+ boolean sharedSphere,
+ double noInterpolationTolerance)
+ throws NotPositiveException {
+ if (exponent < 0) {
+ throw new NotPositiveException(exponent);
+ }
+
+ this.microsphere = microsphere;
+ this.exponent = exponent;
+ this.sharedSphere = sharedSphere;
+ this.noInterpolationTolerance = noInterpolationTolerance;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * @throws DimensionMismatchException if the space dimension of the
+ * given samples does not match the space dimension of the microsphere.
+ */
+ @Override
+ public MultivariateFunction interpolate(final double[][] xval,
+ final double[] yval)
+ throws DimensionMismatchException,
+ NoDataException,
+ NullArgumentException {
+ if (xval == null ||
+ yval == null) {
+ throw new NullArgumentException();
+ }
+ if (xval.length == 0) {
+ throw new NoDataException();
+ }
+ if (xval.length != yval.length) {
+ throw new DimensionMismatchException(xval.length, yval.length);
+ }
+ if (xval[0] == null) {
+ throw new NullArgumentException();
+ }
+ final int dimension = microsphere.getDimension();
+ if (dimension != xval[0].length) {
+ throw new DimensionMismatchException(xval[0].length, dimension);
+ }
+
+ // Microsphere copy.
+ final InterpolatingMicrosphere m = sharedSphere ? microsphere : microsphere.copy();
+
+ return new MultivariateFunction() {
+ /** {inheritDoc} */
+ @Override
+ public double value(double[] point) {
+ return m.value(point,
+ xval,
+ yval,
+ exponent,
+ noInterpolationTolerance);
+ }
+ };
+ }
+}
diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolatorTest.java
new file mode 100644
index 000000000..9eeaba4ba
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolatorTest.java
@@ -0,0 +1,149 @@
+/*
+ * 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.analysis.MultivariateFunction;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test case for the {@link MicrosphereProjectionInterpolator
+ * "microsphere projection"} interpolator.
+ */
+public final class MicrosphereProjectionInterpolatorTest {
+ /**
+ * Test of interpolator for a plane.
+ *
+ * y = 2 x1 - 3 x2 + 5
+ */
+ @Test
+ public void testLinearFunction2D() {
+ MultivariateFunction f = new MultivariateFunction() {
+ @Override
+ public double value(double[] x) {
+ if (x.length != 2) {
+ throw new IllegalArgumentException();
+ }
+ return 2 * x[0] - 3 * x[1] + 5;
+ }
+ };
+
+ final double darkFraction = 0.5;
+ final double darkThreshold = 1e-2;
+ final double background = Double.NaN;
+ final double exponent = 1.1;
+ final boolean shareSphere = true;
+ final double noInterpolationTolerance = Math.ulp(1d);
+
+ // N-dimensional interpolator.
+ final MultivariateInterpolator interpolator
+ = new MicrosphereProjectionInterpolator(2, 500,
+ darkFraction,
+ darkThreshold,
+ background,
+ exponent,
+ shareSphere,
+ noInterpolationTolerance);
+
+ // 2D interpolator.
+ final MultivariateInterpolator interpolator2D
+ = new MicrosphereProjectionInterpolator(new InterpolatingMicrosphere2D(16,
+ darkFraction,
+ darkThreshold,
+ background),
+ exponent,
+ shareSphere,
+ noInterpolationTolerance);
+
+ final double min = -1;
+ final double max = 1;
+ final double range = max - min;
+ final int res = 5;
+ final int n = res * res; // Number of sample points.
+ final int dim = 2;
+ double[][] x = new double[n][dim];
+ double[] y = new double[n];
+ int index = 0;
+ for (int i = 0; i < res; i++) {
+ final double x1Val = toCoordinate(min, range, res, i);
+ for (int j = 0; j < res; j++) {
+ final double x2Val = toCoordinate(min, range, res, j);
+ x[index][0] = x1Val;
+ x[index][1] = x2Val;
+ y[index] = f.value(x[index]);
+ ++index;
+ }
+ }
+
+ final MultivariateFunction p = interpolator.interpolate(x, y);
+ final MultivariateFunction p2D = interpolator2D.interpolate(x, y);
+
+ double[] c = new double[dim];
+ double expected, result, result2D;
+
+ final int sampleIndex = 2;
+ c[0] = x[sampleIndex][0];
+ c[1] = x[sampleIndex][1];
+ expected = f.value(c);
+ result = p.value(c);
+ result2D = p2D.value(c);
+ Assert.assertEquals("on sample point (exact)", expected, result2D, FastMath.ulp(1d));
+ Assert.assertEquals("on sample point (ND vs 2D)", result2D, result, FastMath.ulp(1d));
+
+ // Interpolation.
+ c[0] = 0.654321;
+ c[1] = -0.345678;
+ expected = f.value(c);
+ result = p.value(c);
+ result2D = p2D.value(c);
+ Assert.assertEquals("interpolation (exact)", expected, result2D, 1e-1);
+ Assert.assertEquals("interpolation (ND vs 2D)", result2D, result, 1e-1);
+
+ // Extrapolation.
+ c[0] = 0 - 1e-2;
+ c[1] = 1 + 1e-2;
+ expected = f.value(c);
+ result = p.value(c);
+ result2D = p2D.value(c);
+ Assert.assertFalse(Double.isNaN(result));
+ Assert.assertFalse(Double.isNaN(result2D));
+ Assert.assertEquals("extrapolation (exact)", expected, result2D, 1e-1);
+ Assert.assertEquals("extrapolation (ND vs 2D)", result2D, result, 1e-2);
+
+ // Far away.
+ c[0] = 20;
+ c[1] = -30;
+ result = p.value(c);
+ Assert.assertTrue(result + " should be NaN", Double.isNaN(result));
+ result2D = p2D.value(c);
+ Assert.assertTrue(result2D + " should be NaN", Double.isNaN(result2D));
+ }
+
+ /**
+ * @param min Minimum of the coordinate range.
+ * @param range Extent of the coordinate interval.
+ * @param res Number of pixels.
+ * @param pixel Pixel index.
+ */
+ private static double toCoordinate(double min,
+ double range,
+ int res,
+ int pixel) {
+ return pixel * range / (res - 1) + min;
+ }
+}