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