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 018package org.apache.commons.rng.sampling.shape; 019 020import org.apache.commons.rng.UniformRandomProvider; 021import org.apache.commons.rng.sampling.SharedStateObjectSampler; 022 023/** 024 * Generate points <a href="https://mathworld.wolfram.com/TrianglePointPicking.html"> 025 * uniformly distributed within a triangle</a>. 026 * 027 * <ul> 028 * <li> 029 * Uses the algorithm described in: 030 * <blockquote> 031 * Turk, G. <i>Generating random points in triangles</i>. Glassner, A. S. (ed) (1990).<br> 032 * <strong>Graphic Gems</strong> Academic Press, pp. 24-28. 033 * </blockquote> 034 * </li> 035 * </ul> 036 * 037 * <p>Sampling uses:</p> 038 * 039 * <ul> 040 * <li>{@link UniformRandomProvider#nextDouble()} 041 * </ul> 042 * 043 * @since 1.4 044 */ 045public abstract class TriangleSampler implements SharedStateObjectSampler<double[]> { 046 /** The dimension for 2D sampling. */ 047 private static final int TWO_D = 2; 048 /** The dimension for 3D sampling. */ 049 private static final int THREE_D = 3; 050 /** The source of randomness. */ 051 private final UniformRandomProvider rng; 052 053 // The following code defines a plane as the set of all points r: 054 // r = r0 + sv + tw 055 // where s and t range over all real numbers, v and w are given linearly independent 056 // vectors defining the plane, and r0 is an arbitrary (but fixed) point in the plane. 057 // 058 // Sampling from a triangle (a,b,c) is performed when: 059 // s and t are in [0, 1] and s+t <= 1; 060 // r0 is one triangle vertex (a); 061 // and v (b-a) and w (c-a) are vectors from the other two vertices to r0. 062 // 063 // For robustness with large value coordinates the point r is computed without 064 // the requirement to compute v and w which can overflow: 065 // 066 // a + s(b-a) + t(c-a) == a + sb - sa + tc - ta 067 // == a(1 - s - t) + sb + tc 068 // 069 // Assuming the uniform deviates are from the 2^53 dyadic rationals in [0, 1) if s+t <= 1 070 // then 1 - (s+t) is exact. Sampling is then done using: 071 // 072 // if (s + t <= 1): 073 // p = a(1 - (s + t)) + sb + tc 074 // else: 075 // p = a(1 - (1 - s) - (1 - t)) + (1 - s)b + (1 - t)c 076 // p = a(s - 1 + t) + (1 - s)b + (1 - t)c 077 // 078 // Note do not simplify (1 - (1 - s) - (1 - t)) to (s + t - 1) as s+t > 1 and has potential 079 // loss of a single bit of randomness due to rounding. An exact sum is s - 1 + t. 080 081 /** 082 * Sample uniformly from a triangle in 2D. This is an non-array based specialisation of 083 * {@link TriangleSamplerND} for performance. 084 */ 085 private static final class TriangleSampler2D extends TriangleSampler { 086 /** The x component of vertex a. */ 087 private final double ax; 088 /** The y component of vertex a. */ 089 private final double ay; 090 /** The x component of vertex b. */ 091 private final double bx; 092 /** The y component of vertex b. */ 093 private final double by; 094 /** The x component of vertex c. */ 095 private final double cx; 096 /** The y component of vertex c. */ 097 private final double cy; 098 099 /** 100 * @param rng Source of randomness. 101 * @param a The first vertex. 102 * @param b The second vertex. 103 * @param c The third vertex. 104 */ 105 TriangleSampler2D(UniformRandomProvider rng, double[] a, double[] b, double[] c) { 106 super(rng); 107 ax = a[0]; 108 ay = a[1]; 109 bx = b[0]; 110 by = b[1]; 111 cx = c[0]; 112 cy = c[1]; 113 } 114 115 /** 116 * @param rng Generator of uniformly distributed random numbers 117 * @param source Source to copy. 118 */ 119 TriangleSampler2D(UniformRandomProvider rng, TriangleSampler2D source) { 120 super(rng); 121 ax = source.ax; 122 ay = source.ay; 123 bx = source.bx; 124 by = source.by; 125 cx = source.cx; 126 cy = source.cy; 127 } 128 129 @Override 130 public double[] createSample(double p1msmt, double s, double t) { 131 return new double[] {p1msmt * ax + s * bx + t * cx, 132 p1msmt * ay + s * by + t * cy}; 133 } 134 135 @Override 136 public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) { 137 return new TriangleSampler2D(rng, this); 138 } 139 } 140 141 /** 142 * Sample uniformly from a triangle in 3D. This is an non-array based specialisation of 143 * {@link TriangleSamplerND} for performance. 144 */ 145 private static final class TriangleSampler3D extends TriangleSampler { 146 /** The x component of vertex a. */ 147 private final double ax; 148 /** The y component of vertex a. */ 149 private final double ay; 150 /** The z component of vertex a. */ 151 private final double az; 152 /** The x component of vertex b. */ 153 private final double bx; 154 /** The y component of vertex b. */ 155 private final double by; 156 /** The z component of vertex b. */ 157 private final double bz; 158 /** The x component of vertex c. */ 159 private final double cx; 160 /** The y component of vertex c. */ 161 private final double cy; 162 /** The z component of vertex c. */ 163 private final double cz; 164 165 /** 166 * @param rng Source of randomness. 167 * @param a The first vertex. 168 * @param b The second vertex. 169 * @param c The third vertex. 170 */ 171 TriangleSampler3D(UniformRandomProvider rng, double[] a, double[] b, double[] c) { 172 super(rng); 173 ax = a[0]; 174 ay = a[1]; 175 az = a[2]; 176 bx = b[0]; 177 by = b[1]; 178 bz = b[2]; 179 cx = c[0]; 180 cy = c[1]; 181 cz = c[2]; 182 } 183 184 /** 185 * @param rng Generator of uniformly distributed random numbers 186 * @param source Source to copy. 187 */ 188 TriangleSampler3D(UniformRandomProvider rng, TriangleSampler3D source) { 189 super(rng); 190 ax = source.ax; 191 ay = source.ay; 192 az = source.az; 193 bx = source.bx; 194 by = source.by; 195 bz = source.bz; 196 cx = source.cx; 197 cy = source.cy; 198 cz = source.cz; 199 } 200 201 @Override 202 public double[] createSample(double p1msmt, double s, double t) { 203 return new double[] {p1msmt * ax + s * bx + t * cx, 204 p1msmt * ay + s * by + t * cy, 205 p1msmt * az + s * bz + t * cz}; 206 } 207 208 @Override 209 public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) { 210 return new TriangleSampler3D(rng, this); 211 } 212 } 213 214 /** 215 * Sample uniformly from a triangle in ND. 216 */ 217 private static final class TriangleSamplerND extends TriangleSampler { 218 /** The first vertex. */ 219 private final double[] a; 220 /** The second vertex. */ 221 private final double[] b; 222 /** The third vertex. */ 223 private final double[] c; 224 225 /** 226 * @param rng Source of randomness. 227 * @param a The first vertex. 228 * @param b The second vertex. 229 * @param c The third vertex. 230 */ 231 TriangleSamplerND(UniformRandomProvider rng, double[] a, double[] b, double[] c) { 232 super(rng); 233 // Defensive copy 234 this.a = a.clone(); 235 this.b = b.clone(); 236 this.c = c.clone(); 237 } 238 239 /** 240 * @param rng Generator of uniformly distributed random numbers 241 * @param source Source to copy. 242 */ 243 TriangleSamplerND(UniformRandomProvider rng, TriangleSamplerND source) { 244 super(rng); 245 // Shared state is immutable 246 a = source.a; 247 b = source.b; 248 c = source.c; 249 } 250 251 @Override 252 public double[] createSample(double p1msmt, double s, double t) { 253 final double[] x = new double[a.length]; 254 for (int i = 0; i < x.length; i++) { 255 x[i] = p1msmt * a[i] + s * b[i] + t * c[i]; 256 } 257 return x; 258 } 259 260 @Override 261 public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) { 262 return new TriangleSamplerND(rng, this); 263 } 264 } 265 266 /** 267 * @param rng Source of randomness. 268 */ 269 TriangleSampler(UniformRandomProvider rng) { 270 this.rng = rng; 271 } 272 273 /** 274 * @return a random Cartesian coordinate within the triangle. 275 */ 276 @Override 277 public double[] sample() { 278 final double s = rng.nextDouble(); 279 final double t = rng.nextDouble(); 280 final double spt = s + t; 281 if (spt > 1) { 282 // Transform: s1 = 1 - s; t1 = 1 - t. 283 // Compute: 1 - s1 - t1 284 // Do not assume (1 - (1-s) - (1-t)) is (s + t - 1), i.e. (spt - 1.0), 285 // to avoid loss of a random bit due to rounding when s + t > 1. 286 // An exact sum is (s - 1 + t). 287 return createSample(s - 1.0 + t, 1.0 - s, 1.0 - t); 288 } 289 // Here s + t is exact so can be subtracted to make 1 - s - t 290 return createSample(1.0 - spt, s, t); 291 } 292 293 /** 294 * Creates the sample given the random variates {@code s} and {@code t} in the 295 * interval {@code [0, 1]} and {@code s + t <= 1}. The sum {@code 1 - s - t} is provided. 296 * The sample can be obtained from the triangle abc using: 297 * <pre> 298 * p = a(1 - s - t) + sb + tc 299 * </pre> 300 * 301 * @param p1msmt plus 1 minus s minus t (1 - s - t) 302 * @param s the first variate s 303 * @param t the second variate t 304 * @return the sample 305 */ 306 protected abstract double[] createSample(double p1msmt, double s, double t); 307 308 /** {@inheritDoc} */ 309 // Redeclare the signature to return a TriangleSampler not a SharedStateObjectSampler<double[]> 310 @Override 311 public abstract TriangleSampler withUniformRandomProvider(UniformRandomProvider rng); 312 313 /** 314 * Create a triangle sampler with vertices {@code a}, {@code b} and {@code c}. 315 * Sampled points are uniformly distributed within the triangle. 316 * 317 * <p>Sampling is supported in dimensions of 2 or above. Samples will lie in the 318 * plane (2D Euclidean space) defined by using the three triangle vertices to 319 * create two vectors starting at a point in the plane and orientated in 320 * different directions along the plane. 321 * 322 * <p>No test for collinear points is performed. If the points are collinear the sampling 323 * distribution is undefined. 324 * 325 * @param rng Source of randomness. 326 * @param a The first vertex. 327 * @param b The second vertex. 328 * @param c The third vertex. 329 * @return the sampler 330 * @throws IllegalArgumentException If the vertices do not have the same 331 * dimension; the dimension is less than 2; or vertices have non-finite coordinates 332 */ 333 public static TriangleSampler of(UniformRandomProvider rng, 334 double[] a, 335 double[] b, 336 double[] c) { 337 final int dimension = a.length; 338 if (dimension != b.length || dimension != c.length) { 339 throw new IllegalArgumentException( 340 new StringBuilder("Mismatch of vertex dimensions: ").append(dimension).append(',') 341 .append(b.length).append(',') 342 .append(c.length).toString()); 343 } 344 // Detect non-finite vertices 345 Coordinates.requireFinite(a, "Vertex a"); 346 Coordinates.requireFinite(b, "Vertex b"); 347 Coordinates.requireFinite(c, "Vertex c"); 348 // Low dimension specialisations 349 if (dimension == TWO_D) { 350 return new TriangleSampler2D(rng, a, b, c); 351 } else if (dimension == THREE_D) { 352 return new TriangleSampler3D(rng, a, b, c); 353 } else if (dimension > THREE_D) { 354 return new TriangleSamplerND(rng, a, b, c); 355 } 356 // Less than 2D 357 throw new IllegalArgumentException( 358 "Unsupported dimension: " + dimension); 359 } 360}