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 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import java.util.List;
020import java.util.ArrayList;
021import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
022import org.apache.commons.math3.exception.DimensionMismatchException;
023import org.apache.commons.math3.exception.NotPositiveException;
024import org.apache.commons.math3.exception.NotStrictlyPositiveException;
025import org.apache.commons.math3.exception.MaxCountExceededException;
026import org.apache.commons.math3.exception.OutOfRangeException;
027import org.apache.commons.math3.util.FastMath;
028import org.apache.commons.math3.util.MathArrays;
029
030/**
031 * Utility class for the {@link MicrosphereProjectionInterpolator} algorithm.
032 *
033 * @since 3.6
034 */
035public class InterpolatingMicrosphere {
036    /** Microsphere. */
037    private final List<Facet> microsphere;
038    /** Microsphere data. */
039    private final List<FacetData> microsphereData;
040    /** Space dimension. */
041    private final int dimension;
042    /** Number of surface elements. */
043    private final int size;
044    /** Maximum fraction of the facets that can be dark. */
045    private final double maxDarkFraction;
046    /** Lowest non-zero illumination. */
047    private final double darkThreshold;
048    /** Background value. */
049    private final double background;
050
051    /**
052     * Create an unitialiazed sphere.
053     * Sub-classes are responsible for calling the {@code add(double[]) add}
054     * method in order to initialize all the sphere's facets.
055     *
056     * @param dimension Dimension of the data space.
057     * @param size Number of surface elements of the sphere.
058     * @param maxDarkFraction Maximum fraction of the facets that can be dark.
059     * If the fraction of "non-illuminated" facets is larger, no estimation
060     * of the value will be performed, and the {@code background} value will
061     * be returned instead.
062     * @param darkThreshold Value of the illumination below which a facet is
063     * considered dark.
064     * @param background Value returned when the {@code maxDarkFraction}
065     * threshold is exceeded.
066     * @throws NotStrictlyPositiveException if {@code dimension <= 0}
067     * or {@code size <= 0}.
068     * @throws NotPositiveException if {@code darkThreshold < 0}.
069     * @throws OutOfRangeException if {@code maxDarkFraction} does not
070     * belong to the interval {@code [0, 1]}.
071     */
072    protected InterpolatingMicrosphere(int dimension,
073                                       int size,
074                                       double maxDarkFraction,
075                                       double darkThreshold,
076                                       double background) {
077        if (dimension <= 0) {
078            throw new NotStrictlyPositiveException(dimension);
079        }
080        if (size <= 0) {
081            throw new NotStrictlyPositiveException(size);
082        }
083        if (maxDarkFraction < 0 ||
084            maxDarkFraction > 1) {
085            throw new OutOfRangeException(maxDarkFraction, 0, 1);
086        }
087        if (darkThreshold < 0) {
088            throw new NotPositiveException(darkThreshold);
089        }
090
091        this.dimension = dimension;
092        this.size = size;
093        this.maxDarkFraction = maxDarkFraction;
094        this.darkThreshold = darkThreshold;
095        this.background = background;
096        microsphere = new ArrayList<Facet>(size);
097        microsphereData = new ArrayList<FacetData>(size);
098    }
099
100    /**
101     * Create a sphere from randomly sampled vectors.
102     *
103     * @param dimension Dimension of the data space.
104     * @param size Number of surface elements of the sphere.
105     * @param rand Unit vector generator for creating the microsphere.
106     * @param maxDarkFraction Maximum fraction of the facets that can be dark.
107     * If the fraction of "non-illuminated" facets is larger, no estimation
108     * of the value will be performed, and the {@code background} value will
109     * be returned instead.
110     * @param darkThreshold Value of the illumination below which a facet
111     * is considered dark.
112     * @param background Value returned when the {@code maxDarkFraction}
113     * threshold is exceeded.
114     * @throws DimensionMismatchException if the size of the generated
115     * vectors does not match the dimension set in the constructor.
116     * @throws NotStrictlyPositiveException if {@code dimension <= 0}
117     * or {@code size <= 0}.
118     * @throws NotPositiveException if {@code darkThreshold < 0}.
119     * @throws OutOfRangeException if {@code maxDarkFraction} does not
120     * belong to the interval {@code [0, 1]}.
121     */
122    public InterpolatingMicrosphere(int dimension,
123                                    int size,
124                                    double maxDarkFraction,
125                                    double darkThreshold,
126                                    double background,
127                                    UnitSphereRandomVectorGenerator rand) {
128        this(dimension, size, maxDarkFraction, darkThreshold, background);
129
130        // Generate the microsphere normals, assuming that a number of
131        // randomly generated normals will represent a sphere.
132        for (int i = 0; i < size; i++) {
133            add(rand.nextVector(), false);
134        }
135    }
136
137    /**
138     * Copy constructor.
139     *
140     * @param other Instance to copy.
141     */
142    protected InterpolatingMicrosphere(InterpolatingMicrosphere other) {
143        dimension = other.dimension;
144        size = other.size;
145        maxDarkFraction = other.maxDarkFraction;
146        darkThreshold = other.darkThreshold;
147        background = other.background;
148
149        // Field can be shared.
150        microsphere = other.microsphere;
151
152        // Field must be copied.
153        microsphereData = new ArrayList<FacetData>(size);
154        for (FacetData fd : other.microsphereData) {
155            microsphereData.add(new FacetData(fd.illumination(), fd.sample()));
156        }
157    }
158
159    /**
160     * Perform a copy.
161     *
162     * @return a copy of this instance.
163     */
164    public InterpolatingMicrosphere copy() {
165        return new InterpolatingMicrosphere(this);
166    }
167
168    /**
169     * Get the space dimensionality.
170     *
171     * @return the number of space dimensions.
172     */
173    public int getDimension() {
174        return dimension;
175    }
176
177    /**
178     * Get the size of the sphere.
179     *
180     * @return the number of surface elements of the microspshere.
181     */
182    public int getSize() {
183        return size;
184    }
185
186    /**
187     * Estimate the value at the requested location.
188     * This microsphere is placed at the given {@code point}, contribution
189     * of the given {@code samplePoints} to each sphere facet is computed
190     * (illumination) and the interpolation is performed (integration of
191     * the illumination).
192     *
193     * @param point Interpolation point.
194     * @param samplePoints Sampling data points.
195     * @param sampleValues Sampling data values at the corresponding
196     * {@code samplePoints}.
197     * @param exponent Exponent used in the power law that computes
198     * the weights (distance dimming factor) of the sample data.
199     * @param noInterpolationTolerance When the distance between the
200     * {@code point} and one of the {@code samplePoints} is less than
201     * this value, no interpolation will be performed, and the value
202     * of the sample will just be returned.
203     * @return the estimated value at the given {@code point}.
204     * @throws NotPositiveException if {@code exponent < 0}.
205     */
206    public double value(double[] point,
207                        double[][] samplePoints,
208                        double[] sampleValues,
209                        double exponent,
210                        double noInterpolationTolerance) {
211        if (exponent < 0) {
212            throw new NotPositiveException(exponent);
213        }
214
215        clear();
216
217        // Contribution of each sample point to the illumination of the
218        // microsphere's facets.
219        final int numSamples = samplePoints.length;
220        for (int i = 0; i < numSamples; i++) {
221            // Vector between interpolation point and current sample point.
222            final double[] diff = MathArrays.ebeSubtract(samplePoints[i], point);
223            final double diffNorm = MathArrays.safeNorm(diff);
224
225            if (FastMath.abs(diffNorm) < noInterpolationTolerance) {
226                // No need to interpolate, as the interpolation point is
227                // actually (very close to) one of the sampled points.
228                return sampleValues[i];
229            }
230
231            final double weight = FastMath.pow(diffNorm, -exponent);
232            illuminate(diff, sampleValues[i], weight);
233        }
234
235        return interpolate();
236    }
237
238    /**
239     * Replace {@code i}-th facet of the microsphere.
240     * Method for initializing the microsphere facets.
241     *
242     * @param normal Facet's normal vector.
243     * @param copy Whether to copy the given array.
244     * @throws DimensionMismatchException if the length of {@code n}
245     * does not match the space dimension.
246     * @throws MaxCountExceededException if the method has been called
247     * more times than the size of the sphere.
248     */
249    protected void add(double[] normal,
250                       boolean copy) {
251        if (microsphere.size() >= size) {
252            throw new MaxCountExceededException(size);
253        }
254        if (normal.length > dimension) {
255            throw new DimensionMismatchException(normal.length, dimension);
256        }
257
258        microsphere.add(new Facet(copy ? normal.clone() : normal));
259        microsphereData.add(new FacetData(0d, 0d));
260    }
261
262    /**
263     * Interpolation.
264     *
265     * @return the value estimated from the current illumination of the
266     * microsphere.
267     */
268    private double interpolate() {
269        // Number of non-illuminated facets.
270        int darkCount = 0;
271
272        double value = 0;
273        double totalWeight = 0;
274        for (FacetData fd : microsphereData) {
275            final double iV = fd.illumination();
276            if (iV != 0d) {
277                value += iV * fd.sample();
278                totalWeight += iV;
279            } else {
280                ++darkCount;
281            }
282        }
283
284        final double darkFraction = darkCount / (double) size;
285
286        return darkFraction <= maxDarkFraction ?
287            value / totalWeight :
288            background;
289    }
290
291    /**
292     * Illumination.
293     *
294     * @param sampleDirection Vector whose origin is at the interpolation
295     * point and tail is at the sample location.
296     * @param sampleValue Data value of the sample.
297     * @param weight Weight.
298     */
299    private void illuminate(double[] sampleDirection,
300                            double sampleValue,
301                            double weight) {
302        for (int i = 0; i < size; i++) {
303            final double[] n = microsphere.get(i).getNormal();
304            final double cos = MathArrays.cosAngle(n, sampleDirection);
305
306            if (cos > 0) {
307                final double illumination = cos * weight;
308
309                if (illumination > darkThreshold &&
310                    illumination > microsphereData.get(i).illumination()) {
311                    microsphereData.set(i, new FacetData(illumination, sampleValue));
312                }
313            }
314        }
315    }
316
317    /**
318     * Reset the all the {@link Facet facets} data to zero.
319     */
320    private void clear() {
321        for (int i = 0; i < size; i++) {
322            microsphereData.set(i, new FacetData(0d, 0d));
323        }
324    }
325
326    /**
327     * Microsphere "facet" (surface element).
328     */
329    private static class Facet {
330        /** Normal vector characterizing a surface element. */
331        private final double[] normal;
332
333        /**
334         * @param n Normal vector characterizing a surface element
335         * of the microsphere. No copy is made.
336         */
337        Facet(double[] n) {
338            normal = n;
339        }
340
341        /**
342         * Return a reference to the vector normal to this facet.
343         *
344         * @return the normal vector.
345         */
346        public double[] getNormal() {
347            return normal;
348        }
349    }
350
351    /**
352     * Data associated with each {@link Facet}.
353     */
354    private static class FacetData {
355        /** Illumination received from the sample. */
356        private final double illumination;
357        /** Data value of the sample. */
358        private final double sample;
359
360        /**
361         * @param illumination Illumination.
362         * @param sample Data value.
363         */
364        FacetData(double illumination, double sample) {
365            this.illumination = illumination;
366            this.sample = sample;
367        }
368
369        /**
370         * Get the illumination.
371         * @return the illumination.
372         */
373        public double illumination() {
374            return illumination;
375        }
376
377        /**
378         * Get the data value.
379         * @return the data value.
380         */
381        public double sample() {
382            return sample;
383        }
384    }
385}