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 uniformly distributed within a
025 * <a href="https://en.wikipedia.org/wiki/Tetrahedron">tetrahedron</a>.
026 *
027 * <ul>
028 *  <li>
029 *   Uses the algorithm described in:
030 *   <blockquote>
031 *    Rocchini, C. and Cignoni, P. (2001)<br>
032 *    <i>Generating Random Points in a Tetrahedron</i>.<br>
033 *    <strong>Journal of Graphics Tools</strong> 5(4), pp. 9-12.
034 *   </blockquote>
035 *  </li>
036 * </ul>
037 *
038 * <p>Sampling uses:</p>
039 *
040 * <ul>
041 *   <li>{@link UniformRandomProvider#nextDouble()}
042 * </ul>
043 *
044 * @see <a href="https://doi.org/10.1080/10867651.2000.10487528">
045 *   Rocchini, C. &amp; Cignoni, P. (2001) Journal of Graphics Tools 5, pp. 9-12</a>
046 * @since 1.4
047 */
048public class TetrahedronSampler implements SharedStateObjectSampler<double[]> {
049    /** The dimension for 3D sampling. */
050    private static final int THREE_D = 3;
051    /** The name of vertex a. */
052    private static final String VERTEX_A = "Vertex a";
053    /** The name of vertex b. */
054    private static final String VERTEX_B = "Vertex b";
055    /** The name of vertex c. */
056    private static final String VERTEX_C = "Vertex c";
057    /** The name of vertex d. */
058    private static final String VERTEX_D = "Vertex d";
059
060    /** The first vertex. */
061    private final double[] a;
062    /** The second vertex. */
063    private final double[] b;
064    /** The third vertex. */
065    private final double[] c;
066    /** The fourth vertex. */
067    private final double[] d;
068    /** The source of randomness. */
069    private final UniformRandomProvider rng;
070
071    /**
072     * @param rng Source of randomness.
073     * @param a The first vertex.
074     * @param b The second vertex.
075     * @param c The third vertex.
076     * @param d The fourth vertex.
077     */
078    TetrahedronSampler(UniformRandomProvider rng, double[] a, double[] b, double[] c, double[] d) {
079        // Defensive copy
080        this.a = a.clone();
081        this.b = b.clone();
082        this.c = c.clone();
083        this.d = d.clone();
084        this.rng = rng;
085    }
086
087    /**
088     * @param rng Generator of uniformly distributed random numbers
089     * @param source Source to copy.
090     */
091    TetrahedronSampler(UniformRandomProvider rng, TetrahedronSampler source) {
092        // Shared state is immutable
093        a = source.a;
094        b = source.b;
095        c = source.c;
096        d = source.d;
097        this.rng = rng;
098    }
099
100    /**
101     * @return a random Cartesian point within the tetrahedron.
102     */
103    @Override
104    public double[] sample() {
105        double s = rng.nextDouble();
106        double t = rng.nextDouble();
107        final double u = rng.nextDouble();
108        // Care is taken to ensure the 3 deviates remain in the 2^53 dyadic rationals in [0, 1).
109        // The following are exact for all the 2^53 dyadic rationals:
110        // 1 - u; u in [0, 1]
111        // u - 1; u in [0, 1]
112        // u + 1; u in [-1, 0]
113        // u + v; u in [-1, 0], v in [0, 1]
114        // u + v; u, v in [0, 1], u + v <= 1
115
116        // Cut and fold with the plane s + t = 1
117        if (s + t > 1) {
118            // (s, t, u) = (1 - s, 1 - t, u)                if s + t > 1
119            s = 1 - s;
120            t = 1 - t;
121        }
122        // Now s + t <= 1.
123        // Cut and fold with the planes t + u = 1 and s + t + u = 1.
124        final double tpu = t + u;
125        final double sptpu = s + tpu;
126        if (sptpu > 1) {
127            if (tpu > 1) {
128                // (s, t, u) = (s, 1 - u, 1 - s - t)        if t + u > 1
129                // 1 - s - (1-u) - (1-s-t) == u - 1 + t
130                return createSample(u - 1 + t, s, 1 - u, 1 - s - t);
131            }
132            // (s, t, u) = (1 - t - u, t, s + t + u - 1)    if t + u <= 1
133            // 1 - (1-t-u) - t - (s+t+u-1) == 1 - s - t
134            return createSample(1 - s - t, 1 - tpu, t, s - 1 + tpu);
135        }
136        return createSample(1 - sptpu, s, t, u);
137    }
138
139    /**
140     * Creates the sample given the random variates {@code s}, {@code t} and {@code u} in the
141     * interval {@code [0, 1]} and {@code s + t + u <= 1}. The sum {@code 1 - s - t - u} is
142     * provided. The sample can be obtained from the tetrahedron {@code abcd} using:
143     *
144     * <pre>
145     * p = (1 - s - t - u)a + sb + tc + ud
146     * </pre>
147     *
148     * @param p1msmtmu plus 1 minus s minus t minus u ({@code 1 - s - t - u})
149     * @param s the first variate s
150     * @param t the second variate t
151     * @param u the third variate u
152     * @return the sample
153     */
154    private double[] createSample(double p1msmtmu, double s, double t, double u) {
155        // From the barycentric coordinates s,t,u create the point by moving along
156        // vectors ab, ac and ad.
157        // Here we do not compute the vectors and use the original vertices:
158        // p = a + s(b-a) + t(c-a) + u(d-a)
159        //   = (1-s-t-u)a + sb + tc + ud
160        return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0],
161                             p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1],
162                             p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]};
163    }
164
165    /** {@inheritDoc} */
166    @Override
167    public TetrahedronSampler withUniformRandomProvider(UniformRandomProvider rng) {
168        return new TetrahedronSampler(rng, this);
169    }
170
171    /**
172     * Create a tetrahedron sampler with vertices {@code a}, {@code b}, {@code c} and {@code d}.
173     * Sampled points are uniformly distributed within the tetrahedron.
174     *
175     * <p>No test for a volume is performed. If the vertices are coplanar the sampling
176     * distribution is undefined.
177     *
178     * @param rng Source of randomness.
179     * @param a The first vertex.
180     * @param b The second vertex.
181     * @param c The third vertex.
182     * @param d The fourth vertex.
183     * @return the sampler
184     * @throws IllegalArgumentException If the vertices do not have length 3;
185     * or vertices have non-finite coordinates
186     */
187    public static TetrahedronSampler of(UniformRandomProvider rng,
188                                        double[] a,
189                                        double[] b,
190                                        double[] c,
191                                        double[] d) {
192        // Must be 3D
193        Coordinates.requireLength(a, THREE_D, VERTEX_A);
194        Coordinates.requireLength(b, THREE_D, VERTEX_B);
195        Coordinates.requireLength(c, THREE_D, VERTEX_C);
196        Coordinates.requireLength(d, THREE_D, VERTEX_D);
197        // Detect non-finite vertices
198        Coordinates.requireFinite(a, VERTEX_A);
199        Coordinates.requireFinite(b, VERTEX_B);
200        Coordinates.requireFinite(c, VERTEX_C);
201        Coordinates.requireFinite(d, VERTEX_D);
202        return new TetrahedronSampler(rng, a, b, c, d);
203    }
204}