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.List; 020import java.util.ArrayList; 021import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; 022import org.apache.commons.math3.exception.DimensionMismatchException; 023import org.apache.commons.math3.exception.NotPositiveException; 024import org.apache.commons.math3.exception.NotStrictlyPositiveException; 025import org.apache.commons.math3.exception.MaxCountExceededException; 026import org.apache.commons.math3.exception.OutOfRangeException; 027import org.apache.commons.math3.util.FastMath; 028import org.apache.commons.math3.util.MathArrays; 029 030/** 031 * Utility class for the {@link MicrosphereProjectionInterpolator} algorithm. 032 * 033 * @since 3.6 034 */ 035public class InterpolatingMicrosphere { 036 /** Microsphere. */ 037 private final List<Facet> microsphere; 038 /** Microsphere data. */ 039 private final List<FacetData> microsphereData; 040 /** Space dimension. */ 041 private final int dimension; 042 /** Number of surface elements. */ 043 private final int size; 044 /** Maximum fraction of the facets that can be dark. */ 045 private final double maxDarkFraction; 046 /** Lowest non-zero illumination. */ 047 private final double darkThreshold; 048 /** Background value. */ 049 private final double background; 050 051 /** 052 * Create an unitialiazed sphere. 053 * Sub-classes are responsible for calling the {@code add(double[]) add} 054 * method in order to initialize all the sphere's facets. 055 * 056 * @param dimension Dimension of the data space. 057 * @param size Number of surface elements of the sphere. 058 * @param maxDarkFraction Maximum fraction of the facets that can be dark. 059 * If the fraction of "non-illuminated" facets is larger, no estimation 060 * of the value will be performed, and the {@code background} value will 061 * be returned instead. 062 * @param darkThreshold Value of the illumination below which a facet is 063 * considered dark. 064 * @param background Value returned when the {@code maxDarkFraction} 065 * threshold is exceeded. 066 * @throws NotStrictlyPositiveException if {@code dimension <= 0} 067 * or {@code size <= 0}. 068 * @throws NotPositiveException if {@code darkThreshold < 0}. 069 * @throws OutOfRangeException if {@code maxDarkFraction} does not 070 * belong to the interval {@code [0, 1]}. 071 */ 072 protected InterpolatingMicrosphere(int dimension, 073 int size, 074 double maxDarkFraction, 075 double darkThreshold, 076 double background) { 077 if (dimension <= 0) { 078 throw new NotStrictlyPositiveException(dimension); 079 } 080 if (size <= 0) { 081 throw new NotStrictlyPositiveException(size); 082 } 083 if (maxDarkFraction < 0 || 084 maxDarkFraction > 1) { 085 throw new OutOfRangeException(maxDarkFraction, 0, 1); 086 } 087 if (darkThreshold < 0) { 088 throw new NotPositiveException(darkThreshold); 089 } 090 091 this.dimension = dimension; 092 this.size = size; 093 this.maxDarkFraction = maxDarkFraction; 094 this.darkThreshold = darkThreshold; 095 this.background = background; 096 microsphere = new ArrayList<Facet>(size); 097 microsphereData = new ArrayList<FacetData>(size); 098 } 099 100 /** 101 * Create a sphere from randomly sampled vectors. 102 * 103 * @param dimension Dimension of the data space. 104 * @param size Number of surface elements of the sphere. 105 * @param rand Unit vector generator for creating the microsphere. 106 * @param maxDarkFraction Maximum fraction of the facets that can be dark. 107 * If the fraction of "non-illuminated" facets is larger, no estimation 108 * of the value will be performed, and the {@code background} value will 109 * be returned instead. 110 * @param darkThreshold Value of the illumination below which a facet 111 * is considered dark. 112 * @param background Value returned when the {@code maxDarkFraction} 113 * threshold is exceeded. 114 * @throws DimensionMismatchException if the size of the generated 115 * vectors does not match the dimension set in the constructor. 116 * @throws NotStrictlyPositiveException if {@code dimension <= 0} 117 * or {@code size <= 0}. 118 * @throws NotPositiveException if {@code darkThreshold < 0}. 119 * @throws OutOfRangeException if {@code maxDarkFraction} does not 120 * belong to the interval {@code [0, 1]}. 121 */ 122 public InterpolatingMicrosphere(int dimension, 123 int size, 124 double maxDarkFraction, 125 double darkThreshold, 126 double background, 127 UnitSphereRandomVectorGenerator rand) { 128 this(dimension, size, maxDarkFraction, darkThreshold, background); 129 130 // Generate the microsphere normals, assuming that a number of 131 // randomly generated normals will represent a sphere. 132 for (int i = 0; i < size; i++) { 133 add(rand.nextVector(), false); 134 } 135 } 136 137 /** 138 * Copy constructor. 139 * 140 * @param other Instance to copy. 141 */ 142 protected InterpolatingMicrosphere(InterpolatingMicrosphere other) { 143 dimension = other.dimension; 144 size = other.size; 145 maxDarkFraction = other.maxDarkFraction; 146 darkThreshold = other.darkThreshold; 147 background = other.background; 148 149 // Field can be shared. 150 microsphere = other.microsphere; 151 152 // Field must be copied. 153 microsphereData = new ArrayList<FacetData>(size); 154 for (FacetData fd : other.microsphereData) { 155 microsphereData.add(new FacetData(fd.illumination(), fd.sample())); 156 } 157 } 158 159 /** 160 * Perform a copy. 161 * 162 * @return a copy of this instance. 163 */ 164 public InterpolatingMicrosphere copy() { 165 return new InterpolatingMicrosphere(this); 166 } 167 168 /** 169 * Get the space dimensionality. 170 * 171 * @return the number of space dimensions. 172 */ 173 public int getDimension() { 174 return dimension; 175 } 176 177 /** 178 * Get the size of the sphere. 179 * 180 * @return the number of surface elements of the microspshere. 181 */ 182 public int getSize() { 183 return size; 184 } 185 186 /** 187 * Estimate the value at the requested location. 188 * This microsphere is placed at the given {@code point}, contribution 189 * of the given {@code samplePoints} to each sphere facet is computed 190 * (illumination) and the interpolation is performed (integration of 191 * the illumination). 192 * 193 * @param point Interpolation point. 194 * @param samplePoints Sampling data points. 195 * @param sampleValues Sampling data values at the corresponding 196 * {@code samplePoints}. 197 * @param exponent Exponent used in the power law that computes 198 * the weights (distance dimming factor) of the sample data. 199 * @param noInterpolationTolerance When the distance between the 200 * {@code point} and one of the {@code samplePoints} is less than 201 * this value, no interpolation will be performed, and the value 202 * of the sample will just be returned. 203 * @return the estimated value at the given {@code point}. 204 * @throws NotPositiveException if {@code exponent < 0}. 205 */ 206 public double value(double[] point, 207 double[][] samplePoints, 208 double[] sampleValues, 209 double exponent, 210 double noInterpolationTolerance) { 211 if (exponent < 0) { 212 throw new NotPositiveException(exponent); 213 } 214 215 clear(); 216 217 // Contribution of each sample point to the illumination of the 218 // microsphere's facets. 219 final int numSamples = samplePoints.length; 220 for (int i = 0; i < numSamples; i++) { 221 // Vector between interpolation point and current sample point. 222 final double[] diff = MathArrays.ebeSubtract(samplePoints[i], point); 223 final double diffNorm = MathArrays.safeNorm(diff); 224 225 if (FastMath.abs(diffNorm) < noInterpolationTolerance) { 226 // No need to interpolate, as the interpolation point is 227 // actually (very close to) one of the sampled points. 228 return sampleValues[i]; 229 } 230 231 final double weight = FastMath.pow(diffNorm, -exponent); 232 illuminate(diff, sampleValues[i], weight); 233 } 234 235 return interpolate(); 236 } 237 238 /** 239 * Replace {@code i}-th facet of the microsphere. 240 * Method for initializing the microsphere facets. 241 * 242 * @param normal Facet's normal vector. 243 * @param copy Whether to copy the given array. 244 * @throws DimensionMismatchException if the length of {@code n} 245 * does not match the space dimension. 246 * @throws MaxCountExceededException if the method has been called 247 * more times than the size of the sphere. 248 */ 249 protected void add(double[] normal, 250 boolean copy) { 251 if (microsphere.size() >= size) { 252 throw new MaxCountExceededException(size); 253 } 254 if (normal.length > dimension) { 255 throw new DimensionMismatchException(normal.length, dimension); 256 } 257 258 microsphere.add(new Facet(copy ? normal.clone() : normal)); 259 microsphereData.add(new FacetData(0d, 0d)); 260 } 261 262 /** 263 * Interpolation. 264 * 265 * @return the value estimated from the current illumination of the 266 * microsphere. 267 */ 268 private double interpolate() { 269 // Number of non-illuminated facets. 270 int darkCount = 0; 271 272 double value = 0; 273 double totalWeight = 0; 274 for (FacetData fd : microsphereData) { 275 final double iV = fd.illumination(); 276 if (iV != 0d) { 277 value += iV * fd.sample(); 278 totalWeight += iV; 279 } else { 280 ++darkCount; 281 } 282 } 283 284 final double darkFraction = darkCount / (double) size; 285 286 return darkFraction <= maxDarkFraction ? 287 value / totalWeight : 288 background; 289 } 290 291 /** 292 * Illumination. 293 * 294 * @param sampleDirection Vector whose origin is at the interpolation 295 * point and tail is at the sample location. 296 * @param sampleValue Data value of the sample. 297 * @param weight Weight. 298 */ 299 private void illuminate(double[] sampleDirection, 300 double sampleValue, 301 double weight) { 302 for (int i = 0; i < size; i++) { 303 final double[] n = microsphere.get(i).getNormal(); 304 final double cos = MathArrays.cosAngle(n, sampleDirection); 305 306 if (cos > 0) { 307 final double illumination = cos * weight; 308 309 if (illumination > darkThreshold && 310 illumination > microsphereData.get(i).illumination()) { 311 microsphereData.set(i, new FacetData(illumination, sampleValue)); 312 } 313 } 314 } 315 } 316 317 /** 318 * Reset the all the {@link Facet facets} data to zero. 319 */ 320 private void clear() { 321 for (int i = 0; i < size; i++) { 322 microsphereData.set(i, new FacetData(0d, 0d)); 323 } 324 } 325 326 /** 327 * Microsphere "facet" (surface element). 328 */ 329 private static class Facet { 330 /** Normal vector characterizing a surface element. */ 331 private final double[] normal; 332 333 /** 334 * @param n Normal vector characterizing a surface element 335 * of the microsphere. No copy is made. 336 */ 337 Facet(double[] n) { 338 normal = n; 339 } 340 341 /** 342 * Return a reference to the vector normal to this facet. 343 * 344 * @return the normal vector. 345 */ 346 public double[] getNormal() { 347 return normal; 348 } 349 } 350 351 /** 352 * Data associated with each {@link Facet}. 353 */ 354 private static class FacetData { 355 /** Illumination received from the sample. */ 356 private final double illumination; 357 /** Data value of the sample. */ 358 private final double sample; 359 360 /** 361 * @param illumination Illumination. 362 * @param sample Data value. 363 */ 364 FacetData(double illumination, double sample) { 365 this.illumination = illumination; 366 this.sample = sample; 367 } 368 369 /** 370 * Get the illumination. 371 * @return the illumination. 372 */ 373 public double illumination() { 374 return illumination; 375 } 376 377 /** 378 * Get the data value. 379 * @return the data value. 380 */ 381 public double sample() { 382 return sample; 383 } 384 } 385}