001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.interpolation;
018
019 import java.util.ArrayList;
020 import java.util.HashMap;
021 import java.util.List;
022 import java.util.Map;
023
024 import org.apache.commons.math3.analysis.MultivariateFunction;
025 import org.apache.commons.math3.exception.DimensionMismatchException;
026 import org.apache.commons.math3.exception.NoDataException;
027 import org.apache.commons.math3.exception.NullArgumentException;
028 import org.apache.commons.math3.linear.ArrayRealVector;
029 import org.apache.commons.math3.linear.RealVector;
030 import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
031 import org.apache.commons.math3.util.FastMath;
032
033 /**
034 * Interpolating function that implements the
035 * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>.
036 *
037 * @version $Id: MicrosphereInterpolatingFunction.java 1379904 2012-09-01 23:54:52Z erans $
038 */
039 public class MicrosphereInterpolatingFunction
040 implements MultivariateFunction {
041 /**
042 * Space dimension.
043 */
044 private final int dimension;
045 /**
046 * Internal accounting data for the interpolation algorithm.
047 * Each element of the list corresponds to one surface element of
048 * the microsphere.
049 */
050 private final List<MicrosphereSurfaceElement> microsphere;
051 /**
052 * Exponent used in the power law that computes the weights of the
053 * sample data.
054 */
055 private final double brightnessExponent;
056 /**
057 * Sample data.
058 */
059 private final Map<RealVector, Double> samples;
060
061 /**
062 * Class for storing the accounting data needed to perform the
063 * microsphere projection.
064 */
065 private static class MicrosphereSurfaceElement {
066 /** Normal vector characterizing a surface element. */
067 private final RealVector normal;
068 /** Illumination received from the brightest sample. */
069 private double brightestIllumination;
070 /** Brightest sample. */
071 private Map.Entry<RealVector, Double> brightestSample;
072
073 /**
074 * @param n Normal vector characterizing a surface element
075 * of the microsphere.
076 */
077 MicrosphereSurfaceElement(double[] n) {
078 normal = new ArrayRealVector(n);
079 }
080
081 /**
082 * Return the normal vector.
083 * @return the normal vector
084 */
085 RealVector normal() {
086 return normal;
087 }
088
089 /**
090 * Reset "illumination" and "sampleIndex".
091 */
092 void reset() {
093 brightestIllumination = 0;
094 brightestSample = null;
095 }
096
097 /**
098 * Store the illumination and index of the brightest sample.
099 * @param illuminationFromSample illumination received from sample
100 * @param sample current sample illuminating the element
101 */
102 void store(final double illuminationFromSample,
103 final Map.Entry<RealVector, Double> sample) {
104 if (illuminationFromSample > this.brightestIllumination) {
105 this.brightestIllumination = illuminationFromSample;
106 this.brightestSample = sample;
107 }
108 }
109
110 /**
111 * Get the illumination of the element.
112 * @return the illumination.
113 */
114 double illumination() {
115 return brightestIllumination;
116 }
117
118 /**
119 * Get the sample illuminating the element the most.
120 * @return the sample.
121 */
122 Map.Entry<RealVector, Double> sample() {
123 return brightestSample;
124 }
125 }
126
127 /**
128 * @param xval Arguments for the interpolation points.
129 * {@code xval[i][0]} is the first component of interpolation point
130 * {@code i}, {@code xval[i][1]} is the second component, and so on
131 * until {@code xval[i][d-1]}, the last component of that interpolation
132 * point (where {@code dimension} is thus the dimension of the sampled
133 * space).
134 * @param yval Values for the interpolation points.
135 * @param brightnessExponent Brightness dimming factor.
136 * @param microsphereElements Number of surface elements of the
137 * microsphere.
138 * @param rand Unit vector generator for creating the microsphere.
139 * @throws DimensionMismatchException if the lengths of {@code yval} and
140 * {@code xval} (equal to {@code n}, the number of interpolation points)
141 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
142 * have lengths different from {@code dimension}.
143 * @throws NoDataException if there an array has zero-length.
144 * @throws NullArgumentException if an argument is {@code null}.
145 */
146 public MicrosphereInterpolatingFunction(double[][] xval,
147 double[] yval,
148 int brightnessExponent,
149 int microsphereElements,
150 UnitSphereRandomVectorGenerator rand)
151 throws DimensionMismatchException,
152 NoDataException,
153 NullArgumentException {
154 if (xval == null ||
155 yval == null) {
156 throw new NullArgumentException();
157 }
158 if (xval.length == 0) {
159 throw new NoDataException();
160 }
161 if (xval.length != yval.length) {
162 throw new DimensionMismatchException(xval.length, yval.length);
163 }
164 if (xval[0] == null) {
165 throw new NullArgumentException();
166 }
167
168 dimension = xval[0].length;
169 this.brightnessExponent = brightnessExponent;
170
171 // Copy data samples.
172 samples = new HashMap<RealVector, Double>(yval.length);
173 for (int i = 0; i < xval.length; ++i) {
174 final double[] xvalI = xval[i];
175 if (xvalI == null) {
176 throw new NullArgumentException();
177 }
178 if (xvalI.length != dimension) {
179 throw new DimensionMismatchException(xvalI.length, dimension);
180 }
181
182 samples.put(new ArrayRealVector(xvalI), yval[i]);
183 }
184
185 microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
186 // Generate the microsphere, assuming that a fairly large number of
187 // randomly generated normals will represent a sphere.
188 for (int i = 0; i < microsphereElements; i++) {
189 microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
190 }
191 }
192
193 /**
194 * @param point Interpolation point.
195 * @return the interpolated value.
196 */
197 public double value(double[] point) {
198 final RealVector p = new ArrayRealVector(point);
199
200 // Reset.
201 for (MicrosphereSurfaceElement md : microsphere) {
202 md.reset();
203 }
204
205 // Compute contribution of each sample points to the microsphere elements illumination
206 for (Map.Entry<RealVector, Double> sd : samples.entrySet()) {
207
208 // Vector between interpolation point and current sample point.
209 final RealVector diff = sd.getKey().subtract(p);
210 final double diffNorm = diff.getNorm();
211
212 if (FastMath.abs(diffNorm) < FastMath.ulp(1d)) {
213 // No need to interpolate, as the interpolation point is
214 // actually (very close to) one of the sampled points.
215 return sd.getValue();
216 }
217
218 for (MicrosphereSurfaceElement md : microsphere) {
219 final double w = FastMath.pow(diffNorm, -brightnessExponent);
220 md.store(cosAngle(diff, md.normal()) * w, sd);
221 }
222
223 }
224
225 // Interpolation calculation.
226 double value = 0;
227 double totalWeight = 0;
228 for (MicrosphereSurfaceElement md : microsphere) {
229 final double iV = md.illumination();
230 final Map.Entry<RealVector, Double> sd = md.sample();
231 if (sd != null) {
232 value += iV * sd.getValue();
233 totalWeight += iV;
234 }
235 }
236
237 return value / totalWeight;
238 }
239
240 /**
241 * Compute the cosine of the angle between 2 vectors.
242 *
243 * @param v Vector.
244 * @param w Vector.
245 * @return the cosine of the angle between {@code v} and {@code w}.
246 */
247 private double cosAngle(final RealVector v, final RealVector w) {
248 return v.dotProduct(w) / (v.getNorm() * w.getNorm());
249 }
250 }