MATH-1243

Refactoring of "microsphere interpolation" algorithm.
This commit is contained in:
Gilles 2015-07-06 13:52:02 +02:00
parent ed6c47dd9e
commit 5f9dda6b90
5 changed files with 796 additions and 0 deletions

View File

@ -54,6 +54,14 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
<action dev="erans" type="update" issue="MATH-1243">
Refactored implementation of the "miscrosphere projection"
interpolation algorithm.
New classes: "MicrosphereProjectionInterpolator",
"InterpolatingMicrosphere" and "InterpolatingMicrosphere2D"
replace "MicrosphereInterpolator" and "MicrosphereInterpolatingFunction".
(package "o.a.c.m.analysis.interpolation").
</action>
<action dev="erans" type="add" issue="MATH-1244">
Method "cosAngle" in "o.a.c.m.util.MathArrays".
</action>

View File

@ -0,0 +1,386 @@
/*
* 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.math4.analysis.interpolation;
import java.util.List;
import java.util.ArrayList;
import org.apache.commons.math4.random.UnitSphereRandomVectorGenerator;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NotPositiveException;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.MaxCountExceededException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathArrays;
/**
* Utility class for the {@link MicrosphereProjectionInterpolator} algorithm.
*
* @since 4.0
*/
public class InterpolatingMicrosphere {
/** Microsphere. */
private final List<Facet> microsphere;
/** Microsphere data. */
private final List<FacetData> microsphereData;
/** Space dimension. */
private final int dimension;
/** Number of surface elements. */
private final int size;
/** Maximum fraction of the facets that can be dark. */
private final double maxDarkFraction;
/** Lowest non-zero illumination. */
private final double darkThreshold;
/** Background value. */
private final double background;
/**
* Create an unitialiazed sphere.
* Sub-classes are responsible for calling the {@code add(double[]) add}
* method in order to initialize all the sphere's facets.
*
* @param dimension Dimension of the data space.
* @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 NotStrictlyPositiveException if {@code dimension <= 0}
* or {@code size <= 0}.
* @throws NotPositiveException if {@code darkThreshold < 0}.
* @throws OutOfRangeException if {@code maxDarkFraction} does not
* belong to the interval {@code [0, 1]}.
*/
protected InterpolatingMicrosphere(int dimension,
int size,
double maxDarkFraction,
double darkThreshold,
double background) {
if (dimension <= 0) {
throw new NotStrictlyPositiveException(dimension);
}
if (size <= 0) {
throw new NotStrictlyPositiveException(size);
}
if (maxDarkFraction < 0 ||
maxDarkFraction > 1) {
throw new OutOfRangeException(maxDarkFraction, 0, 1);
}
if (darkThreshold < 0) {
throw new NotPositiveException(darkThreshold);
}
this.dimension = dimension;
this.size = size;
this.maxDarkFraction = maxDarkFraction;
this.darkThreshold = darkThreshold;
this.background = background;
microsphere = new ArrayList<Facet>(size);
microsphereData = new ArrayList<FacetData>(size);
}
/**
* Create a sphere from randomly sampled vectors.
*
* @param dimension Dimension of the data space.
* @param size Number of surface elements of the sphere.
* @param rand Unit vector generator for creating the microsphere.
* @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 DimensionMismatchException if the size of the generated
* vectors does not match the dimension set in the constructor.
* @throws NotStrictlyPositiveException if {@code dimension <= 0}
* or {@code size <= 0}.
* @throws NotPositiveException if {@code darkThreshold < 0}.
* @throws OutOfRangeException if {@code maxDarkFraction} does not
* belong to the interval {@code [0, 1]}.
*/
public InterpolatingMicrosphere(int dimension,
int size,
double maxDarkFraction,
double darkThreshold,
double background,
UnitSphereRandomVectorGenerator rand) {
this(dimension, size, maxDarkFraction, darkThreshold, background);
// Generate the microsphere normals, assuming that a number of
// randomly generated normals will represent a sphere.
for (int i = 0; i < size; i++) {
add(rand.nextVector(), false);
}
}
/**
* Copy constructor.
*
* @param other Instance to copy.
*/
protected InterpolatingMicrosphere(InterpolatingMicrosphere other) {
dimension = other.dimension;
size = other.size;
maxDarkFraction = other.maxDarkFraction;
darkThreshold = other.darkThreshold;
background = other.background;
// Field can be shared.
microsphere = other.microsphere;
// Field must be copied.
microsphereData = new ArrayList<FacetData>(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 {@link Facet 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;
}
}
}

View File

@ -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.math4.analysis.interpolation;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.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 4.0
*/
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.math4.exception.NotStrictlyPositiveException
* if {@code size <= 0}.
* @throws org.apache.commons.math4.exception.NotPositiveException if
* {@code darkThreshold < 0}.
* @throws org.apache.commons.math4.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);
}
}

View File

@ -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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.MultivariateFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NotPositiveException;
import org.apache.commons.math4.exception.NullArgumentException;
import org.apache.commons.math4.random.UnitSphereRandomVectorGenerator;
/**
* Interpolator that implements the algorithm described in
* <em>William Dudziak</em>'s
* <a href="http://www.dudziak.com/microsphere.pdf">MS thesis</a>.
*
* @since 4.0
*/
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 <em>not</em> 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.math4.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.math4.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 <em>not</em> 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);
}
};
}
}

View File

@ -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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.MultivariateFunction;
import org.apache.commons.math4.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.
* <p>
* y = 2 x<sub>1</sub> - 3 x<sub>2</sub> + 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;
}
}