Home | History | Annotate | Download | only in interpolation 
      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  */
     17 package org.apache.commons.math.analysis.interpolation;
     18 
     19 import java.util.ArrayList;
     20 import java.util.HashMap;
     21 import java.util.List;
     22 import java.util.Map;
     23 
     24 import org.apache.commons.math.DimensionMismatchException;
     25 import org.apache.commons.math.analysis.MultivariateRealFunction;
     26 import org.apache.commons.math.exception.NoDataException;
     27 import org.apache.commons.math.linear.ArrayRealVector;
     28 import org.apache.commons.math.linear.RealVector;
     29 import org.apache.commons.math.random.UnitSphereRandomVectorGenerator;
     30 import 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 aot 2010) $
     37  */
     38 public 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