1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.interpolation;
18
19import java.util.ArrayList;
20import java.util.HashMap;
21import java.util.List;
22import java.util.Map;
23
24import org.apache.commons.math.DimensionMismatchException;
25import org.apache.commons.math.analysis.MultivariateRealFunction;
26import org.apache.commons.math.exception.NoDataException;
27import org.apache.commons.math.linear.ArrayRealVector;
28import org.apache.commons.math.linear.RealVector;
29import org.apache.commons.math.random.UnitSphereRandomVectorGenerator;
30import org.apache.commons.math.util.FastMath;
31
32/**
33 * Interpolating function that implements the
34 * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>.
35 *
36 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
37 */
38public class MicrosphereInterpolatingFunction
39    implements MultivariateRealFunction {
40    /**
41     * Space dimension.
42     */
43    private final int dimension;
44    /**
45     * Internal accounting data for the interpolation algorithm.
46     * Each element of the list corresponds to one surface element of
47     * the microsphere.
48     */
49    private final List<MicrosphereSurfaceElement> microsphere;
50    /**
51     * Exponent used in the power law that computes the weights of the
52     * sample data.
53     */
54    private final double brightnessExponent;
55    /**
56     * Sample data.
57     */
58    private final Map<RealVector, Double> samples;
59
60    /**
61     * Class for storing the accounting data needed to perform the
62     * microsphere projection.
63     */
64    private static class MicrosphereSurfaceElement {
65
66        /** Normal vector characterizing a surface element. */
67        private final RealVector normal;
68
69        /** Illumination received from the brightest sample. */
70        private double brightestIllumination;
71
72        /** Brightest sample. */
73        private Map.Entry<RealVector, Double> brightestSample;
74
75        /**
76         * @param n Normal vector characterizing a surface element
77         * of the microsphere.
78         */
79        MicrosphereSurfaceElement(double[] n) {
80            normal = new ArrayRealVector(n);
81        }
82
83        /**
84         * Return the normal vector.
85         * @return the normal vector
86         */
87        RealVector normal() {
88            return normal;
89        }
90
91        /**
92         * Reset "illumination" and "sampleIndex".
93         */
94        void reset() {
95            brightestIllumination = 0;
96            brightestSample = null;
97        }
98
99        /**
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}
245