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