commit d8a52a4bb3bb0cc8653ba75080973c0e51581121
parent 1db74b6f5d469c3ac5f4020b5ded37637c1fe393
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 9 Jul 2024 11:17:46 +0200
BugFix components grouping stage
When grouping connex components, a radius is used to limit the number of
candidates. This radius is determinated in an initial step.
The bug was that the radius calculation could return a too small radius
that prevented considering some components, resulting in wrong grouping.
Diffstat:
1 file changed, 204 insertions(+), 76 deletions(-)
diff --git a/src/senc3d_scene_analyze.c b/src/senc3d_scene_analyze.c
@@ -67,12 +67,14 @@ const struct cc_descriptor CC_DESCRIPTOR_NULL = CC_DESCRIPTOR_NULL__;
* If ray.normal < threshold we suspect accuracy could be a problem */
#define DOT_THRESHOLD 0.0001f
-struct filter_ctx0 {
- component_id_t origin_component;
- struct darray_triangle_comp* triangles_comp;
+enum ctx_type {
+ CTX0,
+ CTX1,
+ CTX2
};
-struct filter_ctx1 {
+struct filter_ctx {
+ enum ctx_type type;
struct senc3d_scene* scn;
struct s3d_scene_view* view;
component_id_t origin_component;
@@ -80,34 +82,14 @@ struct filter_ctx1 {
struct darray_ptr_component_descriptor* components;
/* Tmp data used across filter calls */
double current_6volume;
+ int cpt;
float s;
- /* Result of hit */
+ /* Result of CTX1 hit */
component_id_t hit_component;
float hit_dir[3], hit_dist;
struct s3d_primitive hit_prim;
};
-struct filter_ctx2 {
- struct darray_triangle_comp* triangles_comp;
- int cpt;
- component_id_t component;
-};
-
-enum fctx_type {
- FCTX0,
- FCTX1,
- FCTX2
-};
-
-struct filter_ctx {
- enum fctx_type type;
- union {
- struct filter_ctx0 ctx0;
- struct filter_ctx1 ctx1;
- struct filter_ctx2 ctx2;
- } c;
-};
-
#define HTABLE_NAME overlap
#define HTABLE_KEY trg_id_t
#define HTABLE_DATA char
@@ -129,7 +111,7 @@ static FINLINE int
is_component_inside
(struct cc_descriptor* cc1,
struct cc_descriptor* cc2,
- struct filter_ctx1* ctx)
+ struct filter_ctx* ctx)
{
int i;
side_id_t side;
@@ -184,13 +166,13 @@ is_component_inside
d3_divd(pt, pt, 3);
f3_set_d3(org, pt);
/* Trace a ray and count intersections with component c */
- ctx2.type = FCTX2;
- ctx2.c.ctx2.triangles_comp = ctx->triangles_comp;
- ctx2.c.ctx2.cpt = 0;
- ctx2.c.ctx2.component = cc1->cc_id;
+ ctx2.type = CTX2;
+ ctx2.triangles_comp = ctx->triangles_comp;
+ ctx2.cpt = 0;
+ ctx2.origin_component = cc1->cc_id;
S3D(scene_view_trace_ray(ctx->view, org, dir, rg, &ctx2, &hit));
/* cc2 is not inside cc1 if cpt is even */
- if(ctx2.c.ctx2.cpt % 2 == 0) return 0;
+ if(ctx2.cpt % 2 == 0) return 0;
return 1;
}
@@ -256,68 +238,203 @@ self_hit_filter
void* ray_data,
void* filter_data)
{
- struct filter_ctx* fctx_ = ray_data;
+ struct filter_ctx* ctx = ray_data;
(void)ray_org; (void)ray_range; (void)filter_data;
- ASSERT(fctx_);
+ ASSERT(ctx);
ASSERT(hit->uv[0] == CLAMP(hit->uv[0], 0, 1));
ASSERT(hit->uv[1] == CLAMP(hit->uv[1], 0, 1));
- switch (fctx_->type) {
+ switch (ctx->type) {
default: FATAL("Invalid");
- case FCTX2: {
+ case CTX2: {
/* The filter is used to count the hits on some component along an
* infinite ray */
- struct filter_ctx2* ctx2 = &fctx_->c.ctx2;
const struct triangle_comp* trg_comp;
const component_id_t* hit_comp;
ASSERT(hit->prim.prim_id
- < darray_triangle_comp_size_get(ctx2->triangles_comp));
- trg_comp = darray_triangle_comp_cdata_get(ctx2->triangles_comp);
+ < darray_triangle_comp_size_get(ctx->triangles_comp));
+ trg_comp = darray_triangle_comp_cdata_get(ctx->triangles_comp);
hit_comp = trg_comp[hit->prim.prim_id].component;
- if(hit_comp[SENC3D_FRONT] == ctx2->component
- || hit_comp[SENC3D_BACK] == ctx2->component)
+ if(hit_comp[SENC3D_FRONT] == ctx->origin_component
+ || hit_comp[SENC3D_BACK] == ctx->origin_component)
{
- ctx2->cpt++;
+ ctx->cpt++;
}
return 1; /* Reject to continue counting */
}
- case FCTX0: {
+ case CTX0: {
/* This filter is called from a closest point query from a point belonging
* to origin_component. The returned hit is used to determine the search
- * radius for FCTX1 main computation. */
- struct filter_ctx0* ctx = &fctx_->c.ctx0;
+ * radius for CTX1 main computation. */
const struct triangle_comp*
trg_comp = darray_triangle_comp_cdata_get(ctx->triangles_comp);
const component_id_t*
hit_comp = trg_comp[hit->prim.prim_id].component;
const component_id_t oc = ctx->origin_component;
+ vrtx_id_t other_id;
+ struct cc_descriptor* const* comp_descriptors
+ = darray_ptr_component_descriptor_cdata_get(ctx->components);
+ size_t compsz = darray_ptr_component_descriptor_size_get(ctx->components);
+ const union double3*
+ vertices = darray_position_cdata_get(&ctx->scn->vertices);
+ const double org_z = vertices[comp_descriptors[oc]->max_z_vrtx_id].pos.z;
+ float s = 0, hit_normal[3], rdir[3];
+ enum senc3d_side hit_side;
+ const int log_components =
+ ctx->scn->convention & SENC3D_LOG_COMPONENTS_INFORMATION;
ASSERT(hit->prim.prim_id
< darray_triangle_comp_size_get(ctx->triangles_comp));
+ ASSERT(hit_comp[SENC3D_FRONT] < compsz);
+ ASSERT(hit_comp[SENC3D_BACK] < compsz);
+ /* Hit acceptance must be coherent at CTX0 and FTCX1 stages:
+ * the search radius as found at stage CTX0 must include the hit that
+ * stage CTX1 will select using an infinite radius */
if(hit_comp[SENC3D_FRONT] == oc || hit_comp[SENC3D_BACK] == oc) {
- /* Self hit */
- return 1; /* Reject */
+ return 1; /* Self hit, reject */
}
+ if(hit->distance == 0) {
+ /* origin component is in contact with some other components
+ * We will need further exploration to know if they should be considered
+ * Accepting hit at distance 0 => radius is definitively 0 */
+ int n;
+
+ /* If same component, process only once */
+ FOR_EACH(n, 0, (hit_comp[SENC3D_FRONT] == hit_comp[SENC3D_BACK] ? 1 : 2)) {
+ const enum senc3d_side sides[2] = { SENC3D_FRONT, SENC3D_BACK };
+ component_id_t c = hit_comp[sides[n]];
+ ASSERT(c < darray_ptr_component_descriptor_size_get(ctx->components));
+ if(comp_descriptors[c]->is_outer_border) {
+ /* The inner component we are trying to link can only be linked to
+ * an outer component if it is inside */
+ if(!is_component_inside(comp_descriptors[c],
+ comp_descriptors[oc], ctx))
+ {
+ continue;
+ }
- if(hit->distance > 0 && ray_dir[2] <= 0) {
- return 1; /* Not upward */
+ if(log_components) {
+ #pragma omp critical
+ printf("Component #%u: decreasing search radius "
+ "(R=%g, n=%g,%g,%g, components: %u, %u)\n",
+ oc, hit->distance, SPLIT3(hit->normal),
+ hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
+ }
+ return 0;
+ } else {
+ /* c is an inner component */
+ vrtx_id_t c_z_id;
+ /* The inner component we are trying to link can only be linked to
+ * another inner component if (at least partly) above and not inside */
+ c_z_id = comp_descriptors[c]->max_z_vrtx_id;
+ ASSERT(c_z_id < darray_position_size_get(&ctx->scn->vertices));
+ ASSERT(vertices[c_z_id].pos.z >= org_z);
+ if(vertices[c_z_id].pos.z == org_z) {
+ continue; /* Not above */
+ }
+ if(is_component_inside(comp_descriptors[c],
+ comp_descriptors[oc], ctx))
+ {
+ continue; /* Inside */
+ }
+ if(log_components) {
+ #pragma omp critical
+ printf("Component #%u: decreasing search radius "
+ "(R=%g, n=%g,%g,%g, components: %u, %u)\n",
+ oc, hit->distance, SPLIT3(hit->normal),
+ hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
+ }
+ return 0;
+ }
+ }
+ return 1;
+ }
+
+ ASSERT(hit->distance > 0);
+ if(hit_comp[SENC3D_FRONT] == hit_comp[SENC3D_BACK]) {
+ /* Hit component is known, check if above */
+ other_id = comp_descriptors[hit_comp[SENC3D_FRONT]]->max_z_vrtx_id;
+ ASSERT(other_id < darray_position_size_get(&ctx->scn->vertices));
+ if(vertices[other_id].pos.z <= org_z) {
+ return 1;
+ }
+ if(log_components) {
+ #pragma omp critical
+ printf("Component #%u: decreasing search radius "
+ "(R=%g, n=%g,%g,%g, components: %u, %u)\n",
+ oc, hit->distance, SPLIT3(hit->normal),
+ hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
+ }
+ return 0;
}
- return 0; /* Keep*/
+ /* Compute hit side */
+ /* For s to be comparable, vectors must be normalized */
+ f3_normalize(hit_normal, hit->normal);
+ f3_normalize(rdir, ray_dir);
+ s = f3_dot(rdir, hit_normal); /* Can be NaN for tiny distances */
+ if(isnan(s)) {
+ /* Try to fix it */
+ f3_divf(rdir, ray_dir, hit->distance);
+ f3_normalize(rdir, rdir);
+ s = f3_dot(rdir, hit_normal);
+ ASSERT(!isnan(s));
+ }
+
+ if(fabsf(s) < DOT_THRESHOLD) {
+ /* We cannot know for sure which side to consider */
+ vrtx_id_t i1 = comp_descriptors[hit_comp[SENC3D_FRONT]]->max_z_vrtx_id;
+ vrtx_id_t i2 = comp_descriptors[hit_comp[SENC3D_BACK]]->max_z_vrtx_id;
+ double possible_z = MMIN(vertices[i1].pos.z, vertices[i2].pos.z);
+ if(possible_z > org_z) {
+ /* Both components are above origin component => keep */
+ if(log_components) {
+ #pragma omp critical
+ printf("Component #%u: decreasing search radius "
+ "(R=%g, n=%g,%g,%g, components: %u, %u)\n",
+ oc, hit->distance, SPLIT3(hit->normal),
+ hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
+ }
+ return 0;
+ }
+ /* Cannot be sure => the safest choice is to reject */
+ return 1;
+ }
+ /* Determine which side was hit */
+ hit_side =
+ ((s < 0) /* Facing geometrical normal of hit */
+ == ((ctx->scn->convention & SENC3D_CONVENTION_NORMAL_FRONT) != 0))
+ /* Warning: following Embree 2 convention for geometrical normals,
+ * the Star3D hit normal is left-handed while star-enclosures-3d uses
+ * right-handed convention */
+ ? SENC3D_BACK : SENC3D_FRONT;
+ other_id = comp_descriptors[hit_comp[hit_side]]->max_z_vrtx_id;
+ ASSERT(other_id < darray_position_size_get(&ctx->scn->vertices));
+ if(vertices[other_id].pos.z <= org_z) {
+ return 1;
+ }
+ if(log_components) {
+ #pragma omp critical
+ printf("Component #%u: decreasing search radius "
+ "(R=%g, n=%g,%g,%g, components: %u, %u)\n",
+ oc, hit->distance, SPLIT3(hit->normal),
+ hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
+ }
+ return 0;
}
- case FCTX1: {
+ case CTX1: {
/* This filter is called from a closest point query from a point belonging
* to origin_component. The returned hit is used to determine a component
* to which origin_component is linked. At a later stage the algorithm
* process linked components to determine their relative inclusions.
*
* This filter is called with a search distance that has been ajusted in
- * FCTX0 filter. This distance must be left unchanged to ensure visiting
+ * CTX0 filter. This distance must be left unchanged to ensure visiting
* all the surfaces at the determined distance: allways reject hits to
* avoid decreasing search distance.
*
@@ -330,7 +447,6 @@ self_hit_filter
* (greater volume needed), or they can be disjoint, with (at least) ray_org
* as a common vertex (they can also partially intersect, but this is invalid
* and remains undetected by star enclosures). */
- struct filter_ctx1* ctx = &fctx_->c.ctx1;
struct cc_descriptor* const* comp_descriptors
= darray_ptr_component_descriptor_cdata_get(ctx->components);
const struct triangle_comp*
@@ -362,6 +478,9 @@ self_hit_filter
hit_comp[SENC3D_FRONT], hit_comp[SENC3D_BACK]);
}
+ /* Has CTX1 filter do not reduce search radius, hits can be processed
+ * that are at farther distance than the current best candidate: we need
+ * to reject them here */
if(hit->distance > ctx->hit_dist) {
/* No improvement */
if(log_components) {
@@ -372,6 +491,9 @@ self_hit_filter
return 1;
}
+ /* Hit acceptance must be coherent at CTX0 and FTCX1 stages:
+ * the search radius as found at stage CTX0 must include the hit that
+ * stage CTX1 will select using an infinite radius */
if(hit_comp[SENC3D_FRONT] == oc || hit_comp[SENC3D_BACK] == oc) {
/* Self hit */
if(log_components) {
@@ -425,7 +547,7 @@ self_hit_filter
* we are looking for is the smallest one (to manage outer component
* inside another outer component). */
if((ctx->hit_component != COMPONENT_NULL__
- && !comp_descriptors[ctx->hit_component]->is_outer_border )
+ && !comp_descriptors[ctx->hit_component]->is_outer_border)
|| v < ctx->current_6volume) {
ctx->hit_component = c;
ctx->current_6volume = v;
@@ -462,7 +584,7 @@ self_hit_filter
printf("Component #%u: hit inner component #%u\n", oc, c);
}
if(ctx->hit_component != COMPONENT_NULL__
- && comp_descriptors[ctx->hit_component]->is_outer_border )
+ && comp_descriptors[ctx->hit_component]->is_outer_border)
{
if(log_components) {
#pragma omp critical
@@ -528,8 +650,7 @@ self_hit_filter
if(vertices[other_id].pos.z <= org_z) {
if(log_components) {
#pragma omp critical
- printf("Component #%u: 2 sides, not (even in part) above: reject\n",
- oc);
+ printf("Component #%u: 2 sides, not above: reject\n", oc);
}
return 1;
}
@@ -548,7 +669,7 @@ self_hit_filter
}
/* Compute hit side */
- /* For s to be comparable, vectors must be normailzed */
+ /* For s to be comparable, vectors must be normalized */
f3_normalize(hit_normal, hit->normal);
f3_normalize(rdir, ray_dir);
s = f3_dot(rdir, hit_normal); /* Can be NaN for tiny distances */
@@ -580,9 +701,10 @@ self_hit_filter
vrtx_id_t i2 = comp_descriptors[hit_comp[SENC3D_BACK]]->max_z_vrtx_id;
double possible_z = MMAX(vertices[i1].pos.z, vertices[i2].pos.z);
if(possible_z <= org_z) {
+ /* None of the components are above origin component => reject */
if(log_components) {
#pragma omp critical
- printf("Component #%u: tiny s, not (even in part) above: reject\n",
+ printf("Component #%u: none of the components above: reject\n",
oc);
}
return 1;
@@ -1246,13 +1368,16 @@ group_connex_components
scene_cpt++;
}
- ctx0.type = FCTX0;
- ctx0.c.ctx0.triangles_comp = triangles_comp;
- ctx1.type = FCTX1;
- ctx1.c.ctx1.scn = scn;
- ctx1.c.ctx1.view = s3d_view;
- ctx1.c.ctx1.triangles_comp = triangles_comp;
- ctx1.c.ctx1.components = connex_components;
+ ctx0.type = CTX0;
+ ctx0.scn = scn;
+ ctx0.view = s3d_view;
+ ctx0.triangles_comp = triangles_comp;
+ ctx0.components = connex_components;
+ ctx1.type = CTX1;
+ ctx1.scn = scn;
+ ctx1.view = s3d_view;
+ ctx1.triangles_comp = triangles_comp;
+ ctx1.components = connex_components;
*res = s3d_scene_view_get_aabb(s3d_view, lower, upper);
if(*res != RES_OK) goto end;
#pragma omp for schedule(dynamic)
@@ -1294,8 +1419,11 @@ group_connex_components
rrr[i] = MMIN(origin[i] - lower[i], upper[i] - origin[i]);
}
rrr[2] = upper[2] - origin[2];
- r = f3_len(rrr) + FLT_EPSILON; /* Ensure r > 0 */
- ctx0.c.ctx0.origin_component = cc->cc_id;
+ r = f3_len(rrr) * (1 + FLT_EPSILON) + FLT_EPSILON; /* Ensure r > 0 */
+ ctx0.origin_component = cc->cc_id;
+ ctx0.current_6volume = DBL_MAX;
+ ctx0.hit_dist = FLT_MAX;
+ ctx0.hit_component = COMPONENT_NULL__;
tmp_res = s3d_scene_view_closest_point(s3d_view, origin, r, &ctx0, &hit);
if(tmp_res != RES_OK) {
*res = tmp_res;
@@ -1313,10 +1441,10 @@ group_connex_components
}
/* Second step is to determine which component faces component #c */
- ctx1.c.ctx1.origin_component = cc->cc_id;
- ctx1.c.ctx1.current_6volume = DBL_MAX;
- ctx1.c.ctx1.hit_dist = FLT_MAX;
- ctx1.c.ctx1.hit_component = COMPONENT_NULL__;
+ ctx1.origin_component = cc->cc_id;
+ ctx1.current_6volume = DBL_MAX;
+ ctx1.hit_dist = FLT_MAX;
+ ctx1.hit_component = COMPONENT_NULL__;
/* New search radius is hit.distance + some margin to cope with numerical
* issues (and r==0 is an error) */
r = hit.distance * (1 + FLT_EPSILON) + FLT_EPSILON;
@@ -1331,10 +1459,10 @@ group_connex_components
*res = tmp_res;
continue;
}
- /* As FCTX1 filter rejects any hit, do not rely on hit but use result as
+ /* As CTX1 filter rejects any hit, do not rely on hit but use result as
* stored in ctx1 */
/* If no hit is accepted, the component is facing an infinite medium */
- if(ctx1.c.ctx1.hit_dist == FLT_MAX) {
+ if(ctx1.hit_dist == FLT_MAX) {
cc->cc_group_root = CC_GROUP_ROOT_INFINITE;
cc->enclosure_id = 0;
if(log_components) {
@@ -1343,7 +1471,7 @@ group_connex_components
}
continue;
}
- else if(ctx1.c.ctx1.hit_component == COMPONENT_NULL__) {
+ else if(ctx1.hit_component == COMPONENT_NULL__) {
/* The selected triangle was nearly parallel to the line of sight:
* FRONT/BACK discrimination was not reliable enough and should be done
* differently. */
@@ -1354,7 +1482,7 @@ group_connex_components
continue;
} else {
/* If hit, group this component */
- cc->cc_group_root = ctx1.c.ctx1.hit_component;
+ cc->cc_group_root = ctx1.hit_component;
ASSERT(cc->cc_group_root < cc_count);
if(log_components) {
#pragma omp critical
