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