commit 5cf7afb9eb229070d564830d6f9e32553b996d0f
parent c2d9d411adc7e599b1f9264730b0d28240bc8d01
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 29 Mar 2016 15:39:35 +0200
Improve code readability and comments
Diffstat:
| M | src/s4vs.c | | | 109 | +++++++++++++++++++++++++++++++++++++++++++++++++------------------------------ |
1 file changed, 68 insertions(+), 41 deletions(-)
diff --git a/src/s4vs.c b/src/s4vs.c
@@ -64,63 +64,81 @@ realization(struct ssp_rng* rng, struct context* ctx)
/* sample a starting point on the scene surface to start the random walk from */
for(;;) {
- /* Sample the Scene */
+ /* sample the scene */
float u, v, w;
- do { u = (float)ssp_rng_canonical(rng); } while(u >= 1.f);
- do { v = (float)ssp_rng_canonical(rng); } while(v >= 1.f);
- do { w = (float)ssp_rng_canonical(rng); } while(w >= 1.f);
+ do { u = (float)ssp_rng_canonical(rng); } while(u >= 1);
+ do { v = (float)ssp_rng_canonical(rng); } while(v >= 1);
+ do { w = (float)ssp_rng_canonical(rng); } while(w >= 1);
+ /* get a location: primitive ID and parametric coordinates */
S3D(scene_sample(ctx->scene, u, v, w, &prim, st));
- /* Retrieve the sampled geometric normal and position */
+ /* retrieve the sampled geometric normal and position */
S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
f3_normalize(normal, attrib.value);
S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
f3_set(origin, attrib.value);
- /* Cosine weighted sampling of the hemisphere around the sampled normal */
+ /* cosine weighted sampling of the hemisphere around the sampled normal */
ssp_ran_hemisphere_cos(rng, normal, direction);
- /* Trace ray */
+ /* trace a ray from the sampled location */
range[0] = 1.e-6f; range[1] = FLT_MAX;
- do {
+ for(;;) {
+ /* find first intersection with geometry to determine its distance */
S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
- range[0] += 1.e-6f;
- } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim)); /* Handle auto-intersection */
+ /* handle auto-intersection: hitting the same primitive means auto-intersection */
+ if(S3D_PRIMITIVE_EQ(&hit.prim, &prim)) {
+ /* just increase the minimum distance and retry */
+ range[0] += 1.e-6f;
+ }
+ else break; /* standard case: OK */
+ }
f3_normalize(hit.normal, hit.normal);
- /* Handle precision issues */
+ /* handle precision issues (no intersection means geometry leakage: error) */
if(S3D_HIT_NONE(&hit)) {
if(++nfailures >= MAX_FAILURES) {
- /* Handle infinite loop due to unexpected geometry */
+ /* too many errors: stop */
ctx->exit_failure = 1;
break;
}
}
- else break;
+ else break; /* standard case: OK */
}
if(ctx->exit_failure)
return;
- /* Diffuse random walk */
+ /* here the starting point and a propagation direction have been sampled */
+ /* and we have determined the distance of the geometry in this direction */
+
+ /* start diffuse random walk */
for(;;) {
+ /* sample a free path according to ks */
lambda = ssp_ran_exp(rng, ctx->ks);
if(lambda >= hit.distance) {
+ /* lambda exceeds the geometry distance: the random walk ends on geometry */
sum += hit.distance;
break;
} else {
+ /* lambda doesn't exceed the geometry distance */
+ /* the random walk can continue after a diffusion */
sum += lambda;
+ /* new location after lambda */
f3_add(origin, origin, f3_mulf(direction, direction, (float)lambda));
do {
+ /* sample a new direction */
ssp_ran_sphere_uniform(rng, direction);
- range[0] = 0.f;
+ range[0] = 0;
+ /* find the first intersection with the geometry */
+ /* to refresh the geometry's distance */
S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
- /* Handle precision issues */
+ /* handle precision issues (no intersection means geometry leakage: error) */
if(S3D_HIT_NONE(&hit)) {
if(++nfailures >= MAX_FAILURES) {
- /* Handle infinite loop due to unexpected geometry */
+ /* too many errors: stop */
ctx->exit_failure = 1;
break;
}
@@ -128,6 +146,7 @@ realization(struct ssp_rng* rng, struct context* ctx)
} while(S3D_HIT_NONE(&hit));
}
}
+ /* accumulate results */
#pragma omp atomic update
ctx->sum += sum;
#pragma omp atomic update
@@ -143,41 +162,36 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
struct time begin, end, elapsed;
struct s3d_shape* shape;
struct context ctx;
- double sig, length, tmp;
- float S, V;
- res_T res = RES_OK;
+ double length, var, sig;
+ float S, V, expected;
struct ssp_rng* rng;
struct ssp_rng_proxy* proxy;
unsigned i;
size_t nb_t, me;
+ res_T res = RES_OK;
- ASSERT(max_steps && ks > 0);
+ ASSERT(max_steps > 0 && ks > 0);
S3D(scene_instantiate(scene, &shape));
+ S3D(scene_begin_session(scene, S3D_SAMPLE|S3D_TRACE));
- res = s3d_scene_begin_session(scene, S3D_SAMPLE|S3D_TRACE);
- if(res != RES_OK) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR,
- "Couldn't begin the scene trace/sample session. \n");
- goto error;
- }
-
+ /* compute the expected result using a mesh-based method */
S3D(scene_compute_area(scene, &S));
S3D(scene_compute_volume(scene, &V));
-
- if(eq_epsf(S, 0.f, 1.e-6f) || S < 0.f) {
+ if(eq_epsf(S, 0, 1.e-6f) || S < 0) {
logger_print(LOGGER_DEFAULT, LOG_ERROR,
"No surface to sample. Is the scene empty?\n");
goto error;
}
-
- if(eq_epsf(V, 0.f, 1.e-6f) || V < 0.f) {
+ if(eq_epsf(V, 0, 1.e-6f) || V < 0) {
logger_print(LOGGER_DEFAULT, LOG_ERROR,
"Invalid volume \"%.2f\". The scene might not match the prerequisites:\n"
"it must be closed and its normals must point *into* the volume.\n", V);
goto error;
}
+ expected = 4 * V / S;
+ /* initialize context for MC computation */
ctx.scene = scene;
ctx.nfailures = 0;
ctx.exit_failure = 0;
@@ -189,29 +203,29 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
nb_t = 1;
me = 0;
+ time_current(&begin);
+
#pragma omp parallel private(me, rng)
{
+ /* create a RNG proxy and a derived RNG by thread */
#ifdef _OPENMP
me = (size_t)omp_get_thread_num();
#endif
#pragma omp single
{
nb_t = (size_t)omp_get_num_threads();
- CHECK(ssp_rng_proxy_create(NULL, &ssp_rng_threefry, nb_t, &proxy), RES_OK);
+ SSP(rng_proxy_create(NULL, &ssp_rng_threefry, nb_t, &proxy));
}
#pragma omp barrier
- CHECK(ssp_rng_proxy_create_rng(proxy, me, &rng), RES_OK);
+ SSP(rng_proxy_create_rng(proxy, me, &rng));
- time_current(&begin);
#pragma omp for
for (i = 0; i < max_steps; i++) {
realization(rng, &ctx);
}
}
- length = ctx.sum / (double)max_steps;
- tmp = ctx.sum2 / (double)max_steps - length * length;
- sig = tmp > 0 ? sqrt(tmp / (double)max_steps) : 0;
+
time_current(&end);
time_sub(&elapsed, &end, &begin);
time_dump(&elapsed, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
@@ -225,9 +239,15 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
goto error;
}
+ /* compute mean length and standard error from MC results */
+ length = ctx.sum / (double)max_steps;
+ var = ctx.sum2 / (double)max_steps - length * length;
+ /* numerical errors can trick us: if variance is negative pretend standard error is 0 */
+ sig = var > 0 ? sqrt(var / (double)max_steps) : 0;
+
logger_print(LOGGER_DEFAULT, LOG_OUTPUT,
"\n4V/S = %g ~ %g +/- %g\n# failures: %d/%lu\nElapsed time: %s\n",
- 4.0*V/S, length, sig, ctx.nfailures, max_steps, buf);
+ expected, length, sig, ctx.nfailures, max_steps, buf);
#ifdef _OPENMP
logger_print(LOGGER_DEFAULT, LOG_OUTPUT,
@@ -235,6 +255,7 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
#endif
exit:
+ /* clean data */
if(scene) {
int mask;
S3D(scene_get_session_mask(scene, &mask));
@@ -244,9 +265,11 @@ exit:
return res;
error:
+ res = RES_UNKNOWN_ERR;
goto exit;
}
+/* import shapes in a scene */
static res_T
import_obj(const char *filename, struct s3d_scene **scene) {
struct s3d_device* s3d = NULL;
@@ -288,6 +311,7 @@ int main(int argc, char *argv[]) {
double ks = 0;
res_T res = RES_OK;
+ /* check command arguments */
if(argc < 2 || argc > 4) {
if(argc < 2)
logger_print(LOGGER_DEFAULT, LOG_ERROR, "Argument required\n");
@@ -298,12 +322,14 @@ int main(int argc, char *argv[]) {
return RES_BAD_ARG;
}
+ /* import file's content in the scene */
res = import_obj(argv[1], &scene);
if(res != RES_OK) {
logger_print(LOGGER_DEFAULT, LOG_ERROR, "Invalid file '%s'\n", argv[1]);
return res;
}
+ /* set number of realizations */
if(argc < 3) {
nsteps = 10000; /* default value */
}
@@ -316,12 +342,13 @@ int main(int argc, char *argv[]) {
nsteps = (size_t)ns;
}
+ /* set ks */
if(argc < 4) {
- ks = 1.0;
+ ks = 1; /* default value */
}
else {
ks = atof(argv[3]);
- if(ks <= 0.0) {
+ if(ks <= 0) {
logger_print(LOGGER_DEFAULT, LOG_ERROR,
"Invalid k-scattering value `%s'\n", argv[3]);
return RES_BAD_ARG;
