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