commit 6ff800d59f21902c621b04a86d2ab7ad0b7bc4bc
parent 958cefa602645282b4ac22fcd1fc2f689a6c7f37
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 24 Sep 2019 16:59:06 +0200
Add uniform sampling in a tetrahedron
Diffstat:
5 files changed, 348 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -180,6 +180,7 @@ if(NOT NO_TEST)
new_test(test_ssp_ran_sphere ${MATH_LIB})
new_test(test_ssp_ran_sphere_cap ${MATH_LIB})
new_test(test_ssp_ran_triangle ${MATH_LIB})
+ new_test(test_ssp_ran_tetrahedron ${MATH_LIB})
new_test(test_ssp_ran_uniform_disk)
rcmake_copy_runtime_libraries(test_ssp_rng)
endif()
diff --git a/src/ssp.h b/src/ssp.h
@@ -465,6 +465,63 @@ ssp_ran_triangle_uniform_float_pdf
}
/*******************************************************************************
+ * Tetrahedron distribution
+ ******************************************************************************/
+ /* Uniform sampling of a tetrahedron */
+SSP_API double*
+ssp_ran_tetrahedron_uniform
+ (struct ssp_rng* rng,
+ const double v0[3], /* Position of the 1st vertex */
+ const double v1[3], /* Position of the 2nd vertex */
+ const double v2[3], /* Position of the 3rd vertex */
+ const double v3[3], /* Position of the 4th vertex */
+ double sample[3], /* Sampled position */
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_tetrahedron_uniform_float
+ (struct ssp_rng* rng,
+ const float v0[3], /* Position of the 1st vertex */
+ const float v1[3], /* Position of the 2nd vertex */
+ const float v2[3], /* Position of the 3rd vertex */
+ const float v3[3], /* Position of the 4th vertex */
+ float sample[3], /* Sampled position */
+ float* pdf); /* Can be NULL => pdf not computed */
+
+/* Return the probability distribution for a uniform point sampling in a tetrahedron */
+static INLINE double
+ssp_ran_tetrahedron_uniform_pdf
+ (const double v0[3],
+ const double v1[3],
+ const double v2[3],
+ const double v3[3])
+{
+ double vec0[3], vec1[3], vec2[3], tmp[3];
+ ASSERT(v0 && v1 && v2 && v3);
+
+ d3_sub(vec0, v1, v2);
+ d3_sub(vec1, v3, v2);
+ d3_sub(vec2, v0, v2);
+ return 6 / fabs(d3_dot(d3_cross(tmp, vec0, vec1), vec2));
+}
+
+static INLINE float
+ssp_ran_tetrahedron_uniform_float_pdf
+ (const float v0[3],
+ const float v1[3],
+ const float v2[3],
+ const float v3[3])
+{
+ float vec0[3], vec1[3], vec2[3], tmp[3];
+ ASSERT(v0 && v1 && v2 && v3);
+
+ f3_sub(vec0, v1, v2);
+ f3_sub(vec1, v3, v2);
+ f3_sub(vec2, v0, v2);
+ return 6 / fabsf(f3_dot(f3_cross(tmp, vec0, vec1), vec2));
+}
+
+/*******************************************************************************
* Hemisphere distribution
******************************************************************************/
/* Uniform sampling of an unit hemisphere whose up direction is implicitly the
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -291,6 +291,94 @@ ssp_ran_triangle_uniform_float
}
double*
+ssp_ran_tetrahedron_uniform
+ (struct ssp_rng* rng,
+ const double v0[3],
+ const double v1[3],
+ const double v2[3],
+ const double v3[3],
+ double sample[3],
+ double* pdf)
+{
+ double a, s, t, u; /* Barycentric coordinates of the sampled point. */
+ double tmp0[3], tmp1[3], tmp2[3], tmp3[3], tmp4[3], tmp5[3];
+ ASSERT(rng && v0 && v1 && v2 && v3 && sample);
+
+ s = ssp_rng_canonical(rng);
+ t = ssp_rng_canonical(rng);
+ u = ssp_rng_canonical(rng);
+
+ if (s + t > 1) { /* Cut and fold the cube into a prism */
+ s = 1 - s;
+ t = 1 - t;
+ }
+ if (t + u > 1) { /* Cut and fold the prism into a tetrahedron */
+ double swp = u;
+ u = 1 - s - t;
+ t = 1 - swp;
+ }
+ else if (s + t + u > 1) {
+ double swp = u;
+ u = s + t + u - 1;
+ s = 1 - t - swp;
+ }
+ a = 1 - s - t - u;
+
+ d3_add(sample,
+ d3_add(tmp4, d3_muld(tmp0, v0, a), d3_muld(tmp1, v1, s)),
+ d3_add(tmp5, d3_muld(tmp2, v2, t), d3_muld(tmp3, v3, u)));
+
+ if (pdf)
+ *pdf = ssp_ran_tetrahedron_uniform_pdf(v0, v1, v2, v3);
+
+ return sample;
+}
+
+float*
+ssp_ran_tetrahedron_uniform_float
+ (struct ssp_rng* rng,
+ const float v0[3],
+ const float v1[3],
+ const float v2[3],
+ const float v3[3],
+ float sample[3],
+ float* pdf)
+{
+ float a, s, t, u; /* Barycentric coordinates of the sampled point. */
+ float tmp0[3], tmp1[3], tmp2[3], tmp3[3], tmp4[3], tmp5[3];
+ ASSERT(rng && v0 && v1 && v2 && v3 && sample);
+
+ s = ssp_rng_canonical_float(rng);
+ t = ssp_rng_canonical_float(rng);
+ u = ssp_rng_canonical_float(rng);
+
+ if (s + t > 1) { /* Cut and fold the cube into a prism */
+ s = 1 - s;
+ t = 1 - t;
+ }
+ if (t + u > 1) { /* Cut and fold the prism into a tetrahedron */
+ float swp = u;
+ u = 1 - s - t;
+ t = 1 - swp;
+ }
+ else if (s + t + u > 1) {
+ float swp = u;
+ u = s + t + u - 1;
+ s = 1 - t - swp;
+ }
+ a = 1 - s - t - u;
+
+ f3_add(sample,
+ f3_add(tmp4, f3_mulf(tmp0, v0, a), f3_mulf(tmp1, v1, s)),
+ f3_add(tmp5, f3_mulf(tmp2, v2, t), f3_mulf(tmp3, v3, u)));
+
+ if (pdf)
+ *pdf = ssp_ran_tetrahedron_uniform_float_pdf(v0, v1, v2, v3);
+
+ return sample;
+}
+
+double*
ssp_ran_hemisphere_uniform_local
(struct ssp_rng* rng,
double sample[3],
diff --git a/src/test_ssp_ran_tetrahedron.c b/src/test_ssp_ran_tetrahedron.c
@@ -0,0 +1,42 @@
+/* Copyright (C) |Meso|Star> 2015-2018 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#define TYPE_FLOAT 1
+#include "test_ssp_ran_tetrahedron.h"
+
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_tetrahedron.h"
+
+int
+main(int argc, char** argv)
+{
+ (void)argc, (void)argv;
+ test_float();
+ test_double();
+ return 0;
+}
diff --git a/src/test_ssp_ran_tetrahedron.h b/src/test_ssp_ran_tetrahedron.h
@@ -0,0 +1,160 @@
+/* Copyright (C) |Meso|Star> 2015-2018 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#ifndef TEST_SSP_RAN_TETRA_H
+#define TEST_SSP_RAN_TETRA_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#include <rsys/float4.h>
+
+#define NSAMPS 128
+
+#endif /* TEST_SSP_RAN_TETRA_H */
+
+#if TYPE_FLOAT==0
+ #define REAL double
+ #define TEST test_double
+ #define RNG_UNIFORM_R ssp_rng_uniform_double
+ #define RAN_TETRA_UNIFORM ssp_ran_tetrahedron_uniform
+ #define RAN_TETRA_UNIFORM_PDF ssp_ran_tetrahedron_uniform_pdf
+ #define EQ_EPS_R eq_eps
+ #define FABS_R fabs
+ #define R3 d3
+ #define R3_DOT d3_dot
+ #define R3_SUB d3_sub
+ #define R3_MINUS d3_minus
+ #define R3_CROSS d3_cross
+ #define R3_LEN d3_len
+
+#elif TYPE_FLOAT==1
+ #define REAL float
+ #define TEST test_float
+ #define RNG_UNIFORM_R ssp_rng_uniform_float
+ #define RAN_TETRA_UNIFORM ssp_ran_tetrahedron_uniform_float
+ #define RAN_TETRA_UNIFORM_PDF ssp_ran_tetrahedron_uniform_float_pdf
+ #define EQ_EPS_R eq_epsf
+ #define FABS_R fabsf
+ #define R3 f3
+ #define R3_DOT f3_dot
+ #define R3_SUB f3_sub
+ #define R3_MINUS f3_minus
+ #define R3_CROSS f3_cross
+ #define R3_LEN f3_len
+
+#else
+ #error "TYPE_FLOAT must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ REAL samps[NSAMPS][3];
+ REAL A[3], B[3], C[3], D[3];
+ REAL v0[3], v1[3], v2[3], v3[3], v4[3];
+ REAL in0[3], in1[3], in2[3], in3[3];
+ REAL pdf[NSAMPS];
+ size_t i, j;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng) == RES_OK);
+
+ FOR_EACH(j, 0, 100) {
+ FOR_EACH(i, 0, 3) {
+ A[i] = RNG_UNIFORM_R(rng, 0, 1);
+ B[i] = RNG_UNIFORM_R(rng, 0, 1);
+ C[i] = RNG_UNIFORM_R(rng, 0, 1);
+ D[i] = RNG_UNIFORM_R(rng, 0, 1);
+ }
+
+ R3_SUB(v0, B, A);
+ R3_SUB(v1, C, A);
+ R3_SUB(v2, C, B);
+ R3_SUB(v3, D, A);
+ R3_SUB(v4, D, C);
+ /* Inward vectors perpandicular to faces */
+ R3_CROSS(in0, v0, v1);
+ if (R3_DOT(in0, v3) < 0) R3_MINUS(in0, in0);
+ R3_CROSS(in1, v1, v3);
+ if (R3_DOT(in1, v0) < 0) R3_MINUS(in1, in1);
+ R3_CROSS(in2, v3, v0);
+ if (R3_DOT(in2, v1) < 0) R3_MINUS(in2, in2);
+ R3_CROSS(in3, v2, v4);
+ if (R3_DOT(in3, v3) > 0) R3_MINUS(in3, in3);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ REAL X[3], tmp[3];
+ REAL dot = 0;
+ REAL vol = 0;
+
+ CHK(RAN_TETRA_UNIFORM(rng, A, B, C, D, samps[i], &pdf[i]) == samps[i]);
+
+ /* Check sample is in the tetrahedron */
+ R3_SUB(X, samps[i], A);
+ dot = R3_DOT(X, in0);
+ CHK(sign(dot) == 1);
+ dot = R3_DOT(X, in1);
+ CHK(sign(dot) == 1);
+ dot = R3_DOT(X, in2);
+ CHK(sign(dot) == 1);
+ dot = R3_DOT(X, in3);
+ CHK(sign(dot) == -1);
+
+ /* Check pdf */
+ vol = FABS_R(R3_DOT(R3_CROSS(tmp, v0, v1), v3)) / 6;
+ CHK(EQ_EPS_R(pdf[i], RAN_TETRA_UNIFORM_PDF(A, B, C, D), (REAL)1.e-8) == 1);
+ /* Pdf is 1/vol
+ * But numerical value depends on the vector used to compute it */
+ CHK(EQ_EPS_R(vol * pdf[i], 1, (REAL)1.e-5) == 1);
+ }
+ }
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RNG_UNIFORM_R
+#undef RAN_TETRA_UNIFORM
+#undef RAN_TETRA_UNIFORM_PDF
+#undef EQ_EPS_R
+#undef FABS_R
+#undef R3
+#undef R3_DOT
+#undef R3_SUB
+#undef R3_MINUS
+#undef R3_CROSS
+#undef R3_LEN
+#undef TYPE_FLOAT
