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.math.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.math.DimensionMismatchException;
025 import org.apache.commons.math.MathRuntimeException;
026 import org.apache.commons.math.analysis.MultivariateRealFunction;
027 import org.apache.commons.math.linear.ArrayRealVector;
028 import org.apache.commons.math.linear.RealVector;
029 import org.apache.commons.math.random.UnitSphereRandomVectorGenerator;
030
031 /**
032 * Interpolating function that implements the
033 * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>.
034 *
035 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
036 */
037 public class MicrosphereInterpolatingFunction
038 implements MultivariateRealFunction {
039 /**
040 * Space dimension.
041 */
042 private final int dimension;
043 /**
044 * Internal accounting data for the interpolation algorithm.
045 * Each element of the list corresponds to one surface element of
046 * the microsphere.
047 */
048 private final List<MicrosphereSurfaceElement> microsphere;
049 /**
050 * Exponent used in the power law that computes the weights of the
051 * sample data.
052 */
053 private final double brightnessExponent;
054 /**
055 * Sample data.
056 */
057 private final Map<RealVector, Double> samples;
058
059 /**
060 * Class for storing the accounting data needed to perform the
061 * microsphere projection.
062 */
063 private static class MicrosphereSurfaceElement {
064
065 /** Normal vector characterizing a surface element. */
066 private final RealVector normal;
067
068 /** Illumination received from the brightest sample. */
069 private double brightestIllumination;
070
071 /** Brightest sample. */
072 private Map.Entry<RealVector, Double> brightestSample;
073
074 /**
075 * @param n Normal vector characterizing a surface element
076 * of the microsphere.
077 */
078 MicrosphereSurfaceElement(double[] n) {
079 normal = new ArrayRealVector(n);
080 }
081
082 /**
083 * Return the normal vector.
084 * @return the normal vector
085 */
086 RealVector normal() {
087 return normal;
088 }
089
090 /**
091 * Reset "illumination" and "sampleIndex".
092 */
093 void reset() {
094 brightestIllumination = 0;
095 brightestSample = null;
096 }
097
098 /**
099 * Store the illumination and index of the brightest sample.
100 * @param illuminationFromSample illumination received from sample
101 * @param sample current sample illuminating the element
102 */
103 void store(final double illuminationFromSample,
104 final Map.Entry<RealVector, Double> sample) {
105 if (illuminationFromSample > this.brightestIllumination) {
106 this.brightestIllumination = illuminationFromSample;
107 this.brightestSample = sample;
108 }
109 }
110
111 /**
112 * Get the illumination of the element.
113 * @return the illumination.
114 */
115 double illumination() {
116 return brightestIllumination;
117 }
118
119 /**
120 * Get the sample illuminating the element the most.
121 * @return the sample.
122 */
123 Map.Entry<RealVector, Double> sample() {
124 return brightestSample;
125 }
126 }
127
128 /**
129 * @param xval the arguments for the interpolation points.
130 * {@code xval[i][0]} is the first component of interpolation point
131 * {@code i}, {@code xval[i][1]} is the second component, and so on
132 * until {@code xval[i][d-1]}, the last component of that interpolation
133 * point (where {@code dimension} is thus the dimension of the sampled
134 * space).
135 * @param yval the values for the interpolation points
136 * @param brightnessExponent Brightness dimming factor.
137 * @param microsphereElements Number of surface elements of the
138 * microsphere.
139 * @param rand Unit vector generator for creating the microsphere.
140 * @throws DimensionMismatchException if the lengths of {@code yval} and
141 * {@code xval} (equal to {@code n}, the number of interpolation points)
142 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
143 * have lengths different from {@code dimension}.
144 * @throws IllegalArgumentException if there are no data (xval null or zero length)
145 */
146 public MicrosphereInterpolatingFunction(double[][] xval,
147 double[] yval,
148 int brightnessExponent,
149 int microsphereElements,
150 UnitSphereRandomVectorGenerator rand)
151 throws DimensionMismatchException, IllegalArgumentException {
152 if (xval.length == 0 || xval[0] == null) {
153 throw MathRuntimeException.createIllegalArgumentException("no data");
154 }
155
156 if (xval.length != yval.length) {
157 throw new DimensionMismatchException(xval.length, yval.length);
158 }
159
160 dimension = xval[0].length;
161 this.brightnessExponent = brightnessExponent;
162
163 // Copy data samples.
164 samples = new HashMap<RealVector, Double>(yval.length);
165 for (int i = 0; i < xval.length; ++i) {
166 final double[] xvalI = xval[i];
167 if ( xvalI.length != dimension) {
168 throw new DimensionMismatchException(xvalI.length, dimension);
169 }
170
171 samples.put(new ArrayRealVector(xvalI), yval[i]);
172 }
173
174 microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
175 // Generate the microsphere, assuming that a fairly large number of
176 // randomly generated normals will represent a sphere.
177 for (int i = 0; i < microsphereElements; i++) {
178 microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
179 }
180
181 }
182
183 /**
184 * @param point Interpolation point.
185 * @return the interpolated value.
186 */
187 public double value(double[] point) {
188
189 final RealVector p = new ArrayRealVector(point);
190
191 // Reset.
192 for (MicrosphereSurfaceElement md : microsphere) {
193 md.reset();
194 }
195
196 // Compute contribution of each sample points to the microsphere elements illumination
197 for (Map.Entry<RealVector, Double> sd : samples.entrySet()) {
198
199 // Vector between interpolation point and current sample point.
200 final RealVector diff = sd.getKey().subtract(p);
201 final double diffNorm = diff.getNorm();
202
203 if (Math.abs(diffNorm) < Math.ulp(1d)) {
204 // No need to interpolate, as the interpolation point is
205 // actually (very close to) one of the sampled points.
206 return sd.getValue();
207 }
208
209 for (MicrosphereSurfaceElement md : microsphere) {
210 final double w = Math.pow(diffNorm, -brightnessExponent);
211 md.store(cosAngle(diff, md.normal()) * w, sd);
212 }
213
214 }
215
216 // Interpolation calculation.
217 double value = 0;
218 double totalWeight = 0;
219 for (MicrosphereSurfaceElement md : microsphere) {
220 final double iV = md.illumination();
221 final Map.Entry<RealVector, Double> sd = md.sample();
222 if (sd != null) {
223 value += iV * sd.getValue();
224 totalWeight += iV;
225 }
226 }
227
228 return value / totalWeight;
229
230 }
231
232 /**
233 * Compute the cosine of the angle between 2 vectors.
234 *
235 * @param v Vector.
236 * @param w Vector.
237 * @return cosine of the angle
238 */
239 private double cosAngle(final RealVector v, final RealVector w) {
240 return v.dotProduct(w) / (v.getNorm() * w.getNorm());
241 }
242
243 }