star-enclosures-3d

Extract enclosures from 3D geometry
git clone git://git.meso-star.fr/star-enclosures-3d.git
Log | Files | Refs | README | LICENSE

commit a4909ab21264861af9f526021fe0b5f38fbc5fc2
parent d2e828c0b3c22f97351b11c8c106da710d7d0129
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Fri, 20 Apr 2018 16:10:02 +0200

Fix a numeric accuracy concern.

Could fail to determine wether a component is inner/outer border.

Diffstat:
Msrc/senc_descriptor_c.h | 4+---
Msrc/senc_scene_analyze.c | 114++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
Msrc/senc_scene_analyze_c.h | 3++-
3 files changed, 81 insertions(+), 40 deletions(-)

diff --git a/src/senc_descriptor_c.h b/src/senc_descriptor_c.h @@ -32,7 +32,6 @@ struct triangle_comp { component_id_t component[2]; }; -#ifndef NDEBUG static void triangle_comp_init(struct mem_allocator* alloc, struct triangle_comp* trg) { int i; @@ -40,11 +39,10 @@ triangle_comp_init(struct mem_allocator* alloc, struct triangle_comp* trg) { ASSERT(trg); FOR_EACH(i, 0, 2) trg->component[i] = COMPONENT_NULL__; } -#define DARRAY_FUNCTOR_INIT triangle_comp_init -#endif #define DARRAY_NAME triangle_comp #define DARRAY_DATA struct triangle_comp +#define DARRAY_FUNCTOR_INIT triangle_comp_init #include <rsys/dynamic_array.h> struct triangle_enc { diff --git a/src/senc_scene_analyze.c b/src/senc_scene_analyze.c @@ -44,7 +44,7 @@ #include <stdlib.h> #define CC_DESCRIPTOR_NULL__ {\ - -1, VRTX_NULL__, 0, MEDIUM_NULL__,\ + CHAR_MAX, VRTX_NULL__, 0, MEDIUM_NULL__,\ CC_ID_NONE, CC_GROUP_ROOT_NONE, ENCLOSURE_NULL__,\ { TRG_NULL__, 0}\ } @@ -91,7 +91,7 @@ neighbour_cmp(const void* w1, const void* w2) return (a1 > a2) - (a1 < a2); } -static FINLINE void +static FINLINE res_T add_side_to_stack (const struct senc_scene* scn, struct darray_side_id* stack, @@ -103,8 +103,8 @@ add_side_to_stack && side_id < SIDE_MAX__ && side_id < 2 * scn->nutris); ASSERT((darray_side_id_size_get(stack) > 128) || !find_side_in_list(trgsides, stack, side_id, FLAG_WAITING_STACK)); - darray_side_id_push_back(stack, &side_id); trgsides[side_id].list_id = FLAG_WAITING_STACK; + return darray_side_id_push_back(stack, &side_id); } static side_id_t @@ -213,7 +213,9 @@ extract_connex_components struct mem_allocator* alloc; int64_t mm; struct darray_side_id stack; + struct darray_side_id ids_of_sides_around_max_z_vertex; const union double3* positions; + size_t sz, ii; #ifndef NDEBUG trg_id_t t; component_id_t c; @@ -225,6 +227,7 @@ extract_connex_components scn = desc->scene; positions = darray_position_cdata_get(&scn->vertices); darray_side_id_init(alloc, &stack); + darray_side_id_init(alloc, &ids_of_sides_around_max_z_vertex); #ifndef NDEBUG #pragma omp single @@ -251,6 +254,7 @@ extract_connex_components struct cc_descriptor* cc; /* Any not-already-used side is used as a starting point */ side_id_t first_side_not_in_component; + double max_z_nz; if(*res != RES_OK) continue; first_side_not_in_component @@ -265,7 +269,10 @@ extract_connex_components side_id_t crt_side_id = start_side_id; side_id_t last_side_id = start_side_id; double max_z = -DBL_MAX; + ASSERT(start_side_id == SIDE_NULL__ || start_side_id < 2 * scn->nutris); + + if(*res != RES_OK) break; if(start_side_id == SIDE_NULL__) break; /* start_side_id=SIDE_NULL__ => done! */ ASSERT(trgsides[start_side_id].list_id == FLAG_LIST_SIDE_LIST); @@ -303,12 +310,11 @@ extract_connex_components darray_triangle_comp_data_get(triangles_comp) + crt_trg_id; const struct triangle_tmp* const trg_tmp = darray_triangle_tmp_cdata_get(triangles_tmp_array) + crt_trg_id; - const union double3* vertices = - darray_position_cdata_get(&scn->vertices); - int use_max_z_vrtx = 0; ASSERT(crt_trg_id < scn->nutris); ASSERT(trgsides[crt_side_id].medium == m); + if(*res != RES_OK) break; + /* Record Zmax information * Keep track of the appropriate vertex of the connex component in order * to cast a ray at the component grouping step of the algorithm. @@ -317,44 +323,32 @@ extract_connex_components if(max_z < trg_tmp->max_z) { /* New best vertex */ max_z = trg_tmp->max_z; - cc->max_z_nz = 0; - /* Scan vertices to find one with z == max_z */ + /* New vertex: reset list of sides */ + darray_side_id_clear(&ids_of_sides_around_max_z_vertex); + + /* Select one vertex with z == max_z */ FOR_EACH(i, 0, 3) { if(max_z == positions[trg_in->vertice_id[i]].pos.z) { cc->max_z_vrtx_id = trg_in->vertice_id[i]; break; } - ASSERT(i < 3); /* One of the vertices matched */ } - use_max_z_vrtx = 1; + ASSERT(i < 3); /* Found one */ + /* List of sides using the vertex */ + *res = darray_side_id_push_back(&ids_of_sides_around_max_z_vertex, + &crt_side_id); } else if(max_z == trg_tmp->max_z) { /* Does this triangle use the currently selected max_z vertex? */ FOR_EACH(i, 0, 3) { if(cc->max_z_vrtx_id == trg_in->vertice_id[i]) { - use_max_z_vrtx = 1; + /* List of sides using the vertex */ + *res = darray_side_id_push_back(&ids_of_sides_around_max_z_vertex, + &crt_side_id); break; } } } - if(use_max_z_vrtx) { - /* One of the sides using the current best vertex. - * Contribute to the normal at the vertex. */ - double edge0[3], edge1[3], normal[3], norm; - - d3_sub(edge0, vertices[trg_in->vertice_id[1]].vec, - vertices[trg_in->vertice_id[0]].vec); - d3_sub(edge1, vertices[trg_in->vertice_id[2]].vec, - vertices[trg_in->vertice_id[0]].vec); - d3_cross(normal, edge0, edge1); - norm = d3_normalize(normal, normal); - ASSERT(norm); (void)norm; - - if(TRGSIDE_IS_FRONT(crt_side_id)) { - cc->max_z_nz += normal[2]; - } else { - cc->max_z_nz -= normal[2]; - } - } + if(*res != RES_OK) break; /* Record crt_side both as component and triangle level */ cc->side_count++; @@ -397,7 +391,7 @@ extract_connex_components log_err(desc->scene->dev, "Media: %lu VS %lu\n", (unsigned long)neighbour->medium, (unsigned long)crt_side->medium); *res = RES_BAD_ARG; - continue; + break; } if(neighbour->list_id == FLAG_LIST_COMPONENT) { /* Already processed */ @@ -409,24 +403,71 @@ extract_connex_components if(neighbour->list_id == FLAG_WAITING_STACK) { continue; /* Already processed */ } - add_side_to_stack(scn, &stack, trgsides, neighbour_id); + *res = add_side_to_stack(scn, &stack, trgsides, neighbour_id); + if(*res != RES_OK) break; } + if(*res != RES_OK) break; if(darray_side_id_size_get(&stack) == 0) break; /* Empty stack => connex component is done! */ crt_side_id = get_side_from_stack(trgsides, &stack); last_side_id = MMAX(last_side_id, crt_side_id); } + if(*res != RES_OK) break; /* Keep track of this new connex component */ cc->side_range.last = last_side_id; + + /* Compute the normal at the max_z vertex. */ + max_z_nz = 0; + sz = darray_side_id_size_get(&ids_of_sides_around_max_z_vertex); + FOR_EACH(ii, 0, sz) { + const side_id_t side_id = + darray_side_id_cdata_get(&ids_of_sides_around_max_z_vertex)[ii]; + const trg_id_t trg_id = TRGSIDE_2_TRG(side_id); + const struct triangle_in* trg_in = + darray_triangle_in_cdata_get(&scn->triangles_in) + trg_id; + const struct triangle_comp* trg_comp = + darray_triangle_comp_cdata_get(triangles_comp) + trg_id; + const union double3* vertices = + darray_position_cdata_get(&scn->vertices); + double edge0[3], edge1[3], normal[3], norm; + + /* To garanty that triangles with 2 sides in the component total to 0 + * regardless of numeric accuracy, we need to prevent them to + * contribute (remember than x + y - x - y can be non-zero). */ + ASSERT(trg_comp->component[SIDE_FRONT] == cc->cc_id + || trg_comp->component[SIDE_BACK] == cc->cc_id); + if (trg_comp->component[SIDE_FRONT] == trg_comp->component[SIDE_BACK]) + continue; + + d3_sub(edge0, vertices[trg_in->vertice_id[1]].vec, + vertices[trg_in->vertice_id[0]].vec); + d3_sub(edge1, vertices[trg_in->vertice_id[2]].vec, + vertices[trg_in->vertice_id[0]].vec); + d3_cross(normal, edge0, edge1); + norm = d3_normalize(normal, normal); + ASSERT(norm); (void) norm; + + /* Geometrical normal points toward the front side */ + if (TRGSIDE_IS_FRONT(side_id)) { + max_z_nz += normal[2]; + } + else { + max_z_nz -= normal[2]; + } + } + if(*res != RES_OK) break; + cc->is_outer_border = (max_z_nz < 0); + /* Need to synchronize connex_components growth as this global structure * is accessed by multipe threads */ #pragma omp critical { struct cc_descriptor** components; - size_t sz = darray_ptr_component_descriptor_size_get(connex_components); + sz = darray_ptr_component_descriptor_size_get(connex_components); if(sz <= cc->cc_id) { - res_T tmp_res = darray_ptr_component_descriptor_resize(connex_components, - 1 + cc->cc_id); + res_T tmp_res = + darray_ptr_component_descriptor_resize(connex_components, + 1 + cc->cc_id); if(tmp_res != RES_OK) *res = tmp_res; } if(*res == RES_OK) { @@ -441,6 +482,7 @@ extract_connex_components } /* No barrier here */ darray_side_id_release(&stack); + darray_side_id_release(&ids_of_sides_around_max_z_vertex); /* The first thread here creates the s3d view */ #pragma omp single nowait @@ -549,7 +591,7 @@ group_connex_components ASSERT(cc->cc_id == c); ASSERT(cc->cc_group_root == CC_GROUP_ID_NONE); - if(cc->max_z_nz < 0) { + if(cc->is_outer_border) { int64_t id; /* Don't need to cast a ray */ cc->cc_group_root = cc->cc_id; /* New group with self as root */ diff --git a/src/senc_scene_analyze_c.h b/src/senc_scene_analyze_c.h @@ -109,7 +109,8 @@ struct trgside { #define CC_GROUP_ROOT_INFINITE (COMPONENT_MAX__-1) #define CC_GROUP_ID_NONE COMPONENT_MAX__ struct cc_descriptor { - double max_z_nz; + /* Does this component is an outer border of an enclosure or an inner one? */ + char is_outer_border; vrtx_id_t max_z_vrtx_id; side_id_t side_count; medium_id_t medium;