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. & 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}