commit b7fa585173d41ed54d66c19bfc88ab125b35ded1
parent 1c8fb36dad8b23f69b6677803013e4f894d5a006
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 9 Jul 2021 11:05:13 +0200
Add uniform sampling on a spherical zone
Diffstat:
5 files changed, 350 insertions(+), 1 deletion(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -179,6 +179,7 @@ if(NOT NO_TEST)
new_test(test_ssp_rng_proxy)
new_test(test_ssp_ran_sphere ${MATH_LIB})
new_test(test_ssp_ran_sphere_cap ${MATH_LIB})
+ new_test(test_ssp_ran_spherical_zone ${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)
diff --git a/src/ssp.h b/src/ssp.h
@@ -1049,7 +1049,7 @@ ssp_ran_uniform_disk_float
/*******************************************************************************
* Uniform distribution of a point ont a sphere cap
- * pdf = 1/(2*PI*(1-cos(theta_max))
+ * pdf = 1/(2*PI*(1-cos(aperture))
******************************************************************************/
/* Uniform sampling of unit sphere cap centered in zero whose up direction
* is implicity the Z axis. */
@@ -1115,6 +1115,81 @@ ssp_ran_sphere_cap_uniform_float
return f33_mulf3(sample, f33_basis(basis, up), sample_local);
}
+/*******************************************************************************
+ * Uniform distribution of a point ont a truncated sphere cap
+ * pdf = 1/(2*PI*(cos(aperture_max)-cos(aperture_min))
+ ******************************************************************************/
+/* Uniform sampling of unit sphere cap centered in zero whose up direction
+ * is implicity the Z axis. */
+SSP_API double*
+ssp_ran_spherical_zone_uniform_local
+ (struct ssp_rng* rng,
+ const double height_range[2],
+ double pt[3],
+ double* pdf);
+
+SSP_API float*
+ssp_ran_spherical_zone_uniform_float_local
+ (struct ssp_rng* rng,
+ const float height_range[2],
+ float pt[3],
+ float* pdf);
+
+/* Uniform sampling of unit truncated sphere cap centered in zero whose
+ * up direction is explicitly defined */
+static INLINE double*
+ssp_ran_spherical_zone_uniform
+ (struct ssp_rng* rng,
+ const double up[3],
+ const double height_range[2],
+ double sample[3],
+ double* pdf) /* Can be NULL */
+{
+ double sample_local[3];
+ double basis[9];
+ ASSERT(up);
+ ssp_ran_spherical_zone_uniform_local(rng, height_range, sample_local, pdf);
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
+}
+
+static INLINE float*
+ssp_ran_spherical_zone_uniform_float
+ (struct ssp_rng* rng,
+ const float up[3],
+ const float height_range[2],
+ float sample[3],
+ float* pdf) /* Can be NULL */
+{
+ float sample_local[3];
+ float basis[9];
+ ASSERT(up);
+ ssp_ran_spherical_zone_uniform_float_local(rng, height_range, sample_local, pdf);
+ return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+}
+
+static INLINE double
+ssp_ran_spherical_zone_uniform_pdf(const double height_range[2])
+{
+ ASSERT(height_range && height_range[0] <= height_range[1]
+ && -1 <= height_range[0] && height_range[1] <= 1);
+ if(height_range[0] == height_range[1]) {
+ return height_range[0] == 1 ? INF : ssp_ran_circle_uniform_pdf();
+ }
+ return 1.0/(2.0*PI*(height_range[1]-height_range[0]));
+}
+
+static INLINE float
+ssp_ran_spherical_zone_uniform_float_pdf(const float height_range[2])
+{
+ ASSERT(height_range && height_range[0] <= height_range[1]
+ && -1 <= height_range[0] && height_range[1] <= 1);
+ if(height_range[0] == height_range[1]) {
+ return height_range[0] == 1 ?
+ (float)INF : ssp_ran_circle_uniform_float_pdf();
+ }
+ return 1.f/(2.f*(float)PI*(height_range[1]-height_range[0]));
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -721,3 +721,73 @@ ssp_ran_sphere_cap_uniform_float_local
return pt;
}
+double*
+ssp_ran_spherical_zone_uniform_local
+ (struct ssp_rng* rng,
+ const double height_range[2],
+ double pt[3],
+ double* pdf)
+{
+ ASSERT(rng && pt && height_range
+ && height_range[0] <= height_range[1]
+ && height_range[0] >= -1 && height_range[1] <= 1);
+
+ if(height_range[1] == 1) {
+ ssp_ran_sphere_cap_uniform_local(rng, height_range[0], pt, pdf);
+ }
+ else if(height_range[0] == height_range[1]) {
+ double sin_theta = cos2sin(height_range[0]);
+ ssp_ran_circle_uniform(rng, pt, pdf);
+ pt[0] *= sin_theta;
+ pt[1] *= sin_theta;
+ pt[2] = height_range[0];
+ } else {
+ const double cos_aperture_min = height_range[0];
+ const double cos_aperture_max = height_range[1];
+ double phi, cos_theta, sin_theta;
+ phi = ssp_rng_uniform_double(rng, 0, 2.0*PI);
+ cos_theta = ssp_rng_uniform_double(rng, cos_aperture_min, cos_aperture_max);
+ sin_theta = cos2sin(cos_theta);
+ pt[0] = cos(phi) * sin_theta;
+ pt[1] = sin(phi) * sin_theta;
+ pt[2] = cos_theta;
+ if(pdf) *pdf = 1.0/(2.0*PI*(cos_aperture_max-cos_aperture_min));
+ }
+ return pt;
+}
+
+float*
+ssp_ran_spherical_zone_uniform_float_local
+ (struct ssp_rng* rng,
+ const float height_range[2],
+ float pt[3],
+ float* pdf)
+{
+ ASSERT(rng && pt && height_range
+ && height_range[0] <= height_range[1]
+ && height_range[0] >= -1 && height_range[1] <= 1);
+
+ if(height_range[1] == 1) {
+ ssp_ran_sphere_cap_uniform_float_local(rng, height_range[0], pt, pdf);
+ }
+ else if(height_range[0] == height_range[1]) {
+ float sin_theta = (float)cos2sin(height_range[0]);
+ ssp_ran_circle_uniform_float(rng, pt, pdf);
+ pt[0] *= sin_theta;
+ pt[1] *= sin_theta;
+ pt[2] = height_range[0];
+ } else {
+ const float cos_aperture_min = height_range[0];
+ const float cos_aperture_max = height_range[1];
+ float phi, cos_theta, sin_theta;
+ phi = ssp_rng_uniform_float(rng, 0, 2.f*(float)PI);
+ cos_theta = ssp_rng_uniform_float(rng, cos_aperture_min, cos_aperture_max);
+ sin_theta = (float)cos2sin(cos_theta);
+ pt[0] = cosf(phi) * sin_theta;
+ pt[1] = sinf(phi) * sin_theta;
+ pt[2] = cos_theta;
+ if(pdf) *pdf = 1.f/(2.f*(float)PI*(cos_aperture_max-cos_aperture_min));
+ }
+ return pt;
+}
+
diff --git a/src/test_ssp_ran_spherical_zone.c b/src/test_ssp_ran_spherical_zone.c
@@ -0,0 +1,43 @@
+/* Copyright (C) 2015-2021 |Meso|Star> (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 0
+#include "test_ssp_ran_spherical_zone.h"
+
+#define TYPE_FLOAT 1
+#include "test_ssp_ran_spherical_zone.h"
+
+int
+main(int argc, char** argv)
+{
+ (void)argc, (void)argv;
+ test_float();
+ test_double();
+ return 0;
+}
+
diff --git a/src/test_ssp_ran_spherical_zone.h b/src/test_ssp_ran_spherical_zone.h
@@ -0,0 +1,160 @@
+/* Copyright (C) 2015-2021 |Meso|Star> (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_SPHERICAL_ZONE_H
+#define TEST_SSP_RAN_SPHERICAL_ZONE_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NSAMPS 1024
+
+#endif /* TEST_SSP_RAN_SPHERE_H */
+
+#if TYPE_FLOAT==0
+ #define REAL double
+ #define TEST test_double
+ #define RAN_SPHERICAL_ZONE_UNIFORM_LOCAL ssp_ran_spherical_zone_uniform_local
+ #define RAN_SPHERICAL_ZONE_UNIFORM_PDF ssp_ran_spherical_zone_uniform_pdf
+ #define RAN_SPHERICAL_ZONE_UNIFORM ssp_ran_spherical_zone_uniform
+ #define RAN_SPHERE_UNIFORM ssp_ran_sphere_uniform
+ #define EQ_EPS eq_eps
+ #define R3_EQ_EPS d3_eq_eps
+ #define R3_IS_NORMALIZED d3_is_normalized
+ #define R3_NORMALIZE d3_normalize
+ #define R3_DOT d3_dot
+#elif TYPE_FLOAT==1
+ #define REAL float
+ #define TEST test_float
+ #define RAN_SPHERICAL_ZONE_UNIFORM_LOCAL ssp_ran_spherical_zone_uniform_float_local
+ #define RAN_SPHERICAL_ZONE_UNIFORM_PDF ssp_ran_spherical_zone_uniform_float_pdf
+ #define RAN_SPHERICAL_ZONE_UNIFORM ssp_ran_spherical_zone_uniform_float
+ #define RAN_SPHERE_UNIFORM ssp_ran_sphere_uniform_float
+ #define EQ_EPS eq_epsf
+ #define R3_EQ_EPS f3_eq_eps
+ #define R3_IS_NORMALIZED f3_is_normalized
+ #define R3_NORMALIZE f3_normalize
+ #define R3_DOT f3_dot
+#else
+ #error "TYPE_FLOAT must be defined either 0 or 1"
+#endif
+
+static void
+TEST(void)
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ REAL pdf;
+ REAL samps[NSAMPS][3];
+ REAL up[3];
+ REAL heights[2];
+ int i, j;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+ CHK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng) == RES_OK);
+
+ heights[0] = heights[1] = 1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM_LOCAL(rng, heights, samps[0], &pdf) == samps[0]);
+ CHK(samps[0][0] == 0);
+ CHK(samps[0][1] == 0);
+ CHK(samps[0][2] == 1);
+ CHK(IS_INF(pdf));
+
+ heights[0] = -1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM_LOCAL(rng, heights, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(4*(REAL)PI), (REAL)1.e-6));
+
+ heights[0] = heights[1] = (REAL)0.4;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM_LOCAL(rng, heights, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI), (REAL)1.e-6));
+
+ /* Check NULL pdf */
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM_LOCAL(rng, heights, samps[0], NULL) == samps[0]);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ heights[0] = (REAL)-0.7; heights[1] = 1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM_LOCAL(rng, heights, samps[i], &pdf) == samps[i]);
+ CHK(R3_IS_NORMALIZED(samps[i]));
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI*(heights[1]-heights[0])), (REAL)1.e-6));
+ CHK(EQ_EPS(pdf, RAN_SPHERICAL_ZONE_UNIFORM_PDF(heights), (REAL)1.e-6));
+ CHK(samps[i][2] >= heights[0] && samps[i][2] <= heights[1]);
+ FOR_EACH(j, 0, i) {
+ CHK(!R3_EQ_EPS(samps[j], samps[i], (REAL)1.e-6));
+ }
+ }
+
+ /* Sample an up vector */
+ RAN_SPHERE_UNIFORM(rng, up, NULL);
+
+ heights[0] = heights[1] = 1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM(rng, up, heights, samps[0], &pdf) == samps[0]);
+ CHK(R3_EQ_EPS(samps[0], up, (REAL)1.e-6));
+ CHK(IS_INF(pdf));
+
+ heights[0] = -1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM(rng, up, heights, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(4*(REAL)PI), (REAL)1.e-6));
+
+ heights[0] = heights[1] = (REAL)0.9;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM(rng, up, heights, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI), (REAL)1.e-6));
+
+ /* Check NULL pdf */
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM(rng, up, heights, samps[0], NULL) == samps[0]);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ heights[0] = (REAL)0.3; heights[1] = 1;
+ CHK(RAN_SPHERICAL_ZONE_UNIFORM(rng, up, heights, samps[i], &pdf) == samps[i]);
+ CHK(R3_IS_NORMALIZED(samps[i]));
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI*(heights[1]-heights[0])), (REAL)1.e-6));
+ CHK(EQ_EPS(pdf, RAN_SPHERICAL_ZONE_UNIFORM_PDF(heights), (REAL)1.e-6));
+ CHK(R3_DOT(up, samps[i]) >= heights[0] && R3_DOT(up, samps[i]) <= heights[1]);
+ FOR_EACH(j, 0, i) {
+ CHK(!R3_EQ_EPS(samps[j], samps[i], (REAL)1.e-6));
+ }
+ }
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHK(mem_allocated_size() == 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RAN_SPHERICAL_ZONE_UNIFORM_LOCAL
+#undef RAN_SPHERICAL_ZONE_UNIFORM_PDF
+#undef RAN_SPHERICAL_ZONE_UNIFORM
+#undef RAN_SPHERE_UNIFORM
+#undef EQ_EPS
+#undef R3_EQ_EPS
+#undef R3_IS_NORMALIZED
+#undef R3_NORMALIZE
+#undef R3_DOT
+#undef TYPE_FLOAT
