1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.rng.sampling.shape;
19
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.SharedStateObjectSampler;
22
23 /**
24 * Generate points <a href="https://mathworld.wolfram.com/TrianglePointPicking.html">
25 * uniformly distributed within a triangle</a>.
26 *
27 * <ul>
28 * <li>
29 * Uses the algorithm described in:
30 * <blockquote>
31 * Turk, G. <i>Generating random points in triangles</i>. Glassner, A. S. (ed) (1990).<br>
32 * <strong>Graphic Gems</strong> Academic Press, pp. 24-28.
33 * </blockquote>
34 * </li>
35 * </ul>
36 *
37 * <p>Sampling uses:</p>
38 *
39 * <ul>
40 * <li>{@link UniformRandomProvider#nextDouble()}
41 * </ul>
42 *
43 * @since 1.4
44 */
45 public abstract class TriangleSampler implements SharedStateObjectSampler<double[]> {
46 /** The dimension for 2D sampling. */
47 private static final int TWO_D = 2;
48 /** The dimension for 3D sampling. */
49 private static final int THREE_D = 3;
50 /** The source of randomness. */
51 private final UniformRandomProvider rng;
52
53 // The following code defines a plane as the set of all points r:
54 // r = r0 + sv + tw
55 // where s and t range over all real numbers, v and w are given linearly independent
56 // vectors defining the plane, and r0 is an arbitrary (but fixed) point in the plane.
57 //
58 // Sampling from a triangle (a,b,c) is performed when:
59 // s and t are in [0, 1] and s+t <= 1;
60 // r0 is one triangle vertex (a);
61 // and v (b-a) and w (c-a) are vectors from the other two vertices to r0.
62 //
63 // For robustness with large value coordinates the point r is computed without
64 // the requirement to compute v and w which can overflow:
65 //
66 // a + s(b-a) + t(c-a) == a + sb - sa + tc - ta
67 // == a(1 - s - t) + sb + tc
68 //
69 // Assuming the uniform deviates are from the 2^53 dyadic rationals in [0, 1) if s+t <= 1
70 // then 1 - (s+t) is exact. Sampling is then done using:
71 //
72 // if (s + t <= 1):
73 // p = a(1 - (s + t)) + sb + tc
74 // else:
75 // p = a(1 - (1 - s) - (1 - t)) + (1 - s)b + (1 - t)c
76 // p = a(s - 1 + t) + (1 - s)b + (1 - t)c
77 //
78 // Note do not simplify (1 - (1 - s) - (1 - t)) to (s + t - 1) as s+t > 1 and has potential
79 // loss of a single bit of randomness due to rounding. An exact sum is s - 1 + t.
80
81 /**
82 * Sample uniformly from a triangle in 2D. This is an non-array based specialisation of
83 * {@link TriangleSamplerND} for performance.
84 */
85 private static class TriangleSampler2D extends TriangleSampler {
86 /** The x component of vertex a. */
87 private final double ax;
88 /** The y component of vertex a. */
89 private final double ay;
90 /** The x component of vertex b. */
91 private final double bx;
92 /** The y component of vertex b. */
93 private final double by;
94 /** The x component of vertex c. */
95 private final double cx;
96 /** The y component of vertex c. */
97 private final double cy;
98
99 /**
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 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 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 }