commit 9eb9c4280f89c2d96080229b8082a816ac3faf7e
parent 6f6cfadd6fee0b9a0c0fb0ad3cc4543e73d1e84f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 3 Oct 2016 16:36:24 +0200
Major updates of the random variates API
All random variates use doubles. The API of the spherical,
hemispherical, triangle and Henyey & Greenstein random variates are
consequently updated with respect to this structural choice.
Diffstat:
6 files changed, 242 insertions(+), 235 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -44,7 +44,7 @@ set(Random123_DIR ${_current_source_dir}/)
find_package(Random123 REQUIRED)
find_package(RCMake 0.2.3 REQUIRED)
-find_package(RSys 0.3 REQUIRED)
+find_package(RSys 0.4 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
include(rcmake_runtime)
diff --git a/src/ssp.h b/src/ssp.h
@@ -32,7 +32,7 @@
#ifndef SSP_H
#define SSP_H
-#include <rsys/float33.h>
+#include <rsys/double33.h>
#include <rsys/math.h>
#include <rsys/rsys.h>
@@ -253,43 +253,6 @@ ssp_rng_proxy_get_type
struct ssp_rng_type* type);
/*******************************************************************************
- * General discrete distribution with state
- ******************************************************************************/
-/* Create a discrete random variate data structure from a list of weights.
- * `weights' contain the weights of `nweights' discrete events. Its elements
- * must be positive but they needn't add up to one. */
-SSP_API res_T
-ssp_ranst_discrete_create
- (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
- struct ssp_ranst_discrete** ran);
-
-SSP_API res_T
-ssp_ranst_discrete_setup
- (struct ssp_ranst_discrete* discrete,
- const double* weights,
- const size_t nweights);
-
-SSP_API res_T
-ssp_ranst_discrete_ref_get
- (struct ssp_ranst_discrete* ran);
-
-SSP_API res_T
-ssp_ranst_discrete_ref_put
- (struct ssp_ranst_discrete* ran);
-
-/* Return the index of the sampled discret event. */
-SSP_API size_t
-ssp_ranst_discrete_get
- (struct ssp_rng* rng,
- const struct ssp_ranst_discrete* ran);
-
-/* Return the PDF of the discret event `i'. */
-SSP_API double
-ssp_ranst_discrete_pdf
- (const size_t i,
- const struct ssp_ranst_discrete* ran);
-
-/*******************************************************************************
* Miscellaneous distributions
******************************************************************************/
/* Random variate from the exponential distribution with mean `mu':
@@ -337,8 +300,8 @@ ssp_ran_lognormal_pdf
******************************************************************************/
/* Uniform sampling of an unit sphere. The PDF of the generated sample is
* stored in sample[3] */
-static INLINE float*
-ssp_ran_sphere_uniform(struct ssp_rng* rng, float sample[4])
+static INLINE double*
+ssp_ran_sphere_uniform(struct ssp_rng* rng, double sample[4])
{
double phi, cos_theta, sin_theta, v;
ASSERT(rng && sample);
@@ -346,57 +309,61 @@ ssp_ran_sphere_uniform(struct ssp_rng* rng, float sample[4])
v = ssp_rng_canonical(rng);
cos_theta = 1.0 - 2.0*v;
sin_theta = 2.0 * sqrt(v * (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)(1.0/(4.0*PI));
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ sample[3] = 1.0/(4.0*PI);
return sample;
}
/* Return the probability distribution for a sphere uniform random variate */
-static INLINE float
+static INLINE double
ssp_ran_sphere_uniform_pdf(void)
{
- return (float)(1.0/(4.0*PI));
+ return 1.0/(4.0*PI);
}
/*******************************************************************************
* Triangle distribution
******************************************************************************/
/* Uniform sampling of a triangle */
-static INLINE float*
+static INLINE double*
ssp_ran_triangle_uniform
(struct ssp_rng* rng,
- const float v0[3], const float v1[3], const float v2[3],/*triangle vertices*/
- float sample[3]) /* Sampled position */
+ 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 */
+ double sample[3]) /* Sampled position */
{
double sqrt_u, v, one_minus_u;
- int i;
+ double vec0[3];
+ double vec1[3];
ASSERT(rng && v0 && v1 && v2 && sample);
sqrt_u = sqrt(ssp_rng_canonical(rng));
v = ssp_rng_canonical(rng);
one_minus_u = 1.0 - sqrt_u;
- FOR_EACH(i, 0, 3) {
- sample[i] = (float)
- (v2[i] + one_minus_u * (v0[i] - v2[i]) + (v*sqrt_u) * (v1[i] - v2[i]));
- }
+ d3_sub(vec0, v0, v2);
+ d3_sub(vec1, v1, v2);
+ sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0];
+ sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1];
+ sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2];
return sample;
}
-static INLINE float
+static INLINE double
ssp_ran_triangle_uniform_pdf
- (const float v0[3],
- const float v1[3],
- const float v2[3])
+ (const double v0[3],
+ const double v1[3],
+ const double v2[3])
{
- float vec0[3], vec1[3], tmp[3];
+ double vec0[3], vec1[3], tmp[3];
ASSERT(v0 && v1 && v2);
- f3_sub(vec0, v0, v2);
- f3_sub(vec1, v1, v2);
- return 1.f / (f3_len(f3_cross(tmp, vec0, vec1)) * 0.5f);
+ d3_sub(vec0, v0, v2);
+ d3_sub(vec1, v1, v2);
+ return 1.0 / (d3_len(d3_cross(tmp, vec0, vec1)) * 0.5);
}
/*******************************************************************************
@@ -404,58 +371,58 @@ ssp_ran_triangle_uniform_pdf
******************************************************************************/
/* 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])
+static INLINE double*
+ssp_ran_hemisphere_uniform_local(struct ssp_rng* rng, double 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));
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ sample[3] = 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])
+static INLINE double
+ssp_ran_hemisphere_uniform_local_pdf(const double sample[3])
{
ASSERT(sample);
- return sample[2] < 0.f ? 0.f : (float)(1.0/(2.0*PI));
+ return sample[2] < 0.f ? 0.f : 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*
+static INLINE double*
ssp_ran_hemisphere_uniform
- (struct ssp_rng* rng, const float up[3], float sample[4])
+ (struct ssp_rng* rng, const double up[3], double sample[4])
{
- float basis[9];
- float sample_local[4];
- ASSERT(rng && up && sample && f3_is_normalized(up));
+ double basis[9];
+ double sample_local[4];
+ ASSERT(rng && up && sample && d3_is_normalized(up));
ssp_ran_hemisphere_uniform_local(rng, sample_local);
sample[3] = sample_local[3];
- return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
}
/* 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])
+static INLINE double
+ssp_ran_hemisphere_uniform_pdf(const double up[3], const double sample[3])
{
- ASSERT(up && sample && f3_is_normalized(up));
- return f3_dot(sample, up) < 0.f ? 0.f : (float)(1.0/(2.0*PI));
+ ASSERT(up && sample && d3_is_normalized(up));
+ return d3_dot(sample, up) < 0.f ? 0.f : 1.0/(2.0*PI);
}
/* Cosine weighted 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_cos_local(struct ssp_rng* rng, float sample[4])
+static INLINE double*
+ssp_ran_hemisphere_cos_local(struct ssp_rng* rng, double sample[4])
{
double phi, cos_theta, sin_theta, v;
ASSERT(rng && sample);
@@ -463,44 +430,44 @@ ssp_ran_hemisphere_cos_local(struct ssp_rng* rng, float sample[4])
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);
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ sample[3] = cos_theta / PI;
return sample;
}
/* Return the probability distribution for an hemispheric cosine weighted
* random variate with an implicit up direction in Z */
-static INLINE float
-ssp_ran_hemisphere_cos_local_pdf(const float sample[3])
+static INLINE double
+ssp_ran_hemisphere_cos_local_pdf(const double sample[3])
{
ASSERT(sample);
- return sample[2] < 0.f ? 0.f : (float)(sample[2] / PI);
+ return sample[2] < 0.f ? 0.f : 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])
+static INLINE double*
+ssp_ran_hemisphere_cos(struct ssp_rng* rng, const double up[3], double sample[4])
{
- float sample_local[4];
- float basis[9];
- ASSERT(rng && up && sample && f3_is_normalized(up));
+ double sample_local[4];
+ double basis[9];
+ ASSERT(rng && up && sample && d3_is_normalized(up));
ssp_ran_hemisphere_cos_local(rng, sample_local);
sample[3] = sample_local[3];
- return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
}
/* Return the probability distribution for an hemispheric cosine weighted
* random variate with an explicit `up' direction */
-static INLINE float
-ssp_ran_hemisphere_cos_pdf(const float up[3], const float sample[3])
+static INLINE double
+ssp_ran_hemisphere_cos_pdf(const double up[3], const double sample[3])
{
- float dot;
+ double dot;
ASSERT(up && sample);
- dot = f3_dot(sample, up);
- return dot < 0.f ? 0.f : (float)(dot/PI);
+ dot = d3_dot(sample, up);
+ return dot < 0.f ? 0.f : dot/PI;
}
/*******************************************************************************
@@ -508,11 +475,11 @@ ssp_ran_hemisphere_cos_pdf(const float up[3], const float sample[3])
******************************************************************************/
/* Return the probability distribution for a Henyey-Greenstein random
* variate with an implicit incident direction in Z */
-static INLINE float
-ssp_ran_sphere_hg_local_pdf(const double g, const float sample[3])
+static INLINE double
+ssp_ran_sphere_hg_local_pdf(const double g, const double sample[3])
{
double epsilon_g, epsilon_mu;
- ASSERT(sample && f3_is_normalized(sample));
+ ASSERT(sample && d3_is_normalized(sample));
if (g>0) {
epsilon_g = 1-g;
epsilon_mu = 1-sample[2];
@@ -520,16 +487,16 @@ ssp_ran_sphere_hg_local_pdf(const double g, const float sample[3])
epsilon_g = 1+g;
epsilon_mu = 1+sample[2];
}
- return (float)(1/(4*PI)
+ return 1/(4*PI)
*epsilon_g*(2-epsilon_g)
- *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5));
+ *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5);
}
/* Henyey-Greenstein sampling of an unit sphere for an incident direction
* that is implicitly the Z axis.
* The PDF of the generated sample is stored in sample[3] */
-static INLINE float*
-ssp_ran_sphere_hg_local(struct ssp_rng* rng, const double g, float sample[4])
+static INLINE double*
+ssp_ran_sphere_hg_local(struct ssp_rng* rng, const double g, double sample[4])
{
double phi, R, cos_theta, sin_theta;
ASSERT(rng && sample);
@@ -538,48 +505,88 @@ ssp_ran_sphere_hg_local(struct ssp_rng* rng, const double g, float sample[4])
cos_theta = 2*R*(1+g)*(1+g)*(g*(R-1)+1)
/((1-g*(1-2*R))*(1-g*(1-2*R)))-1;
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[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
sample[3] = ssp_ran_sphere_hg_local_pdf(g,sample);
return sample;
}
/* Return the probability distribution for a Henyey-Greenstein random
* variate with an explicit incident direction that is `up' */
-static INLINE float
-ssp_ran_sphere_hg_pdf(const float up[3], const double g, const float sample[3])
+static INLINE double
+ssp_ran_sphere_hg_pdf
+ (const double up[3],
+ const double g,
+ const double sample[3])
{
double epsilon_g, epsilon_mu;
- ASSERT(up && sample && f3_is_normalized(up) && f3_is_normalized(sample));
+ ASSERT(up && sample && d3_is_normalized(up) && d3_is_normalized(sample));
if (g>0) {
epsilon_g = 1-g;
- epsilon_mu = 1-f3_dot(sample,up);
+ epsilon_mu = 1-d3_dot(sample,up);
} else {
epsilon_g = 1+g;
- epsilon_mu = 1+f3_dot(sample,up);
+ epsilon_mu = 1+d3_dot(sample,up);
}
- return (float)(1/(4*PI)
+ return 1/(4*PI)
*epsilon_g*(2-epsilon_g)
- *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5));
+ *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5);
}
/* Henyey-Greenstein sampling of an unit sphere for an incident direction
* that is `up'.
* The PDF of the generated sample is stored in sample[3] */
-static INLINE float*
+static INLINE double*
ssp_ran_sphere_hg
- (struct ssp_rng* rng, const float up[3], const double g, float sample[4])
+ (struct ssp_rng* rng, const double up[3], const double g, double sample[4])
{
- float sample_local[4];
- float basis[9];
- ASSERT(rng && up && sample && f3_is_normalized(up));
+ double sample_local[4];
+ double basis[9];
+ ASSERT(rng && up && sample && d3_is_normalized(up));
ssp_ran_sphere_hg_local(rng, g, sample_local);
- sample[3]=sample_local[3];
- return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+ sample[3] = sample_local[3];
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
}
/*******************************************************************************
+ * General discrete distribution with state
+ ******************************************************************************/
+/* Create a discrete random variate data structure from a list of weights.
+ * `weights' contain the weights of `nweights' discrete events. Its elements
+ * must be positive but they needn't add up to one. */
+SSP_API res_T
+ssp_ranst_discrete_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ struct ssp_ranst_discrete** ran);
+
+SSP_API res_T
+ssp_ranst_discrete_setup
+ (struct ssp_ranst_discrete* discrete,
+ const double* weights,
+ const size_t nweights);
+
+SSP_API res_T
+ssp_ranst_discrete_ref_get
+ (struct ssp_ranst_discrete* ran);
+
+SSP_API res_T
+ssp_ranst_discrete_ref_put
+ (struct ssp_ranst_discrete* ran);
+
+/* Return the index of the sampled discret event. */
+SSP_API size_t
+ssp_ranst_discrete_get
+ (struct ssp_rng* rng,
+ const struct ssp_ranst_discrete* ran);
+
+/* Return the PDF of the discret event `i'. */
+SSP_API double
+ssp_ranst_discrete_pdf
+ (const size_t i,
+ const struct ssp_ranst_discrete* ran);
+
+/*******************************************************************************
* Gaussian distribution with state
******************************************************************************/
SSP_API res_T
diff --git a/src/test_ssp_ran_hemisphere.c b/src/test_ssp_ran_hemisphere.c
@@ -29,7 +29,7 @@
#include "ssp.h"
#include "test_ssp_utils.h"
-#include <rsys/float4.h>
+#include <rsys/double4.h>
#define NSAMPS 128
@@ -38,12 +38,12 @@ 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];
+ double samps0[NSAMPS][4];
+ double samps1[NSAMPS][4];
+ double samps2[NSAMPS][4];
+ double samps3[NSAMPS][4];
int i = 0, j = 0;
- float* f;
+ double* f;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -53,111 +53,111 @@ main(int argc, char** argv)
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];
+ double frame[9];
+ double up[3] = {0.f, 0.f, 1.f};
+ double xyz[3];
uint64_t 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);
+ CHECK(d3_is_normalized(f), 1);
+ CHECK(eq_eps(f[3], (float)(1.0/(2.0*PI)), 1.e-6), 1);
+ CHECK(eq_eps(f[3], ssp_ran_hemisphere_uniform_local_pdf(f), 1.e-6), 1);
+ CHECK(d3_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);
+ CHECK(d4_eq_eps(f, samps0[i], 1.e-6), 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);
+ up[0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ d3_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);
+ CHECK(d3_eq_eps(samps0[i], samps1[i], 1.e-6), 0);
+ CHECK(d3_is_normalized(f), 1);
+ CHECK(d3_dot(f, up) >= 0.f, 1);
+ CHECK(eq_eps(f[3], 1.0/(2.0*PI), 1.e-6 ), 1);
+ CHECK(eq_eps(f[3], ssp_ran_hemisphere_uniform_pdf(up, f), 1.e-6), 1);
+
+ d33_basis(frame, up);
+ d33_muld3(xyz, frame, samps0[i]);
+ CHECK(d3_eq_eps(samps1[i], xyz, 1.e-6), 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);
+ CHECK(d3_eq_eps(samps0[i], samps0[j], 1.e-6), 0);
+ CHECK(d3_eq_eps(samps1[i], samps1[j], 1.e-6), 0);
}
}
- samps1[1][0] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][1] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][2] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- f3_normalize(samps1[1], samps1[1]);
+ samps1[1][0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ samps1[1][1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ samps1[1][2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ d3_normalize(samps1[1], samps1[1]);
ssp_rng_set(rng0, 0);
ssp_ran_hemisphere_uniform(rng0, samps1[1], samps0[0]);
ssp_rng_set(rng0, 0);
ssp_ran_hemisphere_uniform(rng0, samps1[1], samps1[1]);
- CHECK(f4_eq(samps0[0], samps1[1]), 1);
+ CHECK(d4_eq(samps0[0], samps1[1]), 1);
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];
+ double frame[9];
+ double up[3] = { 0.0, 0.0, 1.0 };
+ double xyz[3];
uint64_t 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);
+ CHECK(d3_eq_eps(samps0[i], samps2[i], 1.e-6), 0);
+ CHECK(d3_is_normalized(f), 1);
+ CHECK(eq_eps(f[3], f[2]/PI, 1.e-6), 1);
+ CHECK(eq_eps(f[3], ssp_ran_hemisphere_cos_local_pdf(f), 1.e-6), 1);
+ CHECK(d3_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);
+ CHECK(d4_eq_eps(f, samps2[i], 1.e-6), 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);
+ up[0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ up[2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ d3_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);
+ CHECK(d3_eq_eps(samps2[i], samps3[i], 1.e-6), 0);
+ CHECK(d3_eq_eps(samps1[i], samps3[i], 1.e-6), 0);
+ CHECK(d3_is_normalized(f), 1);
+ CHECK(d3_dot(f, up) >= 0.f, 1);
+ CHECK(eq_eps(f[3], d3_dot(f, up)/PI, 1.e-6), 1);
+ CHECK(eq_eps(f[3], ssp_ran_hemisphere_cos_pdf(f, up), 1.e-6), 1);
+
+ d33_basis(frame, up);
+ d33_muld3(xyz, frame, samps2[i]);
+ CHECK(d3_eq_eps(samps3[i], xyz, 1.e-6), 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);
+ CHECK(d3_eq_eps(samps2[i], samps2[j], 1.e-6), 0);
+ CHECK(d3_eq_eps(samps3[i], samps3[j], 1.e-6), 0);
}
}
- samps1[1][0] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][1] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][2] = (float)ssp_rng_uniform_double(rng1, -1.0, 1.0);
- f3_normalize(samps1[1], samps1[1]);
+ samps1[1][0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ samps1[1][1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ samps1[1][2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
+ d3_normalize(samps1[1], samps1[1]);
ssp_rng_set(rng0, 0);
ssp_ran_hemisphere_cos(rng0, samps1[1], samps0[0]);
ssp_rng_set(rng0, 0);
ssp_ran_hemisphere_cos(rng0, samps1[1], samps1[1]);
- CHECK(f4_eq(samps0[0], samps1[1]), 1);
+ CHECK(d4_eq(samps0[0], samps1[1]), 1);
ssp_rng_ref_put(rng0);
ssp_rng_ref_put(rng1);
diff --git a/src/test_ssp_ran_hg.c b/src/test_ssp_ran_hg.c
@@ -49,17 +49,17 @@ main(int argc, char** argv)
double g = ssp_rng_uniform_double(rng, -1, +1);
double sum_cos = 0;
double sum_cos_local = 0;
- float dir[4], up[4] = {0, 0, 1};
+ double dir[4], up[4] = {0, 0, 1};
FOR_EACH(j, 0, NSAMPS) {
/* HG relative to the Z axis */
ssp_ran_sphere_hg_local(rng, g, dir);
- sum_cos_local += f3_dot(up, dir);
+ sum_cos_local += d3_dot(up, dir);
}
FOR_EACH(j, 0, NSAMPS) {
/* HG relative to a up uniformaly sampled */
ssp_ran_hemisphere_uniform_local(rng, up);
ssp_ran_sphere_hg(rng, up, g, dir);
- sum_cos += f3_dot(up, dir);
+ sum_cos += d3_dot(up, dir);
}
/* ...on average cos(up, dir) should be g */
CHECK(eq_eps(sum_cos_local / NSAMPS, g, 2.5e-2), 1);
diff --git a/src/test_ssp_ran_sphere.c b/src/test_ssp_ran_sphere.c
@@ -36,8 +36,8 @@ main(int argc, char** argv)
{
struct ssp_rng* rng;
struct mem_allocator allocator;
- float samps[NSAMPS][4];
- float* f = NULL;
+ double samps[NSAMPS][4];
+ double* f = NULL;
int i = 0, j = 0;
(void)argc, (void)argv;
@@ -48,11 +48,11 @@ main(int argc, char** argv)
FOR_EACH(i, 0, NSAMPS) {
f = ssp_ran_sphere_uniform(rng, samps[i]);
CHECK(f, samps[i]);
- CHECK(f3_is_normalized(f), 1);
- CHECK(eq_epsf(samps[i][3], (float)(1.0/(4.0*PI)), 1.e-6f), 1);
- CHECK(eq_epsf(samps[i][3], ssp_ran_sphere_uniform_pdf(), 1.e-6f), 1);
+ CHECK(d3_is_normalized(f), 1);
+ CHECK(eq_eps(samps[i][3], 1.0/(4.0*PI), 1.e-6), 1);
+ CHECK(eq_eps(samps[i][3], ssp_ran_sphere_uniform_pdf(), 1.e-6), 1);
FOR_EACH(j, 0, i) {
- CHECK(f3_eq_eps(samps[j], samps[i], 1.e-6f), 0);
+ CHECK(d3_eq_eps(samps[j], samps[i], 1.e-6), 0);
}
}
diff --git a/src/test_ssp_ran_triangle.c b/src/test_ssp_ran_triangle.c
@@ -38,10 +38,10 @@ main(int argc, char** argv)
{
struct ssp_rng* rng;
struct mem_allocator allocator;
- float samps[NSAMPS][3];
- float A[3], B[3], C[3];
- float v0[3], v1[3], v2[3], v3[3], v4[3], v5[3];
- float plane[4];
+ double samps[NSAMPS][3];
+ double A[3], B[3], C[3];
+ double v0[3], v1[3], v2[3], v3[3], v4[3], v5[3];
+ double plane[4];
size_t counter[2];
size_t nsteps;
size_t i;
@@ -52,50 +52,50 @@ main(int argc, char** argv)
CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
FOR_EACH(i, 0, 3) {
- A[i] = (float)ssp_rng_uniform_double(rng, 0.0, 1.0);
- B[i] = (float)ssp_rng_uniform_double(rng, 0.0, 1.0);
- C[i] = (float)ssp_rng_uniform_double(rng, 0.0, 1.0);
+ A[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
+ B[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
+ C[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
}
- f3_sub(v0, B, A);
- f3_sub(v1, C, A);
- f3_sub(v2, C, B);
- f3_minus(v3, v0);
- f3_minus(v4, v1);
- f3_minus(v5, v2);
- f3_cross(plane, v0, v1);
- plane[3] = -f3_dot(plane, C);
+ d3_sub(v0, B, A);
+ d3_sub(v1, C, A);
+ d3_sub(v2, C, B);
+ d3_minus(v3, v0);
+ d3_minus(v4, v1);
+ d3_minus(v5, v2);
+ d3_cross(plane, v0, v1);
+ plane[3] = -d3_dot(plane, C);
FOR_EACH(i, 0, NSAMPS) {
- float tmp0[3], tmp1[3];
- float dot = 0.f;
- float area = 0.f;
+ double tmp0[3], tmp1[3];
+ double dot = 0.0;
+ double area = 0.0;
CHECK(ssp_ran_triangle_uniform(rng, A, B, C, samps[i]), samps[i]);
- CHECK(eq_epsf(f3_dot(plane, samps[i]), -plane[3], 1.e-6f), 1);
-
- f3_sub(tmp0, samps[i], A);
- dot = f3_dot(f3_cross(tmp0, tmp0, v0), f3_cross(tmp1, v1, v0));
- CHECK(signf(dot), 1.f );
- f3_sub(tmp0, samps[i], B);
- dot = f3_dot(f3_cross(tmp0, tmp0, v2), f3_cross(tmp1, v3, v2));
- CHECK(signf(dot), 1.f);
- f3_sub(tmp0, samps[i], C);
- dot = f3_dot(f3_cross(tmp0, tmp0, v4), f3_cross(tmp1, v5, v4));
- CHECK(signf(dot), 1.f);
-
- area = f3_len(tmp1) * 0.5f;
- CHECK(eq_epsf(1.f / area, ssp_ran_triangle_uniform_pdf(A, B, C), 1.e-6f), 1);
+ CHECK(eq_eps(d3_dot(plane, samps[i]), -plane[3], 1.e-6), 1);
+
+ d3_sub(tmp0, samps[i], A);
+ dot = d3_dot(d3_cross(tmp0, tmp0, v0), d3_cross(tmp1, v1, v0));
+ CHECK(sign(dot), 1.0);
+ d3_sub(tmp0, samps[i], B);
+ dot = d3_dot(d3_cross(tmp0, tmp0, v2), d3_cross(tmp1, v3, v2));
+ CHECK(sign(dot), 1.0);
+ d3_sub(tmp0, samps[i], C);
+ dot = d3_dot(d3_cross(tmp0, tmp0, v4), d3_cross(tmp1, v5, v4));
+ CHECK(sign(dot), 1.0);
+
+ area = d3_len(tmp1) * 0.5f;
+ CHECK(eq_eps(1.0 / area, ssp_ran_triangle_uniform_pdf(A, B, C), 1.e-6), 1);
}
nsteps = 10000;
counter[0] = counter[1] = 0;
- f3(A, -1.f, 0.f, 0.f);
- f3(B, 1.f, 0.f, 0.f);
- f3(C, 0.f, 1.f, 0.f);
+ d3(A, -1.0, 0.0, 0.0);
+ d3(B, 1.0, 0.0, 0.0);
+ d3(C, 0.0, 1.0, 0.0);
FOR_EACH(i, 0, nsteps) {
ssp_ran_triangle_uniform(rng, A, B, C, samps[0]);
- if(samps[0][0] < 0.f)
+ if(samps[0][0] < 0.0)
counter[0] += 1;
else
counter[1] += 1;
