star-gs

Literate program for a geometric sensitivity calculation
git clone git://git.meso-star.fr/star-gs.git
Log | Files | Refs | README | LICENSE

commit 78515cb13eac1d84caa5a91a7e7b1779731b2304
parent 92cd8557e50039ed77de7565d1ffbc4cbef6aec1
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 20 Jul 2022 17:41:59 +0200

Rewrite the source of sensitivity realization

Move scene description into separate functions

Diffstat:
Msrc/sgs_compute_sensitivity_translation.c | 198+++++++++++++++++++++++++++++++++++++++++++++++--------------------------------
1 file changed, 119 insertions(+), 79 deletions(-)

diff --git a/src/sgs_compute_sensitivity_translation.c b/src/sgs_compute_sensitivity_translation.c @@ -55,9 +55,99 @@ struct deformation_2d { double u[3]; }; +struct scene { + /* Upper and lower limit of the scene */ + double box_low[3]; + double box_upp[3]; + + double box_sz[3]; /* Size of the scene */ + + /* Geometry of the receiver on the box bottom plane */ + double recv_min[3]; + double recv_max[3]; + + double emit_e_threshold; /* Geometric upper bound of the radiative source */ + double emit_e_sz[3]; /* Size of the radiative source */ + double emit_s_sz[3]; /* Size of the specular surface */ +}; + /******************************************************************************* * Helper functions ******************************************************************************/ +static void +setup_scene(const struct sgs_geometry* geometry, struct scene* scn) +{ + ASSERT(geometry && scn); + + /* Compute the size of the geometry */ + sgs_geometry_get_aabb(geometry, scn->box_low, scn->box_upp); + d3_sub(scn->box_sz, scn->box_upp, scn->box_low); + + /* Setup the geometry of the receiver */ + scn->recv_min[X] = RECV_MIN[X]*scn->box_sz[X] + scn->box_low[X]; + scn->recv_min[Y] = RECV_MIN[Y]*scn->box_sz[Y] + scn->box_low[Y]; + scn->recv_min[Z] = 0; + scn->recv_max[X] = RECV_MAX[X]*scn->box_sz[X] + scn->box_low[X]; + scn->recv_max[Y] = RECV_MAX[Y]*scn->box_sz[Y] + scn->box_low[Y]; + scn->recv_max[Z] = 0; + + /* Setup the geometry of the radiative source */ + scn->emit_e_threshold = EMIT_E_THRESHOLD*scn->box_sz[Z] + scn->box_low[Z]; + + /* Compute the size of the radiative source/specular surface */ + scn->emit_e_sz[X] = EMIT_E_SZ[X]*scn->box_sz[X] + scn->box_low[X]; + scn->emit_e_sz[Y] = EMIT_E_SZ[Y]*scn->box_sz[Y] + scn->box_low[Y]; + scn->emit_e_sz[Z] = EMIT_E_SZ[Z]*scn->box_sz[Z] + scn->box_low[Z]; + scn->emit_s_sz[X] = EMIT_S_SZ[X]*scn->box_sz[X] + scn->box_low[X]; + scn->emit_s_sz[Y] = EMIT_S_SZ[Y]*scn->box_sz[Y] + scn->box_low[Y]; + scn->emit_s_sz[Z] = EMIT_S_SZ[Z]*scn->box_sz[Z] + scn->box_low[Z]; +} + +/* Return 1 if the `hit' lies on the receiver and 0 otherwise */ +static int +hit_receiver + (const struct scene* scn, + const double ray_org[3], + const double ray_dir[3], + const struct sgs_hit* hit, + double pos_recv[3]) +{ + ASSERT(scn && ray_org && ray_dir && hit && pos_recv); + + pos_recv[X] = ray_org[X] + hit->distance*ray_dir[X]; + pos_recv[Y] = ray_org[Y] + hit->distance*ray_dir[Y]; + pos_recv[Z] = ray_org[Z] + hit->distance*ray_dir[Z]; + + if(hit->surface != SGS_SURFACE_Z_NEG + || pos_recv[X] < scn->recv_min[X] || scn->recv_max[X] < pos_recv[X] + || pos_recv[Y] < scn->recv_min[Y] || scn->recv_max[Y] < pos_recv[Y]) { + return 0; /* The ray does not intersect the receiver */ + } else { + return 1; + } +} + +/* Return 1 if the `hit' lies on the source and 0 otherwise */ +static int +hit_source + (const struct scene* scn, + const double ray_org[3], + const double ray_dir[3], + const struct sgs_hit* hit, + double pos_emit_e[3]) +{ + pos_emit_e[X] = ray_org[X] + hit->distance*ray_dir[X]; + pos_emit_e[Y] = ray_org[Y] + hit->distance*ray_dir[Y]; + pos_emit_e[Z] = ray_org[Z] + hit->distance*ray_dir[Z]; + + if(hit->surface != SGS_SURFACE_X_POS + || pos_emit_e[Z] > scn->emit_e_threshold) { + return 0; /* The ray does not intersect the emitter */ + } else { + return 1; + } +} + static double compute_rho (const double position[2], /* Position sur la surface réfléchissante */ @@ -144,15 +234,9 @@ realisation struct deformation_2d proj_V_e; struct deformation_2d proj_V_s; - double recv_min[2]; - double recv_max[2]; - double emit_e_threshold; - double emit_e_sz[3], emit_e_sz_2d[2]; - double emit_s_sz[3], emit_s_sz_2d[2]; - - double box_low[3]; - double box_upp[3]; - double box_sz[3]; + struct scene scn; + double emit_e_sz_2d[2]; + double emit_s_sz_2d[2]; double normal_e[3]; double normal_s[3]; @@ -164,8 +248,6 @@ realisation double pos_emit_s[3], pos_emit_s_2d[2]; double pos_recv[3]; - double range[2]; - double rho; double grad_rho[3], grad_rho_2d[2]; double I; @@ -176,84 +258,42 @@ realisation double w = 0; double s = 0; - enum sgs_surface_type surf_emit_e; enum sgs_surface_type surf_emit_s; res_T res = RES_OK; (void)ithread; - range[0] = 0, range[1] = INF; - - /* Retrieve the scene data */ - sgs_geometry_get_aabb(sgs->geom, box_low, box_upp); - d3_sub(box_sz, box_upp, box_low); - recv_min[X] = RECV_MIN[X]*box_sz[X] + box_low[X]; - recv_min[Y] = RECV_MIN[Y]*box_sz[Y] + box_low[Y]; - recv_max[X] = RECV_MAX[X]*box_sz[X] + box_low[X]; - recv_max[Y] = RECV_MAX[Y]*box_sz[Y] + box_low[Y]; - emit_e_threshold = EMIT_E_THRESHOLD*box_sz[Z] + box_low[Z]; - emit_e_sz[X] = EMIT_E_SZ[X]*box_sz[X] + box_low[X]; - emit_e_sz[Y] = EMIT_E_SZ[Y]*box_sz[Y] + box_low[Y]; - emit_e_sz[Z] = EMIT_E_SZ[Z]*box_sz[Z] + box_low[Z]; - emit_s_sz[X] = EMIT_S_SZ[X]*box_sz[X] + box_low[X]; - emit_s_sz[Y] = EMIT_S_SZ[Y]*box_sz[Y] + box_low[Y]; - emit_s_sz[Z] = EMIT_S_SZ[Z]*box_sz[Z] + box_low[Z]; + setup_scene(sgs->geom, &scn); /* Sample the sensitivity emissive surface */ sgs_geometry_sample(sgs->geom, rng, &frag); - pos_emit_s[X] = frag.position[X]; - pos_emit_s[Y] = frag.position[Y]; - pos_emit_s[Z] = frag.position[Z]; - normal_s[X] = frag.normal[X]; - normal_s[Y] = frag.normal[Y]; - normal_s[Z] = frag.normal[Z]; + d3_set(pos_emit_s, frag.position); + d3_set(normal_s, frag.normal); surf_emit_s = frag.surface; /* Sample the cosine weighted sampling of the emissive direction */ ssp_ran_hemisphere_cos(rng, normal_s, dir_emit_s, NULL); - /* Trace the ray from the sensitivity emissive surface */ - sgs_geometry_trace_ray - (sgs->geom, pos_emit_s, dir_emit_s, range, surf_emit_s, &hit); - if(SGS_HIT_NONE(&hit)) { - res = RES_BAD_ARG; - goto error; - } - pos_recv[X] = pos_emit_s[X] + hit.distance*dir_emit_s[X]; - pos_recv[Y] = pos_emit_s[Y] + hit.distance*dir_emit_s[Y]; - pos_recv[Z] = pos_emit_s[Z] + hit.distance*dir_emit_s[Z]; - - /* The ray does not reach the receiver, the MC weight is NULL */ - if(hit.surface != SGS_SURFACE_Z_NEG - || pos_recv[X] < recv_min[X] || recv_max[X] < pos_recv[X] - || pos_recv[Y] < recv_min[Y] || recv_max[Y] < pos_recv[Y]) { - goto exit; - } + /* Helper macro used as syntactic sugar */ + #define TRACE_RAY(Org, Dir, StartFrom, Hit) { \ + const double range[2] = {0, DBL_MAX}; \ + sgs_geometry_trace_ray(sgs->geom, (Org), (Dir), range, (StartFrom), (Hit));\ + if(SGS_HIT_NONE(&hit)) { res = RES_BAD_ARG; goto error; } \ + } (void)0 - /* Find the reflected position */ + /* Trance the sampled ray and check that if reaches the receiver? */ + TRACE_RAY(pos_emit_s, dir_emit_s, surf_emit_s, &hit); + if(!hit_receiver(&scn, pos_emit_s, dir_emit_s, &hit, pos_recv)) goto exit; + + /* Trace the reflected ray and check that if reaches the source? */ reflect(dir_spec_s, dir_emit_s, normal_s); - sgs_geometry_trace_ray - (sgs->geom, pos_emit_s, dir_spec_s, range, surf_emit_s, &hit); - if(SGS_HIT_NONE(&hit)) { - res = RES_BAD_ARG; - goto error; - } - pos_emit_e[X] = pos_emit_s[X] + hit.distance*dir_spec_s[X]; - pos_emit_e[Y] = pos_emit_s[Y] + hit.distance*dir_spec_s[Y]; - pos_emit_e[Z] = pos_emit_s[Z] + hit.distance*dir_spec_s[Z]; - normal_e[X] = hit.normal[X]; - normal_e[Y] = hit.normal[Y]; - normal_e[Z] = hit.normal[Z]; - surf_emit_e = hit.surface; - - /* The reflected position is not an the emitter, the MC weight is null */ - if(surf_emit_e != SGS_SURFACE_X_POS || pos_emit_e[Z] > emit_e_threshold) { - goto exit; - } + TRACE_RAY(pos_emit_s, dir_spec_s, surf_emit_s, &hit); + if(!hit_source(&scn, pos_emit_s, dir_spec_s, &hit, pos_emit_e)) goto exit; + + #undef TRACE_RAY - dir_spec_e[X] = -dir_spec_s[X]; - dir_spec_e[Y] = -dir_spec_s[Y]; - dir_spec_e[Z] = -dir_spec_s[Z]; + d3_set(normal_e, hit.normal); + d3_minus(dir_spec_e, dir_spec_s); projection(V, normal_e, dir_spec_e, &proj_V_e); projection(V, normal_s, dir_emit_s, &proj_V_s); @@ -261,8 +301,8 @@ realisation pos_emit_s_2d[X] = pos_emit_s[X]; pos_emit_s_2d[Y] = pos_emit_s[Y]; - emit_s_sz_2d[X] = emit_s_sz[X]; - emit_s_sz_2d[Y] = emit_s_sz[Y]; + emit_s_sz_2d[X] = scn.emit_s_sz[X]; + emit_s_sz_2d[Y] = scn.emit_s_sz[Y]; rho = compute_rho(pos_emit_s_2d, emit_s_sz_2d); compute_grad_rho(pos_emit_s_2d, emit_s_sz_2d, grad_rho_2d); grad_rho[X] = grad_rho_2d[X]; @@ -271,8 +311,8 @@ realisation pos_emit_e_2d[X] = pos_emit_e[Y]; pos_emit_e_2d[Y] = pos_emit_e[Z]; - emit_e_sz_2d[X] = emit_e_sz[Y]; - emit_e_sz_2d[Y] = emit_e_sz[Z]; + emit_e_sz_2d[X] = scn.emit_e_sz[Y]; + emit_e_sz_2d[Y] = scn.emit_e_sz[Z]; I = compute_I(pos_emit_e_2d, emit_e_sz_2d); compute_grad_I(pos_emit_e_2d, emit_e_sz_2d, grad_I_2d); grad_I[X] = 0; @@ -284,8 +324,8 @@ realisation - rho * proj_V_s.beta * proj_u_e.beta * d3_dot(grad_I, proj_u_e.u) + rho * proj_V_e.beta * d3_dot(grad_I, proj_V_e.u); - w = I*emit_s_sz[X]*emit_s_sz[Y]*PI; /* Weight */ - s = S*emit_s_sz[X]*emit_s_sz[Y]*PI; /* Sensib */ + w = I*scn.emit_s_sz[X]*scn.emit_s_sz[Y]*PI; /* Weight */ + s = S*scn.emit_s_sz[X]*scn.emit_s_sz[Y]*PI; /* Sensib */ exit: SMC_DOUBLEN(weights)[WEIGHT] = w;