commit cc3fc4074314b7f03b350f3f3e744ed9f93b048b
parent 1aa2de5d96259fa4324bf8f4fc99782a5aaa345d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Sat, 4 Apr 2015 09:08:36 +0200
Add and test the <uniform|cosine weighted> hemispheric random variates
Diffstat:
3 files changed, 254 insertions(+), 1 deletion(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -94,13 +94,15 @@ if(NOT NO_TEST)
endfunction()
function(new_test _name)
- build_test(_name)
+ build_test(${_name} ${ARGN})
add_test(${_name} ${_name})
endfunction()
build_test(test_ssp_rng)
add_test(test_ssp_rng_kiss test_ssp_rng kiss)
add_test(test_ssp_rng_mt19937_64 test_ssp_rng mt19937_64)
+
+ new_test(test_ssp_ran_hemisphere m)
endif()
################################################################################
diff --git a/src/ssp.h b/src/ssp.h
@@ -32,6 +32,8 @@
#ifndef SSP_H
#define SSP_H
+#include <rsys/float33.h>
+#include <rsys/math.h>
#include <rsys/rsys.h>
/* Library symbol management */
@@ -137,6 +139,106 @@ SSP_API unsigned long
ssp_rng_max
(struct ssp_rng* rng);
+/*******************************************************************************
+ * Hemisphere distribution
+ ******************************************************************************/
+/* Uniform sampling of an unit hemisphere whose up direction is implicitly the
+ * Z axis. The PDF of the generated sample is stored in sample[3] */
+static INLINE float*
+ssp_ran_hemisphere_uniform_local(struct ssp_rng* rng, float sample[4])
+{
+ double phi, cos_theta, sin_theta;
+ ASSERT(rng && sample);
+ phi = 2.0 * PI * ssp_rng_canonical(rng);
+ cos_theta = ssp_rng_canonical(rng);
+ sin_theta = cos2sin(cos_theta);
+ sample[0] = (float)(cos(phi) * sin_theta);
+ sample[1] = (float)(sin(phi) * sin_theta);
+ sample[2] = (float)cos_theta;
+ sample[3] = (float)(1.0/(2.0*PI));
+ return sample;
+}
+
+/* Return the probability distribution for an hemispheric uniform random
+ * variate with an implicit up direction in Z */
+static INLINE float
+ssp_ran_hemisphere_uniform_local_pdf(const float sample[3])
+{
+ ASSERT(sample);
+ return sample[2] < 0.f ? 0.f : (float)(1.0/(2.0*PI));
+}
+
+/* Uniform sampling of an unit hemisphere whose up direction is `up'. The PDF
+ * of the generated sample is stored in sample[3] */
+static INLINE float*
+ssp_ran_hemisphere_uniform
+ (struct ssp_rng* rng, const float up[3], float sample[4])
+{
+ float basis[9];
+ ASSERT(rng && up && sample && f3_is_normalized(up));
+ ssp_ran_hemisphere_uniform_local(rng, sample);
+ return f33_mulf3(sample, f33_basis(basis, up), sample);
+}
+
+/* Return the probability distribution for an hemispheric uniform random
+ * variate with an explicit `up' direction */
+static INLINE float
+ssp_ran_hemisphere_uniform_pdf(const float up[3], const float sample[3])
+{
+ ASSERT(up && sample && f3_is_normalized(up));
+ return f3_dot(sample, up) < 0.f ? 0.f : (float)(1.0/(2.0*PI));
+}
+
+/* Cosine weighted sampling of an unit hemisphere whise up direction is
+ * implicitly the Z axis. The PDF of the generated sample is stored in
+ * sample[3] */
+static INLINE float*
+ssp_ran_hemisphere_cos_local(struct ssp_rng* rng, float sample[4])
+{
+ double phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = 2.0 * PI * ssp_rng_canonical(rng);
+ v = ssp_rng_canonical(rng);
+ cos_theta = sqrt(v);
+ sin_theta = sqrt(1.0 - v);
+ sample[0] = (float)(cos(phi) * sin_theta);
+ sample[1] = (float)(sin(phi) * sin_theta);
+ sample[2] = (float)cos_theta;
+ sample[3] = (float)(cos_theta / PI);
+ return sample;
+}
+
+/* Return the probability distribution for an hemispheric consine weighted
+ * random variate with an implicit up direction in Z */
+static INLINE float
+ssp_ran_hemisphere_cos_local_pdf(const float sample[3])
+{
+ ASSERT(sample);
+ return sample[2] < 0.f ? 0.f : (float)(sample[2] / PI);
+}
+
+/* Cosine weighted sampling of an unit hemisphere whose up direction is `up'.
+ * The PDF of the generated sample is stored in sample[3] */
+static INLINE float*
+ssp_ran_hemisphere_cos(struct ssp_rng* rng, const float up[3], float sample[4])
+{
+ float basis[9];
+ ASSERT(rng && up && sample && f3_is_normalized(up));
+ ssp_ran_hemisphere_cos_local(rng, sample);
+ return f33_mulf3(sample, f33_basis(basis, up), sample);
+}
+
+/* Return the probability distribution for an hemispheric consine weighted
+ * random variate with an explicit `up' direction */
+static INLINE float
+ssp_ran_hemisphere_cos_pdf(const float up[3], const float sample[3])
+{
+ float dot;
+ ASSERT(up && sample);
+ dot = f3_dot(sample, up);
+ return dot < 0.f ? 0.f : (float)(dot/PI);
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/test_ssp_ran_hemisphere.c b/src/test_ssp_ran_hemisphere.c
@@ -0,0 +1,149 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is a program whose purpose is to test the spp library.
+ *
+ * This software is governed by the CeCILL-C 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-C
+ * 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-C license and that you accept its terms. */
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#include <rsys/float4.h>
+
+#define NSAMPS 128
+
+int
+main(int argc, char** argv)
+{
+ struct ssp_rng* rng0, *rng1;
+ struct mem_allocator allocator;
+ float samps0[NSAMPS][4];
+ float samps1[NSAMPS][4];
+ float samps2[NSAMPS][4];
+ float samps3[NSAMPS][4];
+ int i = 0, j = 0;
+ float* f;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng0), RES_OK);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng1), RES_OK);
+
+ ssp_rng_set(rng0, 0);
+ FOR_EACH(i, 0, NSAMPS) {
+ float frame[9];
+ float up[3] = {0.f, 0.f, 1.f};
+ float xyz[3];
+ unsigned long seed = ssp_rng_get(rng0);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_uniform_local(rng1, samps0[i]);
+ CHECK(f, samps0[i]);
+ CHECK(f3_is_normalized(f), 1);
+ CHECK(eq_epsf(f[3], (float)(1.0/(2.0*PI)), 1.e-6f), 1);
+ CHECK(eq_epsf(f[3], ssp_ran_hemisphere_uniform_local_pdf(f), 1.e-6f), 1);
+ CHECK(f3_dot(f, up) >= 0.f, 1);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_uniform(rng1, up, samps1[i]);
+ CHECK(f, samps1[i]);
+ CHECK(f4_eq_eps(f, samps0[i], 1.e-6f), 1);
+
+ up[0] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[1] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[2] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ f3_normalize(up, up);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_uniform(rng1, up, samps1[i]);
+ CHECK(f3_eq_eps(samps0[i], samps1[i], 1.e-6f), 0);
+ CHECK(f3_is_normalized(f), 1);
+ CHECK(f3_dot(f, up) >= 0.f, 1);
+ CHECK(eq_epsf(f[3], (float)(1.0/(2.0*PI)), 1.e-6f ), 1);
+ CHECK(eq_epsf(f[3], ssp_ran_hemisphere_uniform_pdf(up, f), 1.e-6f), 1);
+
+ f33_basis(frame, up);
+ f33_mulf3(xyz, frame, samps0[i]);
+ CHECK(f3_eq_eps(samps1[i], xyz, 1.e-6f), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(f3_eq_eps(samps0[i], samps0[j], 1.e-6f), 0);
+ CHECK(f3_eq_eps(samps1[i], samps1[j], 1.e-6f), 0);
+ }
+ }
+
+ ssp_rng_set(rng0, 0);
+ FOR_EACH(i, 0, NSAMPS) {
+ float frame[9];
+ float up[3] = { 0.f, 0.f, 1.f };
+ float xyz[3];
+ unsigned long seed = ssp_rng_get(rng0);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_cos_local(rng1, samps2[i]);
+ CHECK(f, samps2[i]);
+ CHECK(f3_eq_eps(samps0[i], samps2[i], 1.e-6f), 0);
+ CHECK(f3_is_normalized(f), 1);
+ CHECK(eq_epsf(f[3], (float)(f[2]/PI), 1.e-6f), 1);
+ CHECK(eq_epsf(f[3], ssp_ran_hemisphere_cos_local_pdf(f), 1.e-6f), 1);
+ CHECK(f3_dot(f, up) >= 0.f, 1);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_cos(rng1, up, samps3[i]);
+ CHECK(f, samps3[i]);
+ CHECK(f4_eq_eps(f, samps2[i], 1.e-6f), 1);
+
+ up[0] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[1] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[2] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ f3_normalize(up, up);
+
+ ssp_rng_set(rng1, seed);
+ f = ssp_ran_hemisphere_cos(rng1, up, samps3[i]);
+ CHECK(f3_eq_eps(samps2[i], samps3[i], 1.e-6f), 0);
+ CHECK(f3_eq_eps(samps1[i], samps3[i], 1.e-6f), 0);
+ CHECK(f3_is_normalized(f), 1);
+ CHECK(f3_dot(f, up) >= 0.f, 1);
+ CHECK(eq_epsf(f[3], (float)(f3_dot(f, up)/PI), 1.e-6f), 1);
+ CHECK(eq_epsf(f[3], ssp_ran_hemisphere_cos_pdf(f, up), 1.e-6f), 1);
+
+ f33_basis(frame, up);
+ f33_mulf3(xyz, frame, samps2[i]);
+ CHECK(f3_eq_eps(samps3[i], xyz, 1.e-6f), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(f3_eq_eps(samps2[i], samps2[j], 1.e-6f), 0);
+ CHECK(f3_eq_eps(samps3[i], samps3[j], 1.e-6f), 0);
+ }
+ }
+
+ ssp_rng_ref_put(rng0);
+ ssp_rng_ref_put(rng1);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
