star-sp

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

commit d011c481c4d08e9a7324cd75f4714344fef599ef
parent 69921a28f54332b506b651a1c39f8fec5dcdf763
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 12 Apr 2016 18:31:32 +0200

Added Henyey-Greenstein distribution - not tested yet.

Used Richard's code.

Diffstat:
Msrc/ssp.h | 77+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 77 insertions(+), 0 deletions(-)

diff --git a/src/ssp.h b/src/ssp.h @@ -500,6 +500,83 @@ ssp_ran_hemisphere_cos_pdf(const float up[3], const float sample[3]) return dot < 0.f ? 0.f : (float)(dot/PI); } +/******************************************************************************* + * Henyey Greenstein distribution + ******************************************************************************/ + +/* 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]) +{ + double epsilon_g, epsilon_mu; + ASSERT(sample && f3_is_normalized(sample)); + if (g>0) { + epsilon_g = 1-g; + epsilon_mu = 1-sample[2]; + } else { + epsilon_g = 1+g; + epsilon_mu = 1+sample[2]; + } + return (float)(1/(4*PI) + *epsilon_g*(2-epsilon_g) + *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]) +{ + double phi, R, cos_theta, sin_theta; + ASSERT(rng && sample); + phi = 2 * PI * ssp_rng_canonical(rng); + R = ssp_rng_canonical(rng); + 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[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]) +{ + double epsilon_g, epsilon_mu; + ASSERT(up && sample && f3_is_normalized(up) && f3_is_normalized(sample)); + if (g>0) { + epsilon_g = 1-g; + epsilon_mu = 1-f3_dot(sample,up); + } else { + epsilon_g = 1+g; + epsilon_mu = 1+f3_dot(sample,up); + } + return (float)(1/(4*PI) + *epsilon_g*(2-epsilon_g) + *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* +ssp_ran_sphere_hg + (struct ssp_rng* rng, const float up[3], const double g, float sample[4]) +{ + float sample_local[4]; + float basis[9]; + ASSERT(rng && up && sample && f3_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); +} + END_DECLS #endif /* SSP_H */