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}