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