diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java')
-rw-r--r-- | src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java | 244 |
1 files changed, 244 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java b/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java new file mode 100644 index 0000000..a710e82 --- /dev/null +++ b/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java @@ -0,0 +1,244 @@ +/* + * 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.math.analysis.interpolation; + +import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import org.apache.commons.math.DimensionMismatchException; +import org.apache.commons.math.analysis.MultivariateRealFunction; +import org.apache.commons.math.exception.NoDataException; +import org.apache.commons.math.linear.ArrayRealVector; +import org.apache.commons.math.linear.RealVector; +import org.apache.commons.math.random.UnitSphereRandomVectorGenerator; +import org.apache.commons.math.util.FastMath; + +/** + * Interpolating function that implements the + * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>. + * + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + */ +public class MicrosphereInterpolatingFunction + implements MultivariateRealFunction { + /** + * Space dimension. + */ + private final int dimension; + /** + * Internal accounting data for the interpolation algorithm. + * Each element of the list corresponds to one surface element of + * the microsphere. + */ + private final List<MicrosphereSurfaceElement> microsphere; + /** + * Exponent used in the power law that computes the weights of the + * sample data. + */ + private final double brightnessExponent; + /** + * Sample data. + */ + private final Map<RealVector, Double> samples; + + /** + * Class for storing the accounting data needed to perform the + * microsphere projection. + */ + private static class MicrosphereSurfaceElement { + + /** Normal vector characterizing a surface element. */ + private final RealVector normal; + + /** Illumination received from the brightest sample. */ + private double brightestIllumination; + + /** Brightest sample. */ + private Map.Entry<RealVector, Double> brightestSample; + + /** + * @param n Normal vector characterizing a surface element + * of the microsphere. + */ + MicrosphereSurfaceElement(double[] n) { + normal = new ArrayRealVector(n); + } + + /** + * Return the normal vector. + * @return the normal vector + */ + RealVector normal() { + return normal; + } + + /** + * Reset "illumination" and "sampleIndex". + */ + void reset() { + brightestIllumination = 0; + brightestSample = null; + } + + /** + * Store the illumination and index of the brightest sample. + * @param illuminationFromSample illumination received from sample + * @param sample current sample illuminating the element + */ + void store(final double illuminationFromSample, + final Map.Entry<RealVector, Double> sample) { + if (illuminationFromSample > this.brightestIllumination) { + this.brightestIllumination = illuminationFromSample; + this.brightestSample = sample; + } + } + + /** + * Get the illumination of the element. + * @return the illumination. + */ + double illumination() { + return brightestIllumination; + } + + /** + * Get the sample illuminating the element the most. + * @return the sample. + */ + Map.Entry<RealVector, Double> sample() { + return brightestSample; + } + } + + /** + * @param xval the arguments for the interpolation points. + * {@code xval[i][0]} is the first component of interpolation point + * {@code i}, {@code xval[i][1]} is the second component, and so on + * until {@code xval[i][d-1]}, the last component of that interpolation + * point (where {@code dimension} is thus the dimension of the sampled + * space). + * @param yval the values for the interpolation points + * @param brightnessExponent Brightness dimming factor. + * @param microsphereElements Number of surface elements of the + * microsphere. + * @param rand Unit vector generator for creating the microsphere. + * @throws DimensionMismatchException if the lengths of {@code yval} and + * {@code xval} (equal to {@code n}, the number of interpolation points) + * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]}, + * have lengths different from {@code dimension}. + * @throws NoDataException if there are no data (xval null or zero length) + */ + public MicrosphereInterpolatingFunction(double[][] xval, + double[] yval, + int brightnessExponent, + int microsphereElements, + UnitSphereRandomVectorGenerator rand) + throws DimensionMismatchException, NoDataException { + if (xval.length == 0 || xval[0] == null) { + throw new NoDataException(); + } + + if (xval.length != yval.length) { + throw new DimensionMismatchException(xval.length, yval.length); + } + + dimension = xval[0].length; + this.brightnessExponent = brightnessExponent; + + // Copy data samples. + samples = new HashMap<RealVector, Double>(yval.length); + for (int i = 0; i < xval.length; ++i) { + final double[] xvalI = xval[i]; + if ( xvalI.length != dimension) { + throw new DimensionMismatchException(xvalI.length, dimension); + } + + samples.put(new ArrayRealVector(xvalI), yval[i]); + } + + microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements); + // Generate the microsphere, assuming that a fairly large number of + // randomly generated normals will represent a sphere. + for (int i = 0; i < microsphereElements; i++) { + microsphere.add(new MicrosphereSurfaceElement(rand.nextVector())); + } + + } + + /** + * @param point Interpolation point. + * @return the interpolated value. + */ + public double value(double[] point) { + + final RealVector p = new ArrayRealVector(point); + + // Reset. + for (MicrosphereSurfaceElement md : microsphere) { + md.reset(); + } + + // Compute contribution of each sample points to the microsphere elements illumination + for (Map.Entry<RealVector, Double> sd : samples.entrySet()) { + + // Vector between interpolation point and current sample point. + final RealVector diff = sd.getKey().subtract(p); + final double diffNorm = diff.getNorm(); + + if (FastMath.abs(diffNorm) < FastMath.ulp(1d)) { + // No need to interpolate, as the interpolation point is + // actually (very close to) one of the sampled points. + return sd.getValue(); + } + + for (MicrosphereSurfaceElement md : microsphere) { + final double w = FastMath.pow(diffNorm, -brightnessExponent); + md.store(cosAngle(diff, md.normal()) * w, sd); + } + + } + + // Interpolation calculation. + double value = 0; + double totalWeight = 0; + for (MicrosphereSurfaceElement md : microsphere) { + final double iV = md.illumination(); + final Map.Entry<RealVector, Double> sd = md.sample(); + if (sd != null) { + value += iV * sd.getValue(); + totalWeight += iV; + } + } + + return value / totalWeight; + + } + + /** + * Compute the cosine of the angle between 2 vectors. + * + * @param v Vector. + * @param w Vector. + * @return cosine of the angle + */ + private double cosAngle(final RealVector v, final RealVector w) { + return v.dotProduct(w) / (v.getNorm() * w.getNorm()); + } + +} |