    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.math3.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.math3.analysis.MultivariateFunction;
    025    import org.apache.commons.math3.exception.DimensionMismatchException;
    026    import org.apache.commons.math3.exception.NoDataException;
    027    import org.apache.commons.math3.exception.NullArgumentException;
    028    import org.apache.commons.math3.linear.ArrayRealVector;
    029    import org.apache.commons.math3.linear.RealVector;
    030    import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
    031    import org.apache.commons.math3.util.FastMath;
    032
    033    /**
    034     * Interpolating function that implements the
    035     * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>.
    036     *
    037     * @version $Id: MicrosphereInterpolatingFunction.java 1379904 2012-09-01 23:54:52Z erans $
    038     */
    039    public class MicrosphereInterpolatingFunction
    040        implements MultivariateFunction {
    041        /**
    042         * Space dimension.
    043         */
    044        private final int dimension;
    045        /**
    046         * Internal accounting data for the interpolation algorithm.
    047         * Each element of the list corresponds to one surface element of
    048         * the microsphere.
    049         */
    050        private final List<MicrosphereSurfaceElement> microsphere;
    051        /**
    052         * Exponent used in the power law that computes the weights of the
    053         * sample data.
    054         */
    055        private final double brightnessExponent;
    056        /**
    057         * Sample data.
    058         */
    059        private final Map<RealVector, Double> samples;
    060
    061        /**
    062         * Class for storing the accounting data needed to perform the
    063         * microsphere projection.
    064         */
    065        private static class MicrosphereSurfaceElement {
    066            /** Normal vector characterizing a surface element. */
    067            private final RealVector normal;
    068            /** Illumination received from the brightest sample. */
    069            private double brightestIllumination;
    070            /** Brightest sample. */
    071            private Map.Entry<RealVector, Double> brightestSample;
    072
    073            /**
    074             * @param n Normal vector characterizing a surface element
    075             * of the microsphere.
    076             */
    077            MicrosphereSurfaceElement(double[] n) {
    078                normal = new ArrayRealVector(n);
    079            }
    080
    081            /**
    082             * Return the normal vector.
    083             * @return the normal vector
    084             */
    085            RealVector normal() {
    086                return normal;
    087            }
    088
    089            /**
    090             * Reset "illumination" and "sampleIndex".
    091             */
    092            void reset() {
    093                brightestIllumination = 0;
    094                brightestSample = null;
    095            }
    096
    097            /**
    098             * Store the illumination and index of the brightest sample.
    099             * @param illuminationFromSample illumination received from sample
    100             * @param sample current sample illuminating the element
    101             */
    102            void store(final double illuminationFromSample,
    103                       final Map.Entry<RealVector, Double> sample) {
    104                if (illuminationFromSample > this.brightestIllumination) {
    105                    this.brightestIllumination = illuminationFromSample;
    106                    this.brightestSample = sample;
    107                }
    108            }
    109
    110            /**
    111             * Get the illumination of the element.
    112             * @return the illumination.
    113             */
    114            double illumination() {
    115                return brightestIllumination;
    116            }
    117
    118            /**
    119             * Get the sample illuminating the element the most.
    120             * @return the sample.
    121             */
    122            Map.Entry<RealVector, Double> sample() {
    123                return brightestSample;
    124            }
    125        }
    126
    127        /**
    128         * @param xval Arguments for the interpolation points.
    129         * {@code xval[i][0]} is the first component of interpolation point
    130         * {@code i}, {@code xval[i][1]} is the second component, and so on
    131         * until {@code xval[i][d-1]}, the last component of that interpolation
    132         * point (where {@code dimension} is thus the dimension of the sampled
    133         * space).
    134         * @param yval Values for the interpolation points.
    135         * @param brightnessExponent Brightness dimming factor.
    136         * @param microsphereElements Number of surface elements of the
    137         * microsphere.
    138         * @param rand Unit vector generator for creating the microsphere.
    139         * @throws DimensionMismatchException if the lengths of {@code yval} and
    140         * {@code xval} (equal to {@code n}, the number of interpolation points)
    141         * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
    142         * have lengths different from {@code dimension}.
    143         * @throws NoDataException if there an array has zero-length.
    144         * @throws NullArgumentException if an argument is {@code null}.
    145         */
    146        public MicrosphereInterpolatingFunction(double[][] xval,
    147                                                double[] yval,
    148                                                int brightnessExponent,
    149                                                int microsphereElements,
    150                                                UnitSphereRandomVectorGenerator rand)
    151            throws DimensionMismatchException,
    152                   NoDataException,
    153                   NullArgumentException {
    154            if (xval == null ||
    155                yval == null) {
    156                throw new NullArgumentException();
    157            }
    158            if (xval.length == 0) {
    159                throw new NoDataException();
    160            }
    161            if (xval.length != yval.length) {
    162                throw new DimensionMismatchException(xval.length, yval.length);
    163            }
    164            if (xval[0] == null) {
    165                throw new NullArgumentException();
    166            }
    167
    168            dimension = xval[0].length;
    169            this.brightnessExponent = brightnessExponent;
    170
    171            // Copy data samples.
    172            samples = new HashMap<RealVector, Double>(yval.length);
    173            for (int i = 0; i < xval.length; ++i) {
    174                final double[] xvalI = xval[i];
    175                if (xvalI == null) {
    176                    throw new NullArgumentException();
    177                }
    178                if (xvalI.length != dimension) {
    179                    throw new DimensionMismatchException(xvalI.length, dimension);
    180                }
    181
    182                samples.put(new ArrayRealVector(xvalI), yval[i]);
    183            }
    184
    185            microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
    186            // Generate the microsphere, assuming that a fairly large number of
    187            // randomly generated normals will represent a sphere.
    188            for (int i = 0; i < microsphereElements; i++) {
    189                microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
    190            }
    191        }
    192
    193        /**
    194         * @param point Interpolation point.
    195         * @return the interpolated value.
    196         */
    197        public double value(double[] point) {
    198            final RealVector p = new ArrayRealVector(point);
    199
    200            // Reset.
    201            for (MicrosphereSurfaceElement md : microsphere) {
    202                md.reset();
    203            }
    204
    205            // Compute contribution of each sample points to the microsphere elements illumination
    206            for (Map.Entry<RealVector, Double> sd : samples.entrySet()) {
    207
    208                // Vector between interpolation point and current sample point.
    209                final RealVector diff = sd.getKey().subtract(p);
    210                final double diffNorm = diff.getNorm();
    211
    212                if (FastMath.abs(diffNorm) < FastMath.ulp(1d)) {
    213                    // No need to interpolate, as the interpolation point is
    214                    // actually (very close to) one of the sampled points.
    215                    return sd.getValue();
    216                }
    217
    218                for (MicrosphereSurfaceElement md : microsphere) {
    219                    final double w = FastMath.pow(diffNorm, -brightnessExponent);
    220                    md.store(cosAngle(diff, md.normal()) * w, sd);
    221                }
    222
    223            }
    224
    225            // Interpolation calculation.
    226            double value = 0;
    227            double totalWeight = 0;
    228            for (MicrosphereSurfaceElement md : microsphere) {
    229                final double iV = md.illumination();
    230                final Map.Entry<RealVector, Double> sd = md.sample();
    231                if (sd != null) {
    232                    value += iV * sd.getValue();
    233                    totalWeight += iV;
    234                }
    235            }
    236
    237            return value / totalWeight;
    238        }
    239
    240        /**
    241         * Compute the cosine of the angle between 2 vectors.
    242         *
    243         * @param v Vector.
    244         * @param w Vector.
    245         * @return the cosine of the angle between {@code v} and {@code w}.
    246         */
    247        private double cosAngle(final RealVector v, final RealVector w) {
    248            return v.dotProduct(w) / (v.getNorm() * w.getNorm());
    249        }
    250    }