star-sp

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

commit b92e881cc40a68ac16c5bf354488e0064fb24f28
parent 5b4bfe8c675f554652e9becc7a7c88a05262b772
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Wed, 25 Oct 2017 17:02:45 +0200

Change the way pdf is returned when sampling points/vectors.

Use an additional optional arg instead of using the last element of the
vector receiving the sampled value.

Diffstat:
Msrc/ssp.h | 141+++++++++++++++++++++++++++++++++++++++++--------------------------------------
Msrc/ssp_ran.c | 64+++++++++++++++++++++++++++++++++++++++++-----------------------
Msrc/test_ssp_ran_hemisphere.h | 24++++++++++++++----------
Msrc/test_ssp_ran_hg.h | 8++++----
Msrc/test_ssp_ran_sphere.h | 2+-
Msrc/test_ssp_ran_triangle.h | 7++++---
Msrc/test_ssp_ran_uniform_disk.h | 2+-
7 files changed, 139 insertions(+), 109 deletions(-)

diff --git a/src/ssp.h b/src/ssp.h @@ -354,24 +354,26 @@ ssp_ran_lognormal_float_pdf SSP_API double* ssp_ran_sphere_uniform (struct ssp_rng* rng, - double sample[4]); + double sample[3], + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_sphere_uniform_float (struct ssp_rng* rng, - float sample[4]); + float sample[3], + float* pdf); /* Can be NULL => pdf not computed */ /* Return the probability distribution for a sphere uniform random variate */ static INLINE double ssp_ran_sphere_uniform_pdf(void) { - return 1 / (4 * PI); + return 1 / (4*PI); } static INLINE float ssp_ran_sphere_uniform_float_pdf(void) { - return 1/(4*(float)PI); + return 1 / (4*(float)PI); } /******************************************************************************* @@ -384,7 +386,8 @@ ssp_ran_triangle_uniform 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 sample[3], /* Sampled position */ + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_triangle_uniform_float @@ -392,7 +395,8 @@ ssp_ran_triangle_uniform_float const float v0[3], /* Position of the 1st vertex */ const float v1[3], /* Position of the 2nd vertex */ const float v2[3], /* Position of the 3rd vertex */ - float sample[3]); /* Sampled position */ + float sample[3], /* Sampled position */ + float* pdf); /* Can be NULL => pdf not computed */ /* Return the probability distribution for a uniform point sampling on a triangle */ static INLINE double @@ -427,16 +431,18 @@ ssp_ran_triangle_uniform_float_pdf * 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] */ + * Z axis. */ SSP_API double* ssp_ran_hemisphere_uniform_local (struct ssp_rng* rng, - double sample[4]); + double sample[3], + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_hemisphere_uniform_float_local (struct ssp_rng* rng, - float sample[4]); + float sample[3], + float* pdf); /* Can be NULL => pdf not computed */ /* Return the probability distribution for an hemispheric uniform random * variate with an implicit up direction in Z */ @@ -451,22 +457,21 @@ static INLINE float ssp_ran_hemisphere_uniform_float_local_pdf(const float sample[3]) { ASSERT(sample); - return sample[2] < 0 ? 0 : 1 / (2 * (float)PI); + return sample[2] < 0 ? 0 : 1/(2*(float)PI); } -/* Uniform sampling of an unit hemisphere whose up direction is `up'. The PDF - * of the generated sample is stored in sample[3] */ +/* Uniform sampling of an unit hemisphere whose up direction is `up'. */ static INLINE double* ssp_ran_hemisphere_uniform (struct ssp_rng* rng, const double up[3], - double sample[4]) + double sample[3], + double* pdf) /* Can be NULL => pdf not computed */ { double basis[9]; - double sample_local[4]; + double sample_local[3]; ASSERT(rng && up && sample && d3_is_normalized(up)); - ssp_ran_hemisphere_uniform_local(rng, sample_local); - sample[3] = sample_local[3]; + ssp_ran_hemisphere_uniform_local(rng, sample_local, pdf); return d33_muld3(sample, d33_basis(basis, up), sample_local); } @@ -474,13 +479,13 @@ static INLINE float* ssp_ran_hemisphere_uniform_float (struct ssp_rng* rng, const float up[3], - float sample[4]) + float sample[3], + float* pdf) /* Can be NULL => pdf not computed */ { float basis[9]; - float sample_local[4]; + float sample_local[3]; ASSERT(rng && up && sample && f3_is_normalized(up)); - ssp_ran_hemisphere_uniform_float_local(rng, sample_local); - sample[3] = sample_local[3]; + ssp_ran_hemisphere_uniform_float_local(rng, sample_local, pdf); return f33_mulf3(sample, f33_basis(basis, up), sample_local); } @@ -501,7 +506,7 @@ ssp_ran_hemisphere_uniform_float_pdf const float sample[3]) { ASSERT(up && sample && f3_is_normalized(up)); - return f3_dot(sample, up) < 0 ? 0 : 1 / (2*(float)PI); + return f3_dot(sample, up) < 0 ? 0 : 1/(2*(float)PI); } /* Cosine weighted sampling of an unit hemisphere whose up direction is @@ -510,12 +515,14 @@ ssp_ran_hemisphere_uniform_float_pdf SSP_API double* ssp_ran_hemisphere_cos_local (struct ssp_rng* rng, - double sample[4]); + double sample[3], + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_hemisphere_cos_float_local (struct ssp_rng* rng, - float sample[4]); + float sample[3], + float* pdf); /* Can be NULL => pdf not computed */ /* Return the probability distribution for an hemispheric cosine weighted * random variate with an implicit up direction in Z */ @@ -523,29 +530,28 @@ static INLINE double ssp_ran_hemisphere_cos_local_pdf(const double sample[3]) { ASSERT(sample); - return sample[2] < 0 ? 0 : sample[2] / PI; + return sample[2] < 0 ? 0 : sample[2]/PI; } static INLINE float ssp_ran_hemisphere_cos_float_local_pdf(const float sample[3]) { ASSERT(sample); - return sample[2] < 0 ? 0 : sample[2] / (float)PI; + return sample[2] < 0 ? 0 : sample[2]/(float)PI; } -/* Cosine weighted sampling of an unit hemisphere whose up direction is `up'. - * The PDF of the generated sample is stored in sample[3] */ +/* Cosine weighted sampling of an unit hemisphere whose up direction is `up'. */ static INLINE double* ssp_ran_hemisphere_cos (struct ssp_rng* rng, const double up[3], - double sample[4]) + double sample[3], + double* pdf) /* Can be NULL => pdf not computed */ { - double sample_local[4]; + double sample_local[3]; double basis[9]; ASSERT(rng && up && sample && d3_is_normalized(up)); - ssp_ran_hemisphere_cos_local(rng, sample_local); - sample[3] = sample_local[3]; + ssp_ran_hemisphere_cos_local(rng, sample_local, pdf); return d33_muld3(sample, d33_basis(basis, up), sample_local); } @@ -553,13 +559,13 @@ static INLINE float* ssp_ran_hemisphere_cos_float (struct ssp_rng* rng, const float up[3], - float sample[4]) + float sample[3], + float* pdf) /* Can be NULL => pdf not computed */ { - float sample_local[4]; + float sample_local[3]; float basis[9]; ASSERT(rng && up && sample && f3_is_normalized(up)); - ssp_ran_hemisphere_cos_float_local(rng, sample_local); - sample[3] = sample_local[3]; + ssp_ran_hemisphere_cos_float_local(rng, sample_local, pdf); return f33_mulf3(sample, f33_basis(basis, up), sample_local); } @@ -584,26 +590,27 @@ ssp_ran_hemisphere_cos_float_pdf float dot; ASSERT(up && sample); dot = f3_dot(sample, up); - return dot < 0 ? 0 : dot / (float)PI; + return dot < 0 ? 0 : dot/(float)PI; } /******************************************************************************* * Henyey Greenstein distribution ******************************************************************************/ /* 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] */ + * that is implicitly the Z axis. */ SSP_API double* ssp_ran_sphere_hg_local (struct ssp_rng* rng, const double g, - double sample[4]); + double sample[3], + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_sphere_hg_float_local (struct ssp_rng* rng, const float g, - float sample[4]); + float sample[3], + float* pdf); /* Can be NULL => pdf not computed */ /* Return the probability distribution for a Henyey-Greenstein random * variate with an implicit incident direction in Z */ @@ -617,35 +624,20 @@ ssp_ran_sphere_hg_float_local_pdf (const float g, const float sample[3]); -/* Return the probability distribution for a Henyey-Greenstein random - * variate with an explicit incident direction that is `up' */ -SSP_API double -ssp_ran_sphere_hg_pdf - (const double up[3], - const double g, - const double sample[3]); - -SSP_API float -ssp_ran_sphere_hg_float_pdf - (const float up[3], - const float g, - const float sample[3]); - /* 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] */ + * that is `up'. */ static INLINE double* ssp_ran_sphere_hg (struct ssp_rng* rng, const double up[3], const double g, - double sample[4]) + double sample[3], + double* pdf) /* Can be NULL => pdf not computed */ { - double sample_local[4]; + double sample_local[3]; double basis[9]; ASSERT(-1 <= g && g <= +1 && rng && up && sample && d3_is_normalized(up)); - ssp_ran_sphere_hg_local(rng, g, sample_local); - sample[3] = sample_local[3]; + ssp_ran_sphere_hg_local(rng, g, sample_local, pdf); return d33_muld3(sample, d33_basis(basis, up), sample_local); } @@ -654,16 +646,30 @@ ssp_ran_sphere_hg_float (struct ssp_rng* rng, const float up[3], const float g, - float sample[4]) + float sample[3], + float* pdf) /* Can be NULL => pdf not computed */ { - float sample_local[4]; + float sample_local[3]; float basis[9]; ASSERT(-1 <= g && g <= +1 && rng && up && sample && f3_is_normalized(up)); - ssp_ran_sphere_hg_float_local(rng, g, sample_local); - sample[3] = sample_local[3]; + ssp_ran_sphere_hg_float_local(rng, g, sample_local, pdf); return f33_mulf3(sample, f33_basis(basis, up), sample_local); } +/* Return the probability distribution for a Henyey-Greenstein random +* variate with an explicit incident direction that is `up' */ +SSP_API double +ssp_ran_sphere_hg_pdf + (const double up[3], + const double g, + const double sample[3]); + +SSP_API float +ssp_ran_sphere_hg_float_pdf + (const float up[3], + const float g, + const float sample[3]); + /******************************************************************************* * General discrete distribution with state ******************************************************************************/ @@ -828,20 +834,21 @@ ssp_ranst_piecewise_linear_float_pdf /******************************************************************************* * Uniform disk distribution. * - * TODO: return the pdf in addition to the position. * TODO: provide a world space version. ******************************************************************************/ SSP_API double* ssp_ran_uniform_disk (struct ssp_rng* rng, const double radius, - double pt[2]); + double pt[2], + double* pdf); /* Can be NULL => pdf not computed */ SSP_API float* ssp_ran_uniform_disk_float (struct ssp_rng* rng, const float radius, - float pt[2]); + float pt[2], + float* pdf); /* Can be NULL => pdf not computed */ static INLINE double ssp_ran_uniform_disk_pdf(const double radius) diff --git a/src/ssp_ran.c b/src/ssp_ran.c @@ -169,7 +169,8 @@ ssp_ran_lognormal_float_pdf double* ssp_ran_sphere_uniform (struct ssp_rng* rng, - double sample[4]) + double sample[3], + double* pdf) { double phi, cos_theta, sin_theta, v; ASSERT(rng && sample); @@ -180,14 +181,15 @@ ssp_ran_sphere_uniform sample[0] = cos(phi) * sin_theta; sample[1] = sin(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = 1 / (4 * PI); + if(pdf) *pdf = 1 / (4*PI); return sample; } float* ssp_ran_sphere_uniform_float (struct ssp_rng* rng, - float sample[4]) + float sample[3], + float* pdf) { float phi, cos_theta, sin_theta, v; ASSERT(rng && sample); @@ -198,21 +200,23 @@ ssp_ran_sphere_uniform_float sample[0] = cosf(phi) * sin_theta; sample[1] = sinf(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = 1 / (4 * (float)PI); + if (pdf) *pdf = 1 / (4*(float)PI); return sample; } double* ssp_ran_triangle_uniform (struct ssp_rng* rng, - 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 */ + const double v0[3], + const double v1[3], + const double v2[3], + double sample[3], + double* pdf) { double sqrt_u, v, one_minus_u; double vec0[3]; double vec1[3]; + double tmp[3]; ASSERT(rng && v0 && v1 && v2 && sample); sqrt_u = sqrt(ssp_rng_canonical(rng)); @@ -224,6 +228,7 @@ ssp_ran_triangle_uniform 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]; + if (pdf) *pdf = 2 / d3_len(d3_cross(tmp, vec0, vec1)); return sample; } @@ -233,11 +238,13 @@ ssp_ran_triangle_uniform_float const float v0[3], const float v1[3], const float v2[3], - float sample[3]) + float sample[3], + float* pdf) { float sqrt_u, v, one_minus_u; float vec0[3]; float vec1[3]; + float tmp[3]; ASSERT(rng && v0 && v1 && v2 && sample); sqrt_u = sqrtf(ssp_rng_canonical_float(rng)); @@ -249,13 +256,15 @@ ssp_ran_triangle_uniform_float 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]; + if (pdf) *pdf = 2 / f3_len(f3_cross(tmp, vec0, vec1)); return sample; } double* ssp_ran_hemisphere_uniform_local (struct ssp_rng* rng, - double sample[4]) + double sample[3], + double* pdf) { double phi, cos_theta, sin_theta; ASSERT(rng && sample); @@ -265,14 +274,15 @@ ssp_ran_hemisphere_uniform_local sample[0] = cos(phi) * sin_theta; sample[1] = sin(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = 1 / (2 * PI); + if (pdf) *pdf = 1 / (2*PI); return sample; } float* ssp_ran_hemisphere_uniform_float_local (struct ssp_rng* rng, - float sample[4]) + float sample[3], + float* pdf) { float phi, cos_theta, sin_theta; ASSERT(rng && sample); @@ -282,14 +292,15 @@ ssp_ran_hemisphere_uniform_float_local sample[0] = cosf(phi) * sin_theta; sample[1] = sinf(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = 1 / (2 * (float)PI); + if (pdf) *pdf = 1 / (2*(float)PI); return sample; } double* ssp_ran_hemisphere_cos_local (struct ssp_rng* rng, - double sample[4]) + double sample[3], + double* pdf) { double phi, cos_theta, sin_theta, v; ASSERT(rng && sample); @@ -300,14 +311,15 @@ ssp_ran_hemisphere_cos_local sample[0] = cos(phi) * sin_theta; sample[1] = sin(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = cos_theta / PI; + if (pdf) *pdf = cos_theta / PI; return sample; } float* ssp_ran_hemisphere_cos_float_local (struct ssp_rng* rng, - float sample[4]) + float sample[3], + float* pdf) { float phi, cos_theta, sin_theta, v; ASSERT(rng && sample); @@ -318,7 +330,7 @@ ssp_ran_hemisphere_cos_float_local sample[0] = cosf(phi) * sin_theta; sample[1] = sinf(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = cos_theta / (float)PI; + if (pdf) *pdf = cos_theta / (float)PI; return sample; } @@ -364,7 +376,8 @@ double* ssp_ran_sphere_hg_local (struct ssp_rng* rng, const double g, - double sample[4]) + double sample[3], + double* pdf) { double phi, R, cos_theta, sin_theta; ASSERT(fabs(g) <= 1 && rng && sample); @@ -376,7 +389,7 @@ ssp_ran_sphere_hg_local 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); + if (pdf) *pdf = ssp_ran_sphere_hg_local_pdf(g, sample); return sample; } @@ -384,7 +397,8 @@ float* ssp_ran_sphere_hg_float_local (struct ssp_rng* rng, const float g, - float sample[4]) + float sample[3], + float* pdf) { float phi, R, cos_theta, sin_theta; ASSERT(fabsf(g) <= 1 && rng && sample); @@ -396,7 +410,7 @@ ssp_ran_sphere_hg_float_local sample[0] = cosf(phi) * sin_theta; sample[1] = sinf(phi) * sin_theta; sample[2] = cos_theta; - sample[3] = ssp_ran_sphere_hg_float_local_pdf(g, sample); + if (pdf) *pdf = ssp_ran_sphere_hg_float_local_pdf(g, sample); return sample; } @@ -446,7 +460,8 @@ double* ssp_ran_uniform_disk (struct ssp_rng* rng, const double radius, - double pt[2]) + double pt[2], + double* pdf) { double theta, r; ASSERT(rng && pt && radius > 0); @@ -454,6 +469,7 @@ ssp_ran_uniform_disk r = radius * sqrt(ssp_rng_canonical(rng)); pt[0] = r * cos(theta); pt[1] = r * sin(theta); + if (pdf) *pdf = 1 / (radius * radius); return pt; } @@ -461,7 +477,8 @@ float* ssp_ran_uniform_disk_float (struct ssp_rng* rng, const float radius, - float pt[2]) + float pt[2], + float* pdf) { float theta, r; ASSERT(rng && pt && radius > 0); @@ -469,5 +486,6 @@ ssp_ran_uniform_disk_float r = radius * sqrtf(ssp_rng_canonical_float(rng)); pt[0] = r * cosf(theta); pt[1] = r * sinf(theta); + if (pdf) *pdf = 1 / (radius * radius); return pt; } \ No newline at end of file diff --git a/src/test_ssp_ran_hemisphere.h b/src/test_ssp_ran_hemisphere.h @@ -104,6 +104,10 @@ TEST() CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng0), RES_OK); CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng1), RES_OK); + samps0[0][0] = 0; samps0[0][1] = 0; samps0[0][2] = 1; + f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps1[0], NULL); + f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps2[0], &samps2[0][3]); + ssp_rng_set(rng0, 0); FOR_EACH(i, 0, NSAMPS) { REAL frame[9]; @@ -112,7 +116,7 @@ TEST() uint64_t seed = ssp_rng_get(rng0); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_UNIFORM_LOCAL(rng1, samps0[i]); + f = RAN_HEMISPHERE_UNIFORM_LOCAL(rng1, samps0[i], &samps0[i][3]); CHECK(f, samps0[i]); CHECK(R3_IS_NORMALIZED(f), 1); CHECK(EQ_EPS_R(f[3], (1/(2*(REAL)PI)), (REAL)1.e-6), 1); @@ -120,7 +124,7 @@ TEST() CHECK(R3_DOT(f, up) >= 0, 1); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i]); + f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]); CHECK(f, samps1[i]); CHECK(R4_EQ_EPS(f, samps0[i], (REAL)1.e-6), 1); @@ -130,7 +134,7 @@ TEST() R3_NORMALIZE(up, up); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i]); + f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]); CHECK(R3_EQ_EPS(samps0[i], samps1[i], (REAL)1.e-6), 0); CHECK(R3_IS_NORMALIZED(f), 1); CHECK(R3_DOT(f, up) >= 0, 1); @@ -152,9 +156,9 @@ TEST() R3_NORMALIZE(samps1[1], samps1[1]); ssp_rng_set(rng0, 0); - RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps0[0]); + RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps0[0], &samps0[0][3]); ssp_rng_set(rng0, 0); - RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps1[1]); + RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps1[1], &samps1[1][3]); CHECK(R4_EQ(samps0[0], samps1[1]), 1); ssp_rng_set(rng0, 0); @@ -165,7 +169,7 @@ TEST() uint64_t seed = ssp_rng_get(rng0); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_COS_LOCAL(rng1, samps2[i]); + f = RAN_HEMISPHERE_COS_LOCAL(rng1, samps2[i], &samps2[i][3]); CHECK(f, samps2[i]); CHECK(R3_EQ_EPS(samps0[i], samps2[i], (REAL)1.e-6), 0); CHECK(R3_IS_NORMALIZED(f), 1); @@ -174,7 +178,7 @@ TEST() CHECK(R3_DOT(f, up) >= 0, 1); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_COS(rng1, up, samps3[i]); + f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]); CHECK(f, samps3[i]); CHECK(R4_EQ_EPS(f, samps2[i], (REAL)1.e-6), 1); @@ -184,7 +188,7 @@ TEST() R3_NORMALIZE(up, up); ssp_rng_set(rng1, seed); - f = RAN_HEMISPHERE_COS(rng1, up, samps3[i]); + f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]); CHECK(R3_EQ_EPS(samps2[i], samps3[i], (REAL)1.e-6), 0); CHECK(R3_EQ_EPS(samps1[i], samps3[i], (REAL)1.e-6), 0); CHECK(R3_IS_NORMALIZED(f), 1); @@ -207,9 +211,9 @@ TEST() R3_NORMALIZE(samps1[1], samps1[1]); ssp_rng_set(rng0, 0); - RAN_HEMISPHERE_COS(rng0, samps1[1], samps0[0]); + RAN_HEMISPHERE_COS(rng0, samps1[1], samps0[0], &samps0[0][3]); ssp_rng_set(rng0, 0); - RAN_HEMISPHERE_COS(rng0, samps1[1], samps1[1]); + RAN_HEMISPHERE_COS(rng0, samps1[1], samps1[1], &samps1[1][3]); CHECK(R4_EQ(samps0[0], samps1[1]), 1); ssp_rng_ref_put(rng0); diff --git a/src/test_ssp_ran_hg.h b/src/test_ssp_ran_hg.h @@ -76,16 +76,16 @@ TEST() REAL g = RNG_UNIFORM_R(rng, -1, +1); REAL sum_cos = 0; REAL sum_cos_local = 0; - REAL dir[4], up[4] = {0, 0, 1}; + REAL dir[3], pdf, up[4] = {0, 0, 1}; FOR_EACH(j, 0, NSAMPS) { /* HG relative to the Z axis */ - RAN_SPHERE_HG_LOCAL(rng, g, dir); + RAN_SPHERE_HG_LOCAL(rng, g, dir, &pdf); sum_cos_local += R3_DOT(up, dir); } FOR_EACH(j, 0, NSAMPS) { /* HG relative to a up uniformaly sampled */ - RAN_HEMISPHERE_UNIFORM_LOCAL(rng, up); - RAN_SPHERE_HG(rng, up, g, dir); + RAN_HEMISPHERE_UNIFORM_LOCAL(rng, up, &pdf); + RAN_SPHERE_HG(rng, up, g, dir, &pdf); sum_cos += R3_DOT(up, dir); } /* ...on average cos(up, dir) should be g */ diff --git a/src/test_ssp_ran_sphere.h b/src/test_ssp_ran_sphere.h @@ -72,7 +72,7 @@ TEST() CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK); FOR_EACH(i, 0, NSAMPS) { - f = RAN_SPHERE_UNIFORM(rng, samps[i]); + f = RAN_SPHERE_UNIFORM(rng, samps[i], &samps[i][3]); CHECK(f, samps[i]); CHECK(R3_IS_NORMALIZED(f), 1); CHECK(EQ_EPS(samps[i][3], 1/(4*(REAL)PI), (REAL)1.e-6), 1); diff --git a/src/test_ssp_ran_triangle.h b/src/test_ssp_ran_triangle.h @@ -78,6 +78,7 @@ TEST() REAL A[3], B[3], C[3]; REAL v0[3], v1[3], v2[3], v3[3], v4[3], v5[3]; REAL plane[4]; + REAL pdf; size_t counter[2]; size_t nsteps; size_t i; @@ -106,7 +107,7 @@ TEST() REAL dot = 0; REAL area = 0; - CHECK(RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[i]), samps[i]); + CHECK(RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[i], &pdf), samps[i]); CHECK(EQ_EPS_R(R3_DOT(plane, samps[i]), -plane[3], (REAL)1.e-6), 1); R3_SUB(tmp0, samps[i], A); @@ -129,7 +130,7 @@ TEST() R3(B, 1, 0, 0); R3(C, 0, 1, 0); FOR_EACH(i, 0, nsteps) { - RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0]); + RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0], NULL); if(samps[0][0] < 0.0) counter[0] += 1; else @@ -139,7 +140,7 @@ TEST() counter[0] = counter[1] = 0; FOR_EACH(i, 0, nsteps) { - RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0]); + RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0], NULL); if(samps[0][1] < 1 - 1/sqrt(2)) counter[0] += 1; else diff --git a/src/test_ssp_ran_uniform_disk.h b/src/test_ssp_ran_uniform_disk.h @@ -89,7 +89,7 @@ TEST() memset(counts, 0, sizeof(counts)); for(i = 0; i < NBS; i++) { - RNG_UNIFORM_DISK(rng, 100, pt); + RNG_UNIFORM_DISK(rng, 100, pt, NULL); /* locate pt in a 10x10 grid */ r = (unsigned)((100 + pt[0]) / 20); c = (unsigned)((100 + pt[1]) / 20);