commit 543cf67f1416ff96ed4be215b044d9efda8c3913
parent 746f4b8ea0ebb9117c74b33cd577aa54518453dd
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 25 Oct 2017 17:50:32 +0200
Add world version of uniform_disk distribution.
Diffstat:
3 files changed, 95 insertions(+), 31 deletions(-)
diff --git a/src/ssp.h b/src/ssp.h
@@ -833,35 +833,63 @@ ssp_ranst_piecewise_linear_float_pdf
/*******************************************************************************
* Uniform disk distribution.
- *
- * TODO: provide a world space version.
******************************************************************************/
SSP_API double*
-ssp_ran_uniform_disk
+ssp_ran_uniform_disk_local
(struct ssp_rng* rng,
const double radius,
- double pt[2],
+ double pt[3], /* pt[2] <= 0 */
double* pdf); /* Can be NULL => pdf not computed */
SSP_API float*
-ssp_ran_uniform_disk_float
+ssp_ran_uniform_disk_float_local
(struct ssp_rng* rng,
const float radius,
- float pt[2],
+ float pt[3], /* pt[2] <= 0 */
float* pdf); /* Can be NULL => pdf not computed */
static INLINE double
-ssp_ran_uniform_disk_pdf(const double radius)
+ssp_ran_uniform_disk_local_pdf(const double radius)
{
return 1 / (radius * radius);
}
static INLINE float
-ssp_ran_uniform_disk_float_pdf(const float radius)
+ssp_ran_uniform_disk_float_local_pdf(const float radius)
{
return 1 / (radius * radius);
}
+static INLINE double*
+ssp_ran_uniform_disk
+ (struct ssp_rng* rng,
+ const double radius,
+ const double up[3],
+ double pt[3],
+ double* pdf) /* Can be NULL => pdf not computed */
+{
+ double sample_local[3];
+ double basis[9];
+ ASSERT(rng && up && pt && radius > 0);
+ ssp_ran_uniform_disk_local(rng, radius, sample_local, pdf);
+ return d33_muld3(pt, d33_basis(basis, up), sample_local);
+}
+
+static INLINE float*
+ssp_ran_uniform_disk_float
+ (struct ssp_rng* rng,
+ const float radius,
+ const float up[3],
+ float pt[3],
+ float* pdf) /* Can be NULL => pdf not computed */
+{
+ float sample_local[3];
+ float basis[9];
+ ASSERT(rng && up && pt && radius > 0);
+ ssp_ran_uniform_disk_float_local(rng, radius, sample_local, pdf);
+ return f33_mulf3(pt, f33_basis(basis, up), sample_local);
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -457,10 +457,10 @@ ssp_ran_sphere_hg_float_pdf
}
double*
-ssp_ran_uniform_disk
+ssp_ran_uniform_disk_local
(struct ssp_rng* rng,
const double radius,
- double pt[2],
+ double pt[3],
double* pdf)
{
double theta, r;
@@ -469,15 +469,16 @@ ssp_ran_uniform_disk
r = radius * sqrt(ssp_rng_canonical(rng));
pt[0] = r * cos(theta);
pt[1] = r * sin(theta);
+ pt[2] = 0;
if (pdf) *pdf = 1 / (radius * radius);
return pt;
}
float*
-ssp_ran_uniform_disk_float
+ssp_ran_uniform_disk_float_local
(struct ssp_rng* rng,
const float radius,
- float pt[2],
+ float pt[3],
float* pdf)
{
float theta, r;
@@ -486,6 +487,7 @@ ssp_ran_uniform_disk_float
r = radius * sqrtf(ssp_rng_canonical_float(rng));
pt[0] = r * cosf(theta);
pt[1] = r * sinf(theta);
+ pt[2] = 0;
if (pdf) *pdf = 1 / (radius * radius);
return pt;
}
\ No newline at end of file
diff --git a/src/test_ssp_ran_uniform_disk.h b/src/test_ssp_ran_uniform_disk.h
@@ -44,13 +44,25 @@
#if TYPE_FLOAT==0
#define REAL double
#define TEST test_double
+#define RNG_UNIFORM_R ssp_rng_uniform_double
+#define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_local
#define RNG_UNIFORM_DISK ssp_ran_uniform_disk
+#define R33_BASIS d33_basis
+#define R33_MULR3 d33_muld3
+#define R3_EQ_EPS d3_eq_eps
+#define R3_NORMALIZE d3_normalize
#define SQRT sqrt
#elif TYPE_FLOAT==1
#define REAL float
#define TEST test_float
+#define RNG_UNIFORM_R ssp_rng_uniform_float
+#define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_float_local
#define RNG_UNIFORM_DISK ssp_ran_uniform_disk_float
+#define R33_BASIS f33_basis
+#define R33_MULR3 f33_mulf3
+#define R3_EQ_EPS f3_eq_eps
+#define R3_NORMALIZE f3_normalize
#define SQRT sqrtf
#else
#error "TYPE_SUFFIX must be defined either 0 or 1"
@@ -59,25 +71,40 @@
static void
TEST()
{
- struct ssp_rng_proxy* proxy;
- struct ssp_rng* rng;
+ struct ssp_rng* rng0, *rng1;
struct mem_allocator allocator;
int i;
- unsigned r, c, nb = 0;
- REAL pt[2];
- REAL x = 0, x2 = 0;
+ int r, c, nb = 0;
+ REAL pt[3], pt2[3], up[3];
+ REAL frame[9];
+ REAL x_sum = 0, x2_sum = 0;
REAL mean, std;
- REAL exp_mean = 20 * 20 * NBS / ((REAL)PI * 100 * 100), exp_std = 0;
+ REAL exp_mean = NBS * (20 * 20 / ((REAL)PI * 100 * 100)), exp_std = 0;
unsigned counts[10][10];
+ uint64_t seed = 1234;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
- CHECK(ssp_rng_proxy_create(&allocator, &ssp_rng_threefry, 1, &proxy), RES_OK);
- CHECK(ssp_rng_proxy_create_rng(proxy, 0, &rng), RES_OK);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng0), RES_OK);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng1), RES_OK);
+
+ up[0] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[1] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(up, up);
+ R33_BASIS(frame, up);
+
+ ssp_rng_set(rng0, seed);
+ ssp_rng_set(rng1, seed);
memset(counts, 0, sizeof(counts));
for(i = 0; i < NBS; i++) {
- RNG_UNIFORM_DISK(rng, 100, pt, NULL);
+ REAL tmp[3];
+ RNG_UNIFORM_DISK_LOCAL(rng0, 100, pt, NULL);
+ RNG_UNIFORM_DISK(rng1, 100, up, pt2, NULL);
+ R33_MULR3(tmp, frame, pt);
+ CHECK(R3_EQ_EPS(tmp, pt2, (REAL)1.e-6), 1);
+ ASSERT(pt[2] == 0);
/* locate pt in a 10x10 grid */
r = (unsigned)((100 + pt[0]) / 20);
c = (unsigned)((100 + pt[1]) / 20);
@@ -85,26 +112,27 @@ TEST()
}
for(r = 0; r < 10; r++) {
- unsigned _x = (r >= 5 ? r - 4 : r - 5) * 20;
+ int x = (r >= 5 ? r - 4 : r - 5) * 20;
for(c = 0; c < 10; c++) {
- unsigned _y = (c >= 5 ? c - 4 : c - 5) * 20;
- unsigned r2 = _x * _x + _y * _y;
+ int y = (c >= 5 ? c - 4 : c - 5) * 20;
+ unsigned r2 = x * x + y * y;
if(r2 > 100 * 100)
/* this square is not (fully) in the disk */
continue;
++nb;
- x += counts[r][c];
- x2 += counts[r][c] * counts[r][c];
+ x_sum += counts[r][c];
+ x2_sum += counts[r][c] * counts[r][c];
}
}
- mean = x / nb;
- std = SQRT(x2 / nb - x / nb*x / nb);
+ mean = x_sum / nb;
+ std = x2_sum / nb - mean*mean;
+ std = std > 0 ? SQRT(std) : 0;
printf("%g %g\n", mean, std);
CHECK(fabs(mean - exp_mean) < 10, 1);
- CHECK(fabs(std - exp_std) < 100, 1);
+ CHECK(fabs(std - exp_std) < 200, 1);
- CHECK(ssp_rng_ref_put(rng), RES_OK);
- CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK);
+ CHECK(ssp_rng_ref_put(rng0), RES_OK);
+ CHECK(ssp_rng_ref_put(rng1), RES_OK);
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
@@ -113,6 +141,12 @@ TEST()
#undef REAL
#undef TEST
+#undef RNG_UNIFORM_R
+#undef RNG_UNIFORM_DISK_LOCAL
#undef RNG_UNIFORM_DISK
+#undef R33_BASIS
+#undef R33_MULR3
+#undef R3_EQ_EPS
+#undef R3_NORMALIZE
#undef SQRT
#undef TYPE_FLOAT
