    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    }