commit 70f003a188d55e51a16e4b6ce97a16b8f45e8c19
parent 33a383a157ae3108e766d6bd79d39964dbd7b147
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 8 Oct 2018 17:24:10 +0200
Merge branch 'release_0.8'
Diffstat:
7 files changed, 347 insertions(+), 14 deletions(-)
diff --git a/README.md b/README.md
@@ -34,10 +34,18 @@ variable to the directory that contains the `boost` include directory.
## Release notes
+### Version 0.8
+
+- Add the `ssp_ran_sphere_cap_uniform[_local][_float]` random variates that
+ uniformly distributes a point onto a unit sphere cap centered in zero.
+- Fix the allocation of the `ssp_rng` data structure: the `ssp_rng` structure
+ contains a long double attribute that might be not correctly aligned on 16
+ bytes.
+
### Version 0.7
-- Add the `ssp_ran_circle_uniform` random variate that uniformly distribute a 2
- dimensional position onto a unit circle.
+- Add the `ssp_ran_circle_uniform` random variate that uniformly distributes a
+ 2 dimensional position onto a unit circle.
### Version 0.6
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,7 +29,7 @@
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.0)
project(star-sp C CXX)
enable_testing()
@@ -54,7 +54,7 @@ include(rcmake_runtime)
if(MSVC)
# Currently, the C++11 random library on MSVC is bugged and consequently we
# use the Boost library instead
- find_package(Boost 1.58 REQUIRED COMPONENTS random)
+ find_package(Boost 1.67 REQUIRED COMPONENTS random)
add_definitions(-DUSE_BOOST_RANDOM)
include_directories(${RSys_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${Random123_INCLUDE_DIR})
else()
@@ -71,7 +71,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys ${Boost_LIBRARY_DIRS})
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 7)
+set(VERSION_MINOR 8)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -98,6 +98,12 @@ rcmake_prepend_path(SSP_FILES_INC_API ${SSP_SOURCE_DIR})
rcmake_prepend_path(SSP_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
set_source_files_properties(${SSP_FILES_SRC} PROPERTIES LANGUAGE CXX)
+if(CMAKE_COMPILER_IS_GNUCXX)
+ # Random123 warning that cannot be disabled by a GCC pragma
+ set_source_files_properties(${SSP_SOURCE_DIR}/ssp_rng.c PROPERTIES
+ COMPILE_FLAGS "-Wno-expansion-to-defined")
+endif()
+
add_library(ssp SHARED ${SSP_FILES_SRC} ${SSP_FILES_INC} ${SSP_FILES_INC_API})
set_target_properties(ssp PROPERTIES
DEFINE_SYMBOL SSP_SHARED_BUILD
@@ -161,15 +167,16 @@ if(NOT NO_TEST)
if(BUILD_R123_AES)
register_test(test_ssp_rng_aes test_ssp_rng aes)
endif()
+ new_test(test_ssp_ran_circle ${MATH_LIB})
new_test(test_ssp_ran_discrete)
+ new_test(test_ssp_ran_gaussian)
new_test(test_ssp_ran_hemisphere ${MATH_LIB})
- new_test(test_ssp_ran_sphere ${MATH_LIB})
- new_test(test_ssp_ran_circle ${MATH_LIB})
new_test(test_ssp_ran_hg ${MATH_LIB})
- new_test(test_ssp_ran_triangle ${MATH_LIB})
- new_test(test_ssp_rng_proxy)
- new_test(test_ssp_ran_gaussian)
new_test(test_ssp_ran_piecewise_linear)
+ 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_triangle ${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
@@ -74,7 +74,7 @@ struct ssp_rng_type {
res_T (*set)(void* state, const uint64_t seed);
uint64_t (*get)(void* state);
- res_T(*discard)(void* state, uint64_t n);
+ res_T (*discard)(void* state, uint64_t n);
res_T (*read)(void* state, FILE* file);
res_T (*write)(const void* state, FILE* file);
double (*entropy)(const void* state);
@@ -303,7 +303,7 @@ ssp_ran_exp_float_pdf
/* Stateless random variate from the gaussian (or normal) distribution
* with mean `mu':
- * P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(1/2*((x-mu)/sigma)^2) dx */
+ * P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(-1/2*((x-mu)/sigma)^2) dx */
SSP_API double
ssp_ran_gaussian
(struct ssp_rng* rng,
@@ -927,6 +927,74 @@ ssp_ran_uniform_disk_float
return f33_mulf3(pt, f33_basis(basis, up), sample_local);
}
+/*******************************************************************************
+ * Uniform distribution of a point ont a sphere cap
+ * pdf = 1/(2*PI*(1-cos(theta_max))
+ ******************************************************************************/
+/* Uniform sampling of unit sphere cap centered in zero whose up direction
+ * is implicity the Z axis. */
+SSP_API double*
+ssp_ran_sphere_cap_uniform_local
+ (struct ssp_rng* rng,
+ /* In [-1, 1]. Height where the sphere cap begins. It equal cos(theta_max)
+ * where theta_max is the aperture angle from the up axis */
+ const double height,
+ double sample[3],
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_sphere_cap_uniform_float_local
+ (struct ssp_rng* rng,
+ const float height, /* In [-1, 1]. Height of the sphere cap. */
+ float sample[3],
+ float* pdf); /* Can be NULL => pdf not computed */
+
+static INLINE double
+ssp_ran_sphere_cap_uniform_pdf(const double height)
+{
+ ASSERT(height <= 1.0 && height >= -1.0);
+ return height == 1.0 ? INF : 1.0/(2.0*PI*(1-height));
+}
+
+static INLINE float
+ssp_ran_sphere_cap_uniform_float_pdf(const float height)
+{
+ ASSERT(height <= 1.f && height >= -1.f);
+ return height == 1.f ? (float)INF : 1.f/(2.f*(float)PI*(1.f-height));
+}
+
+/* Uniform sampling of unit sphere cap centered in zero whose up direction
+ * is explicitly defined */
+static INLINE double*
+ssp_ran_sphere_cap_uniform
+ (struct ssp_rng* rng,
+ const double up[3],
+ const double height,
+ double sample[3],
+ double* pdf) /* Can be NULL */
+{
+ double sample_local[3];
+ double basis[9];
+ ASSERT(up);
+ ssp_ran_sphere_cap_uniform_local(rng, height, sample_local, pdf);
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
+}
+
+static INLINE float*
+ssp_ran_sphere_cap_uniform_float
+ (struct ssp_rng* rng,
+ const float up[3],
+ const float height,
+ float sample[3],
+ float* pdf)
+{
+ float sample_local[3];
+ float basis[9];
+ ASSERT(up);
+ ssp_ran_sphere_cap_uniform_float_local(rng, height, sample_local, pdf);
+ return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -219,7 +219,7 @@ ssp_ran_circle_uniform
return sample;
}
-float*
+float*
ssp_ran_circle_uniform_float
(struct ssp_rng* rng,
float sample[2],
@@ -521,3 +521,59 @@ ssp_ran_uniform_disk_float_local
if (pdf) *pdf = 1 / (radius * radius);
return pt;
}
+
+double*
+ssp_ran_sphere_cap_uniform_local
+ (struct ssp_rng* rng,
+ const double height,
+ double pt[3],
+ double* pdf)
+{
+ ASSERT(rng && pt && (height >= -1 || height <= 1));
+
+ if(height == 1) {
+ d3(pt, 0, 0, 1);
+ if(pdf) *pdf = INF;
+ } else if(height == -1) {
+ ssp_ran_sphere_uniform(rng, pt, pdf);
+ } else {
+ const double cos_aperture = height;
+ 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, 1);
+ 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*(1.0-cos_aperture));
+ }
+ return pt;
+}
+
+float*
+ssp_ran_sphere_cap_uniform_float_local
+ (struct ssp_rng* rng,
+ const float height,
+ float pt[3],
+ float* pdf)
+{
+ ASSERT(rng && pt && (height >= -1 || height <= 1));
+ if(height == 1) {
+ f3(pt, 0, 0, 1);
+ if(pdf) *pdf = (float)INF;
+ } else if(height == -1) {
+ ssp_ran_sphere_uniform_float(rng, pt, pdf);
+ } else {
+ const float cos_aperture = height;
+ 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, 1);
+ 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*(1.f-cos_aperture));
+ }
+ return pt;
+}
+
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -479,11 +479,15 @@ ssp_rng_create
}
allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
- rng = (struct ssp_rng*)MEM_CALLOC(allocator, 1, sizeof(struct ssp_rng));
+
+ /* Align the rng on 16 bytes since its 'r' attrib is a long double that on
+ * Linux 64-bits systems must be aligned on 16 bytes. */
+ rng = (struct ssp_rng*)MEM_ALLOC_ALIGNED(allocator, sizeof(*rng), 16);
if(!rng) {
res = RES_MEM_ERR;
goto error;
}
+ memset(rng, 0, sizeof(struct ssp_rng));
rng->allocator = allocator;
ref_init(&rng->ref);
diff --git a/src/test_ssp_ran_sphere_cap.c b/src/test_ssp_ran_sphere_cap.c
@@ -0,0 +1,43 @@
+/* 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 0
+#include "test_ssp_ran_sphere_cap.h"
+
+#define TYPE_FLOAT 1
+#include "test_ssp_ran_sphere_cap.h"
+
+int
+main(int argc, char** argv)
+{
+ (void)argc, (void)argv;
+ test_float();
+ test_double();
+ return 0;
+}
+
diff --git a/src/test_ssp_ran_sphere_cap.h b/src/test_ssp_ran_sphere_cap.h
@@ -0,0 +1,147 @@
+/* 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_SPHERE_CAP_H
+#define TEST_SSP_RAN_SPHERE_CAP_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_SPHERE_CAP_UNIFORM_LOCAL ssp_ran_sphere_cap_uniform_local
+ #define RAN_SPHERE_CAP_UNIFORM_PDF ssp_ran_sphere_cap_uniform_pdf
+ #define RAN_SPHERE_CAP_UNIFORM ssp_ran_sphere_cap_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_SPHERE_CAP_UNIFORM_LOCAL ssp_ran_sphere_cap_uniform_float_local
+ #define RAN_SPHERE_CAP_UNIFORM_PDF ssp_ran_sphere_cap_uniform_float_pdf
+ #define RAN_SPHERE_CAP_UNIFORM ssp_ran_sphere_cap_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];
+ int i, j;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+ CHK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng) == RES_OK);
+
+ CHK(RAN_SPHERE_CAP_UNIFORM_LOCAL(rng, 1, samps[0], &pdf) == samps[0]);
+ CHK(samps[0][0] == 0);
+ CHK(samps[0][1] == 0);
+ CHK(samps[0][2] == 1);
+ CHK(IS_INF(pdf));
+
+ CHK(RAN_SPHERE_CAP_UNIFORM_LOCAL(rng,-1, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(4*(REAL)PI), (REAL)1.e-6));
+
+ /* Check NULL pdf */
+ CHK(RAN_SPHERE_CAP_UNIFORM_LOCAL(rng, (REAL)0.2, samps[0], NULL) == samps[0]);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ const REAL height = (REAL)-0.7;
+ CHK(RAN_SPHERE_CAP_UNIFORM_LOCAL(rng, height, samps[i], &pdf) == samps[i]);
+ CHK(R3_IS_NORMALIZED(samps[i]));
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI*(1-height)), (REAL)1.e-6));
+ CHK(EQ_EPS(pdf, RAN_SPHERE_CAP_UNIFORM_PDF(height), (REAL)1.e-6));
+ CHK(samps[i][2] >= height);
+ 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);
+
+ CHK(RAN_SPHERE_CAP_UNIFORM(rng, up, 1, samps[0], &pdf) == samps[0]);
+ CHK(R3_EQ_EPS(samps[0], up, (REAL)1.e-6));
+ CHK(IS_INF(pdf));
+
+ CHK(RAN_SPHERE_CAP_UNIFORM(rng, up, -1, samps[0], &pdf) == samps[0]);
+ CHK(EQ_EPS(pdf, 1/(4*(REAL)PI), (REAL)1.e-6));
+
+ /* Check NULL pdf */
+ CHK(RAN_SPHERE_CAP_UNIFORM(rng, up, (REAL)0.2, samps[0], NULL) == samps[0]);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ const REAL height = (REAL)0.3;
+ CHK(RAN_SPHERE_CAP_UNIFORM(rng, up, height, samps[i], &pdf) == samps[i]);
+ CHK(R3_IS_NORMALIZED(samps[i]));
+ CHK(EQ_EPS(pdf, 1/(2*(REAL)PI*(1-height)), (REAL)1.e-6));
+ CHK(EQ_EPS(pdf, RAN_SPHERE_CAP_UNIFORM_PDF(height), (REAL)1.e-6));
+ CHK(R3_DOT(up, samps[i]) >= height);
+ 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_SPHERE_CAP_UNIFORM_LOCAL
+#undef RAN_SPHERE_CAP_UNIFORM_PDF
+#undef RAN_SPHERE_CAP_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
