commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TetrahedronSampler.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.rng.sampling.shape;

import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.SharedStateSampler;

/**
 * Generate points uniformly distributed within a
 * <a href="https://en.wikipedia.org/wiki/Tetrahedron">tetrahedron</a>.
 *
 * <ul>
 *  <li>
 *   Uses the algorithm described in:
 *   <blockquote>
 *    Rocchini, C. and Cignoni, P. (2001)<br>
 *    <i>Generating Random Points in a Tetrahedron</i>.<br>
 *    <strong>Journal of Graphics Tools</strong> 5(4), pp. 9-12.
 *   </blockquote>
 *  </li>
 * </ul>
 *
 * <p>Sampling uses:</p>
 *
 * <ul>
 *   <li>{@link UniformRandomProvider#nextDouble()}
 * </ul>
 *
 * @see <a href="https://doi.org/10.1080/10867651.2000.10487528">
 *   Rocchini, C. &amp; Cignoni, P. (2001) Journal of Graphics Tools 5, pp. 9-12</a>
 * @since 1.4
 */
public class TetrahedronSampler implements SharedStateSampler<TetrahedronSampler> {
    /** The dimension for 3D sampling. */
    private static final int THREE_D = 3;
    /** The name of vertex a. */
    private static final String VERTEX_A = "Vertex a";
    /** The name of vertex b. */
    private static final String VERTEX_B = "Vertex b";
    /** The name of vertex c. */
    private static final String VERTEX_C = "Vertex c";
    /** The name of vertex d. */
    private static final String VERTEX_D = "Vertex d";

    /** The first vertex. */
    private final double[] a;
    /** The second vertex. */
    private final double[] b;
    /** The third vertex. */
    private final double[] c;
    /** The fourth vertex. */
    private final double[] d;
    /** The source of randomness. */
    private final UniformRandomProvider rng;

    /**
     * @param a The first vertex.
     * @param b The second vertex.
     * @param c The third vertex.
     * @param d The fourth vertex.
     * @param rng Source of randomness.
     */
    TetrahedronSampler(double[] a, double[] b, double[] c, double[] d, UniformRandomProvider rng) {
        // Defensive copy
        this.a = a.clone();
        this.b = b.clone();
        this.c = c.clone();
        this.d = d.clone();
        this.rng = rng;
    }

    /**
     * @param rng Generator of uniformly distributed random numbers
     * @param source Source to copy.
     */
    TetrahedronSampler(UniformRandomProvider rng, TetrahedronSampler source) {
        // Shared state is immutable
        a = source.a;
        b = source.b;
        c = source.c;
        d = source.d;
        this.rng = rng;
    }

    /**
     * @return a random Cartesian point within the tetrahedron.
     */
    public double[] sample() {
        double s = rng.nextDouble();
        double t = rng.nextDouble();
        final double u = rng.nextDouble();
        // Care is taken to ensure the 3 deviates remain in the 2^53 dyadic rationals in [0, 1).
        // The following are exact for all the 2^53 dyadic rationals:
        // 1 - u; u in [0, 1]
        // u - 1; u in [0, 1]
        // u + 1; u in [-1, 0]
        // u + v; u in [-1, 0], v in [0, 1]
        // u + v; u, v in [0, 1], u + v <= 1

        // Cut and fold with the plane s + t = 1
        if (s + t > 1) {
            // (s, t, u) = (1 - s, 1 - t, u)                if s + t > 1
            s = 1 - s;
            t = 1 - t;
        }
        // Now s + t <= 1.
        // Cut and fold with the planes t + u = 1 and s + t + u = 1.
        final double tpu = t + u;
        final double sptpu = s + tpu;
        if (sptpu > 1) {
            if (tpu > 1) {
                // (s, t, u) = (s, 1 - u, 1 - s - t)        if t + u > 1
                // 1 - s - (1-u) - (1-s-t) == u - 1 + t
                return createSample(u - 1 + t, s, 1 - u, 1 - s - t);
            }
            // (s, t, u) = (1 - t - u, t, s + t + u - 1)    if t + u <= 1
            // 1 - (1-t-u) - t - (s+t+u-1) == 1 - s - t
            return createSample(1 - s - t, 1 - tpu, t, s - 1 + tpu);
        }
        return createSample(1 - sptpu, s, t, u);
    }

    /**
     * Creates the sample given the random variates {@code s}, {@code t} and {@code u} in the
     * interval {@code [0, 1]} and {@code s + t + u <= 1}. The sum {@code 1 - s - t - u} is
     * provided. The sample can be obtained from the tetrahedron {@code abcd} using:
     *
     * <pre>
     * p = (1 - s - t - u)a + sb + tc + ud
     * </pre>
     *
     * @param p1msmtmu plus 1 minus s minus t minus u ({@code 1 - s - t - u})
     * @param s the first variate s
     * @param t the second variate t
     * @param u the third variate u
     * @return the sample
     */
    private double[] createSample(double p1msmtmu, double s, double t, double u) {
        // From the barycentric coordinates s,t,u create the point by moving along
        // vectors ab, ac and ad.
        // Here we do not compute the vectors and use the original vertices:
        // p = a + s(b-a) + t(c-a) + u(d-a)
        //   = (1-s-t-u)a + sb + tc + ud
        return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0],
                             p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1],
                             p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]};
    }

    /** {@inheritDoc} */
    @Override
    public TetrahedronSampler withUniformRandomProvider(UniformRandomProvider rng) {
        return new TetrahedronSampler(rng, this);
    }

    /**
     * Create a tetrahedron sampler with vertices {@code a}, {@code b}, {@code c} and {@code d}.
     * Sampled points are uniformly distributed within the tetrahedron.
     *
     * <p>No test for a volume is performed. If the vertices are coplanar the sampling
     * distribution is undefined.
     *
     * @param a The first vertex.
     * @param b The second vertex.
     * @param c The third vertex.
     * @param d The fourth vertex.
     * @param rng Source of randomness.
     * @return the sampler
     * @throws IllegalArgumentException If the vertices do not have length 3;
     * or vertices have non-finite coordinates
     */
    public static TetrahedronSampler of(double[] a,
                                        double[] b,
                                        double[] c,
                                        double[] d,
                                        UniformRandomProvider rng) {
        // Must be 3D
        Coordinates.requireLength(a, THREE_D, VERTEX_A);
        Coordinates.requireLength(b, THREE_D, VERTEX_B);
        Coordinates.requireLength(c, THREE_D, VERTEX_C);
        Coordinates.requireLength(d, THREE_D, VERTEX_D);
        // Detect non-finite vertices
        Coordinates.requireFinite(a, VERTEX_A);
        Coordinates.requireFinite(b, VERTEX_B);
        Coordinates.requireFinite(c, VERTEX_C);
        Coordinates.requireFinite(d, VERTEX_D);
        return new TetrahedronSampler(a, b, c, d, rng);
    }
}