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    }