star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

commit bbec97c367a8ffb4e3c6ced66e6042ecaf381c04
parent 61de7b92c788bd5549c8be4655327df5ba42d134
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 12 Jul 2021 11:31:54 +0200

Merge branch 'release_0.11'

Diffstat:
MREADME.md | 14++++++++++----
Mcmake/CMakeLists.txt | 25+++++++++++++------------
Msrc/ssp.h | 79+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
Msrc/ssp_ran.c | 70++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/test_ssp_ran_spherical_zone.c | 43+++++++++++++++++++++++++++++++++++++++++++
Asrc/test_ssp_ran_spherical_zone.h | 160+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 373 insertions(+), 18 deletions(-)

diff --git a/README.md b/README.md @@ -1,8 +1,8 @@ # Star SamPling -The purpose of this library is to generate [pseudo] random numbers and random -variates. While it is partly based on C++11 random generators, its API remains -pure C. +The purpose of this library is to generate [pseudo] random numbers and +distributions. While it is partly based on C++11 random generators, its API +remains pure C. ## How to build @@ -34,9 +34,15 @@ variable to the directory that contains the `boost` include directory. ## Release notes +### Version 0.11 + +- Add the `ssp_ran_spherical_zone_uniform` distribution that uniformly samples + a position on a truncated spherical cap. +- Rename the library to `star-sp` to avoid conflicts with GCC's `libssp`. + ### Version 0.10 -- Add the `ssp_ran_exp_truncated` random variates. +- Add the `ssp_ran_exp_truncated` distribution. ### Version 0.9 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 3.0) +cmake_minimum_required(VERSION 3.1) project(star-sp C CXX) enable_testing() @@ -74,7 +74,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys ${Boost_LIBRARY_DIRS}) # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 10) +set(VERSION_MINOR 11) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) @@ -107,27 +107,27 @@ if(CMAKE_COMPILER_IS_GNUCXX) 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 +add_library(star-sp SHARED ${SSP_FILES_SRC} ${SSP_FILES_INC} ${SSP_FILES_INC_API}) +set_target_properties(star-sp PROPERTIES DEFINE_SYMBOL SSP_SHARED_BUILD LINKER_LANGUAGE CXX VERSION ${VERSION} SOVERSION ${VERSION_MAJOR}) -target_link_libraries(ssp RSys) +target_link_libraries(star-sp RSys) if(MSVC) - target_link_libraries(ssp ${Boost_LIBRARIES}) + target_link_libraries(star-sp ${Boost_LIBRARIES}) # disable autolink - set_target_properties(ssp PROPERTIES COMPILE_FLAGS "/DBOOST_ALL_NO_LIB") + set_target_properties(star-sp PROPERTIES COMPILE_FLAGS "/DBOOST_ALL_NO_LIB") elseif(CMAKE_COMPILER_IS_GNUCXX) if(BUILD_R123_AES) - set_target_properties(ssp PROPERTIES COMPILE_FLAGS "-std=c++11 -maes") + set_target_properties(star-sp PROPERTIES COMPILE_FLAGS "-std=c++11 -maes") else() - set_target_properties(ssp PROPERTIES COMPILE_FLAGS "-std=c++11") + set_target_properties(star-sp PROPERTIES COMPILE_FLAGS "-std=c++11") endif() endif() -rcmake_setup_devel(ssp StarSP ${VERSION} star/ssp_version.h) +rcmake_setup_devel(star-sp StarSP ${VERSION} star/ssp_version.h) ################################################################################ # Add tests @@ -139,7 +139,7 @@ if(NOT NO_TEST) ${SSP_SOURCE_DIR}/${_name}.h ${SSP_SOURCE_DIR}/test_ssp_utils.h) - target_link_libraries(${_name} ssp RSys ${MATH_LIB}) + target_link_libraries(${_name} star-sp RSys ${MATH_LIB}) set(_libraries ${ARGN}) foreach(_lib ${_libraries}) target_link_libraries(${_name} ${_lib}) @@ -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) @@ -188,7 +189,7 @@ endif() ################################################################################ # Define output & install directories ################################################################################ -install(TARGETS ssp +install(TARGETS star-sp ARCHIVE DESTINATION bin LIBRARY DESTINATION lib RUNTIME DESTINATION bin) diff --git a/src/ssp.h b/src/ssp.h @@ -1048,8 +1048,8 @@ ssp_ran_uniform_disk_float } /******************************************************************************* - * Uniform distribution of a point ont a sphere cap - * pdf = 1/(2*PI*(1-cos(theta_max)) + * Uniform distribution of a point onto a sphere cap + * 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 onto 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