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 */ 017package org.apache.commons.math3.analysis.interpolation; 018 019import java.util.ArrayList; 020import java.util.HashMap; 021import java.util.List; 022import java.util.Map; 023 024import org.apache.commons.math3.analysis.MultivariateFunction; 025import org.apache.commons.math3.exception.DimensionMismatchException; 026import org.apache.commons.math3.exception.NoDataException; 027import org.apache.commons.math3.exception.NullArgumentException; 028import org.apache.commons.math3.linear.ArrayRealVector; 029import org.apache.commons.math3.linear.RealVector; 030import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; 031import 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 */ 038public class MicrosphereInterpolatingFunction 039 implements MultivariateFunction { 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 /** Normal vector characterizing a surface element. */ 066 private final RealVector normal; 067 /** Illumination received from the brightest sample. */ 068 private double brightestIllumination; 069 /** Brightest sample. */ 070 private Map.Entry<RealVector, Double> brightestSample; 071 072 /** 073 * @param n Normal vector characterizing a surface element 074 * of the microsphere. 075 */ 076 MicrosphereSurfaceElement(double[] n) { 077 normal = new ArrayRealVector(n); 078 } 079 080 /** 081 * Return the normal vector. 082 * @return the normal vector 083 */ 084 RealVector normal() { 085 return normal; 086 } 087 088 /** 089 * Reset "illumination" and "sampleIndex". 090 */ 091 void reset() { 092 brightestIllumination = 0; 093 brightestSample = null; 094 } 095 096 /** 097 * Store the illumination and index of the brightest sample. 098 * @param illuminationFromSample illumination received from sample 099 * @param sample current sample illuminating the element 100 */ 101 void store(final double illuminationFromSample, 102 final Map.Entry<RealVector, Double> sample) { 103 if (illuminationFromSample > this.brightestIllumination) { 104 this.brightestIllumination = illuminationFromSample; 105 this.brightestSample = sample; 106 } 107 } 108 109 /** 110 * Get the illumination of the element. 111 * @return the illumination. 112 */ 113 double illumination() { 114 return brightestIllumination; 115 } 116 117 /** 118 * Get the sample illuminating the element the most. 119 * @return the sample. 120 */ 121 Map.Entry<RealVector, Double> sample() { 122 return brightestSample; 123 } 124 } 125 126 /** 127 * @param xval Arguments for the interpolation points. 128 * {@code xval[i][0]} is the first component of interpolation point 129 * {@code i}, {@code xval[i][1]} is the second component, and so on 130 * until {@code xval[i][d-1]}, the last component of that interpolation 131 * point (where {@code dimension} is thus the dimension of the sampled 132 * space). 133 * @param yval Values for the interpolation points. 134 * @param brightnessExponent Brightness dimming factor. 135 * @param microsphereElements Number of surface elements of the 136 * microsphere. 137 * @param rand Unit vector generator for creating the microsphere. 138 * @throws DimensionMismatchException if the lengths of {@code yval} and 139 * {@code xval} (equal to {@code n}, the number of interpolation points) 140 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]}, 141 * have lengths different from {@code dimension}. 142 * @throws NoDataException if there an array has zero-length. 143 * @throws NullArgumentException if an argument is {@code null}. 144 */ 145 public MicrosphereInterpolatingFunction(double[][] xval, 146 double[] yval, 147 int brightnessExponent, 148 int microsphereElements, 149 UnitSphereRandomVectorGenerator rand) 150 throws DimensionMismatchException, 151 NoDataException, 152 NullArgumentException { 153 if (xval == null || 154 yval == null) { 155 throw new NullArgumentException(); 156 } 157 if (xval.length == 0) { 158 throw new NoDataException(); 159 } 160 if (xval.length != yval.length) { 161 throw new DimensionMismatchException(xval.length, yval.length); 162 } 163 if (xval[0] == null) { 164 throw new NullArgumentException(); 165 } 166 167 dimension = xval[0].length; 168 this.brightnessExponent = brightnessExponent; 169 170 // Copy data samples. 171 samples = new HashMap<RealVector, Double>(yval.length); 172 for (int i = 0; i < xval.length; ++i) { 173 final double[] xvalI = xval[i]; 174 if (xvalI == null) { 175 throw new NullArgumentException(); 176 } 177 if (xvalI.length != dimension) { 178 throw new DimensionMismatchException(xvalI.length, dimension); 179 } 180 181 samples.put(new ArrayRealVector(xvalI), yval[i]); 182 } 183 184 microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements); 185 // Generate the microsphere, assuming that a fairly large number of 186 // randomly generated normals will represent a sphere. 187 for (int i = 0; i < microsphereElements; i++) { 188 microsphere.add(new MicrosphereSurfaceElement(rand.nextVector())); 189 } 190 } 191 192 /** 193 * @param point Interpolation point. 194 * @return the interpolated value. 195 * @throws DimensionMismatchException if point dimension does not math sample 196 */ 197 public double value(double[] point) throws DimensionMismatchException { 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}