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 9b4f6e83f6ff79f2bdc2788b73583ffcf0a96022
parent 32b385e608b266f3de9d0114c5bba1ef6bc27560
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon,  3 May 2021 19:15:48 +0200

Fix the ray-tracing filtering function

Diffstat:
Msrc/sgs_compute_4v_s.c | 24+++++++++++++-----------
Msrc/sgs_geometry.c | 96++++++++++++++++++++++---------------------------------------------------------
Msrc/sgs_geometry.h | 27+++++++--------------------
Msrc/sgs_geometry_box.c | 6+++---
Msrc/sgs_geometry_c.h | 3++-
Msrc/sgs_geometry_slope.c | 6+++---
Msrc/sgs_geometry_step.c | 6+++---
7 files changed, 58 insertions(+), 110 deletions(-)

diff --git a/src/sgs_compute_4v_s.c b/src/sgs_compute_4v_s.c @@ -55,24 +55,25 @@ realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx double vec[3] = {0, 0, 0}; double pos[3] = {0, 0, 0}; double dir[3] = {0, 0, 0}; - double range[2] = {0, 0}; + double range[2] = {0, DBL_MAX}; - struct sgs_s3d_position* pos_from = NULL; + enum sgs_surface_type surface_from = SGS_SURFACE_NONE; struct sgs* sgs = ctx; double len = 0; - int sampling_mask = 0; + int emissive_surfaces = 0; res_T res = RES_OK; ASSERT(weight && rng && ctx); (void)ithread; - sampling_mask = sgs_geometry_get_sampling_mask(sgs->geom); + /* Retrieve the mask of the emissive surfaces */ + emissive_surfaces = sgs_geometry_get_sampling_mask(sgs->geom); /* Sample the emissive surface */ sgs_geometry_sample(sgs->geom, rng, &frag); - pos_from = &frag.s3d_pos; + surface_from = frag.surface; /* Cosine weighted sampling of the hemisphere around the sampled normal */ ssp_ran_hemisphere_cos(rng, frag.normal, dir, NULL); @@ -85,15 +86,16 @@ realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx range[1] = INF; for(;;) { - sgs_geometry_trace_ray(sgs->geom, pos, dir, range, pos_from, &hit); + sgs_geometry_trace_ray(sgs->geom, pos, dir, range, surface_from, &hit); if(SGS_HIT_NONE(&hit)) { /* Unexpected error */ res = RES_BAD_OP; goto error; } + /* Update the path length */ len += hit.distance; - if(((int)hit.surface_type & sampling_mask) != 0) + if(((int)hit.surface & emissive_surfaces) != 0) break; /* Stop random walk */ /* Move to the intersection */ @@ -104,7 +106,7 @@ realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx d3_minus(dir, dir); /* Ensure reflect convention */ reflect(dir, dir, hit.normal); - pos_from = &hit.s3d_pos; + surface_from = hit.surface; } exit: SMC_DOUBLE(weight) = len; @@ -153,9 +155,9 @@ compute_4v_s(struct sgs* sgs) goto error; } - /* Retrieve the estimated result */ + /* Retrieve the estimated value */ SMC(estimator_get_status(estimator, &status)); - if(status.NF == integrator.max_failures) { + if(status.NF >= integrator.max_failures) { sgs_log_err(sgs, "Too many failures: %lu\n", (unsigned long)status.NF); res = RES_BAD_OP; goto error; @@ -164,7 +166,7 @@ compute_4v_s(struct sgs* sgs) /* Fetch the reference values */ V = sgs_geometry_compute_volume(sgs->geom); S = sgs_geometry_compute_sampling_area(sgs->geom); - sgs_log(sgs, "V = %g; S = %g\n", V, S); + sgs_log(sgs, "V = %g; S = %g\n", V, S); /* Print the result */ sgs_log(sgs, "4V/S = %g ~ %g +/- %g; #failures = %lu/%lu\n", diff --git a/src/sgs_geometry.c b/src/sgs_geometry.c @@ -31,9 +31,10 @@ struct hit_filter_context { float range[2]; - struct sgs_s3d_position pos_from; + enum sgs_surface_type surface_from; + const int* prim2surface; /* Map a primitive id to its surface type */ }; -#define HIT_FILTER_CONTEXT_NULL__ {{0,FLT_MAX}, SGS_S3D_POSITION_NULL__} +#define HIT_FILTER_CONTEXT_NULL__ {{0,FLT_MAX}, SGS_SURFACE_NONE, NULL} static const struct hit_filter_context HIT_FILTER_CONTEXT_NULL = HIT_FILTER_CONTEXT_NULL__; @@ -48,55 +49,6 @@ static const struct mesh MESH_NULL = {NULL, NULL, 0, 0}; /******************************************************************************* * Helper functions ******************************************************************************/ -/* Check that the barycentric coordinates defines a position near of an edge. - * To do so test that at least one barycentric coordinate is roughly equal to 0 - * or 1. */ -static FINLINE int -is_on_edge(const float uv[2]) -{ - const float on_edge_eps = 1.e-4f; - float w; - ASSERT(uv); - w = 1.f - uv[0] - uv[1]; - return eq_epsf(uv[0], 0.f, on_edge_eps) - || eq_epsf(uv[0], 1.f, on_edge_eps) - || eq_epsf(uv[1], 0.f, on_edge_eps) - || eq_epsf(uv[1], 1.f, on_edge_eps) - || eq_epsf(w, 0.f, on_edge_eps) - || eq_epsf(w, 1.f, on_edge_eps); -} - -/* Return 1 if the submitted hit is actually a self hit, i.e. if the ray starts - * the primitive from which it starts */ -static INLINE int -self_hit - (const struct s3d_hit* hit, - const float ray_org[3], - const float ray_dir[3], - const float ray_range[2], - const struct sgs_s3d_position* pos_from) -{ - ASSERT(hit && pos_from); - (void)ray_org, (void)ray_dir, (void)ray_range; - - /* The hit primitive is the one from which the the ray starts. Discard it */ - if(S3D_PRIMITIVE_EQ(&hit->prim, &pos_from->prim)) - return 1; - - /* If the targeted point is near of the origin, check that it lies on an - * edge/vertex shared by the 2 primitives. In such situation, we assume that - * it is a self intersection and discard this hit */ - if(!S3D_PRIMITIVE_EQ(&pos_from->prim, &S3D_PRIMITIVE_NULL) - && eq_epsf(hit->distance, 0, 1.e-1f) - && is_on_edge(pos_from->st) - && is_on_edge(hit->uv)) { - return 1; - } - - /* No self hit */ - return 0; -} - static int hit_filter (const struct s3d_hit* hit, @@ -106,7 +58,8 @@ hit_filter void* filter_data) { const struct hit_filter_context* ctx = ray_data; - (void)filter_data; + enum sgs_surface_type surface; + (void)ray_org, (void)ray_dir, (void)filter_data; if(!ctx) /* Nothing to do */ return 0; @@ -117,8 +70,10 @@ hit_filter if(hit->distance <= ctx->range[0] || hit->distance >= ctx->range[1]) return 1; - if(self_hit(hit, ray_org, ray_dir, ctx->range, &ctx->pos_from)) - return 1; /* Discard this hit */ + /* The hit surface is the one from which the the ray starts. Discard it */ + surface = ctx->prim2surface[hit->prim.prim_id]; + if(surface == ctx->surface_from) + return 1; return 0; /* No filtering */ } @@ -202,7 +157,8 @@ release_geometry(ref_T* ref) if(geom->view_rt) S3D(scene_view_ref_put(geom->view_rt)); darray_double_release(&geom->verts); darray_size_t_release(&geom->tris); - darray_int_release(&geom->tris_type); + darray_int_release(&geom->tris_rt_type); + darray_int_release(&geom->tris_sp_type); sgs = geom->sgs; MEM_RM(sgs_get_allocator(sgs), geom); @@ -256,7 +212,7 @@ sgs_geometry_trace_ray const double ray_org[3], const double ray_dir[3], const double ray_range[2], - const struct sgs_s3d_position* pos_from, + const enum sgs_surface_type surface_from, struct sgs_hit* out_hit) { struct hit_filter_context filter_ctx = HIT_FILTER_CONTEXT_NULL; @@ -278,7 +234,8 @@ sgs_geometry_trace_ray /* Setup the hit filter data */ filter_ctx.range[0] = range[0]; filter_ctx.range[1] = range[1]; - filter_ctx.pos_from = pos_from ? *pos_from : SGS_S3D_POSITION_NULL; + filter_ctx.surface_from = surface_from; + filter_ctx.prim2surface = darray_int_cdata_get(&geom->tris_rt_type); /* Trace the ray */ S3D(scene_view_trace_ray(geom->view_rt, org, dir, range, &filter_ctx, &hit)); @@ -289,10 +246,7 @@ sgs_geometry_trace_ray d3_set_f3(out_hit->normal, hit.normal); d3_normalize(out_hit->normal, out_hit->normal); out_hit->distance = hit.distance; - out_hit->surface_type = darray_int_cdata_get(&geom->tris_type)[hit.prim.prim_id]; - out_hit->s3d_pos.prim = hit.prim; - out_hit->s3d_pos.st[0] = hit.uv[0]; - out_hit->s3d_pos.st[1] = hit.uv[1]; + out_hit->surface = darray_int_cdata_get(&geom->tris_rt_type)[hit.prim.prim_id]; } } @@ -323,10 +277,7 @@ sgs_geometry_sample d3_set_f3(frag->position, position.value); d3_set_f3(frag->normal, normal.value); d3_normalize(frag->normal, frag->normal); - frag->surface_type = darray_int_cdata_get(&geom->tris_type)[prim.prim_id]; - frag->s3d_pos.prim = prim; - frag->s3d_pos.st[0] = st[0]; - frag->s3d_pos.st[1] = st[1]; + frag->surface = darray_int_cdata_get(&geom->tris_sp_type)[prim.prim_id]; } res_T @@ -378,7 +329,8 @@ sgs_geometry_dump_vtk(const struct sgs_geometry* geom, FILE* stream) FPRINTF("SCALARS Triangle_Type integer 1\n", ARG0()); FPRINTF("LOOKUP_TABLE default\n", ARG0()); FOR_EACH(i, 0, ntris) { - const enum sgs_surface_type type = darray_int_cdata_get(&geom->tris_type)[i]; + const enum sgs_surface_type type = + darray_int_cdata_get(&geom->tris_rt_type)[i]; FPRINTF("%i\n", ARG1(type)); } @@ -386,7 +338,7 @@ sgs_geometry_dump_vtk(const struct sgs_geometry* geom, FILE* stream) FPRINTF("SCALARS Sampling_Surface integer 1\n", ARG0()); FPRINTF("LOOKUP_TABLE Sampling_Type\n", ARG0()); FOR_EACH(i, 0, ntris) { - const int type = darray_int_cdata_get(&geom->tris_type)[i]; + const int type = darray_int_cdata_get(&geom->tris_rt_type)[i]; FPRINTF("%i\n", ARG1((type & geom->sampling_mask) != 0)); } FPRINTF("LOOKUP_TABLE Sampling_Type 2\n", ARG0()); @@ -428,7 +380,8 @@ geometry_create geom->sgs = sgs; darray_double_init(allocator, &geom->verts); darray_size_t_init(allocator, &geom->tris); - darray_int_init(allocator, &geom->tris_type); + darray_int_init(allocator, &geom->tris_rt_type); + darray_int_init(allocator, &geom->tris_sp_type); geom->type = type; geom->sampling_mask = sampling_mask; @@ -504,11 +457,12 @@ geometry_setup_view_sp(struct sgs_geometry* geom, const int sampling_mask) /* Fetch the geometry data */ ntris = geometry_get_ntriangles(geom); ids = darray_size_t_cdata_get(&geom->tris); - types = darray_int_cdata_get(&geom->tris_type); + types = darray_int_cdata_get(&geom->tris_rt_type); /* Find the triangles to sample */ FOR_EACH(itri, 0, ntris) { if((types[itri] & sampling_mask) != 0) { + /* Register the triangle as a triangle to sample */ res = darray_size_t_push_back(&tris_sp, &ids[itri*3+0]); if(res != RES_OK) goto error; @@ -516,6 +470,10 @@ geometry_setup_view_sp(struct sgs_geometry* geom, const int sampling_mask) if(res != RES_OK) goto error; res = darray_size_t_push_back(&tris_sp, &ids[itri*3+2]); if(res != RES_OK) goto error; + + /* Save the surface type of the triangle */ + res = darray_int_push_back(&geom->tris_sp_type, &types[itri]); + if(res != RES_OK) goto error; } } diff --git a/src/sgs_geometry.h b/src/sgs_geometry.h @@ -41,7 +41,8 @@ enum sgs_surface_type { SGS_SURFACE_Z_SLOPE = BIT(6), /* +Z face representing a slope */ SGS_SURFACE_Z_LVL0 = BIT(7), /* +Z face representing the 1st level along X */ SGS_SURFACE_Z_LVL1 = BIT(8), /* +Z face representing the 2nd level along X */ - SGS_SURFACE_Z_STEP = BIT(9) /* +Z face representing the step between the 2 levels */ + SGS_SURFACE_Z_STEP = BIT(9), /* +Z face representing the step between the 2 levels */ + SGS_SURFACE_NONE = 0 }; /* A simple axis aligned bounding box @@ -115,37 +116,23 @@ struct sgs_geometry_step_args { static const struct sgs_geometry_step_args SGS_GEOMETRY_STEP_ARGS_DEFAULT = SGS_GEOMETRY_STEP_ARGS_DEFAULT__; -/* Position defined relatively to a Star-3D primitive */ -struct sgs_s3d_position { - struct s3d_primitive prim; - float st[2]; /* Parametric coordinates onto the primitive */ -}; -#define SGS_S3D_POSITION_NULL__ {S3D_PRIMITIVE_NULL__, {0,0}} -static const struct sgs_s3d_position SGS_S3D_POSITION_NULL = - SGS_S3D_POSITION_NULL__; - /* Position onto a geometry surface */ struct sgs_fragment { double normal[3]; /* Normalized */ double position[3]; - enum sgs_surface_type surface_type; + enum sgs_surface_type surface; /* Surface to which the fragment belongs */ - /* Position defined wrt the corresponding Star-3D primitive */ - struct sgs_s3d_position s3d_pos; }; -#define SGS_FRAGMENT_NULL__ {{0,0,0},{0,0,0},0,SGS_S3D_POSITION_NULL__} +#define SGS_FRAGMENT_NULL__ {{0,0,0},{0,0,0},0} static const struct sgs_fragment SGS_FRAGMENT_NULL = SGS_FRAGMENT_NULL__; /* Intersection along a ray */ struct sgs_hit { double normal[3]; /* Normalized */ double distance; /* Distance up to the intersection */ - enum sgs_surface_type surface_type; - - /* Position defined wrt the corresponding Star-3D primitive */ - struct sgs_s3d_position s3d_pos; + enum sgs_surface_type surface; /* Surface to which the hit belongs */ }; -#define SGS_HIT_NULL__ {{0,0,0},DBL_MAX,0,SGS_S3D_POSITION_NULL__} +#define SGS_HIT_NULL__ {{0,0,0},DBL_MAX,0} static const struct sgs_hit SGS_HIT_NULL = SGS_HIT_NULL__; /* Helper macro */ @@ -203,7 +190,7 @@ sgs_geometry_trace_ray const double ray_org[3], const double ray_dir[3], const double ray_range[2], - const struct sgs_s3d_position* pos_from, /* To avoid self hit. May be NULL */ + const enum sgs_surface_type surface_from, /* To avoid self hit */ struct sgs_hit* hit); /* Uniformly sample a point onto the faces of the geometry defined in the diff --git a/src/sgs_geometry_box.c b/src/sgs_geometry_box.c @@ -58,13 +58,13 @@ setup_box_mesh if(res != RES_OK) goto error; res = darray_size_t_resize(&geom->tris, 12/*#triangles*/*3/*#ids per triangle*/); if(res != RES_OK) goto error; - res = darray_int_resize(&geom->tris_type, 12/*#triangles*/); + res = darray_int_resize(&geom->tris_rt_type, 12/*#triangles*/); if(res != RES_OK) goto error; /* Fetch allocated memory space */ vtx = darray_double_data_get(&geom->verts); ids = darray_size_t_data_get(&geom->tris); - types = darray_int_data_get(&geom->tris_type); + types = darray_int_data_get(&geom->tris_rt_type); /* Fetch input args */ low = args->lower; @@ -107,7 +107,7 @@ exit: error: darray_double_clear(&geom->verts); darray_size_t_clear(&geom->tris); - darray_int_clear(&geom->tris_type); + darray_int_clear(&geom->tris_rt_type); goto exit; } diff --git a/src/sgs_geometry_c.h b/src/sgs_geometry_c.h @@ -32,7 +32,8 @@ struct s3d_scene_view; struct sgs_geometry { struct darray_double verts; struct darray_size_t tris; - struct darray_int tris_type; /* List of enum sgs_surface_type */ + struct darray_int tris_rt_type; /* List of enum sgs_surface_type */ + struct darray_int tris_sp_type; /* List of enum sgs_surface_type */ struct s3d_scene_view* view_rt; struct s3d_scene_view* view_sp; diff --git a/src/sgs_geometry_slope.c b/src/sgs_geometry_slope.c @@ -61,13 +61,13 @@ setup_slope_mesh if(res != RES_OK) goto error; res = darray_size_t_resize(&geom->tris, 12/*#triangles*/*3/*#ids per triangle*/); if(res != RES_OK) goto error; - res = darray_int_resize(&geom->tris_type, 12/*#triangles*/); + res = darray_int_resize(&geom->tris_rt_type, 12/*#triangles*/); if(res != RES_OK) goto error; /* Fetch allocated memory space */ vtx = darray_double_data_get(&geom->verts); ids = darray_size_t_data_get(&geom->tris); - types = darray_int_data_get(&geom->tris_type); + types = darray_int_data_get(&geom->tris_rt_type); /* Fetch input args */ low = args->lower; @@ -111,7 +111,7 @@ exit: error: darray_double_clear(&geom->verts); darray_size_t_clear(&geom->tris); - darray_int_clear(&geom->tris_type); + darray_int_clear(&geom->tris_rt_type); goto exit; } diff --git a/src/sgs_geometry_step.c b/src/sgs_geometry_step.c @@ -66,13 +66,13 @@ setup_step_mesh if(res != RES_OK) goto error; res = darray_size_t_resize(&geom->tris, 20/*#triangles*/*3/*#ids per triangle*/); if(res != RES_OK) goto error; - res = darray_int_resize(&geom->tris_type, 20/*#triangles*/); + res = darray_int_resize(&geom->tris_rt_type, 20/*#triangles*/); if(res != RES_OK) goto error; /* Fetch allocated memory space */ vtx = darray_double_data_get(&geom->verts); ids = darray_size_t_data_get(&geom->tris); - types = darray_int_data_get(&geom->tris_type); + types = darray_int_data_get(&geom->tris_rt_type); /* Fetch input args */ low = args->lower; @@ -160,7 +160,7 @@ exit: error: darray_double_clear(&geom->verts); darray_size_t_clear(&geom->tris); - darray_int_clear(&geom->tris_type); + darray_int_clear(&geom->tris_rt_type); goto exit; }