diff --git a/xs/src/admesh/connect.c b/xs/src/admesh/connect.c index 6ce775f9f..56ebfa144 100644 --- a/xs/src/admesh/connect.c +++ b/xs/src/admesh/connect.c @@ -106,53 +106,42 @@ stl_check_facets_exact(stl_file *stl) { } } stl_free_edges(stl); + +#if 0 + printf("Number of faces: %d, number of manifold edges: %d, number of connected edges: %d, number of unconnected edges: %d\r\n", + stl->stats.number_of_facets, stl->stats.number_of_facets * 3, + stl->stats.connected_edges, stl->stats.number_of_facets * 3 - stl->stats.connected_edges); +#endif } static void stl_load_edge_exact(stl_file *stl, stl_hash_edge *edge, stl_vertex *a, stl_vertex *b) { - float diff_x; - float diff_y; - float diff_z; - float max_diff; - if (stl->error) return; - diff_x = ABS(a->x - b->x); - diff_y = ABS(a->y - b->y); - diff_z = ABS(a->z - b->z); - max_diff = STL_MAX(diff_x, diff_y); - max_diff = STL_MAX(diff_z, max_diff); - stl->stats.shortest_edge = STL_MIN(max_diff, stl->stats.shortest_edge); + { + float diff_x = ABS(a->x - b->x); + float diff_y = ABS(a->y - b->y); + float diff_z = ABS(a->z - b->z); + float max_diff = STL_MAX(diff_x, diff_y); + max_diff = STL_MAX(diff_z, max_diff); + stl->stats.shortest_edge = STL_MIN(max_diff, stl->stats.shortest_edge); + } - if(diff_x == max_diff) { - if(a->x > b->x) { - memcpy(&edge->key[0], a, sizeof(stl_vertex)); - memcpy(&edge->key[3], b, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], b, sizeof(stl_vertex)); - memcpy(&edge->key[3], a, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } - } else if(diff_y == max_diff) { - if(a->y > b->y) { - memcpy(&edge->key[0], a, sizeof(stl_vertex)); - memcpy(&edge->key[3], b, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], b, sizeof(stl_vertex)); - memcpy(&edge->key[3], a, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } + // Ensure identical vertex ordering of equal edges. + // This method is numerically robust. + if ((a->x != b->x) ? + (a->x < b->x) : + ((a->y != b->y) ? + (a->y < b->y) : + (a->z < b->z))) { + memcpy(&edge->key[0], a, sizeof(stl_vertex)); + memcpy(&edge->key[3], b, sizeof(stl_vertex)); } else { - if(a->z > b->z) { - memcpy(&edge->key[0], a, sizeof(stl_vertex)); - memcpy(&edge->key[3], b, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], b, sizeof(stl_vertex)); - memcpy(&edge->key[3], a, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } + memcpy(&edge->key[0], b, sizeof(stl_vertex)); + memcpy(&edge->key[3], a, sizeof(stl_vertex)); + edge->which_edge += 3; /* this edge is loaded backwards */ } } @@ -309,26 +298,17 @@ stl_check_facets_nearby(stl_file *stl, float tolerance) { static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge, stl_vertex *a, stl_vertex *b, float tolerance) { - float diff_x; - float diff_y; - float diff_z; - float max_diff; - unsigned vertex1[3]; - unsigned vertex2[3]; - - - diff_x = ABS(a->x - b->x); - diff_y = ABS(a->y - b->y); - diff_z = ABS(a->z - b->z); - max_diff = STL_MAX(diff_x, diff_y); - max_diff = STL_MAX(diff_z, max_diff); - - vertex1[0] = (unsigned)((a->x - stl->stats.min.x) / tolerance); - vertex1[1] = (unsigned)((a->y - stl->stats.min.y) / tolerance); - vertex1[2] = (unsigned)((a->z - stl->stats.min.z) / tolerance); - vertex2[0] = (unsigned)((b->x - stl->stats.min.x) / tolerance); - vertex2[1] = (unsigned)((b->y - stl->stats.min.y) / tolerance); - vertex2[2] = (unsigned)((b->z - stl->stats.min.z) / tolerance); + // Index of a grid cell spaced by tolerance. + uint32_t vertex1[3] = { + (uint32_t)((a->x - stl->stats.min.x) / tolerance), + (uint32_t)((a->y - stl->stats.min.y) / tolerance), + (uint32_t)((a->z - stl->stats.min.z) / tolerance) + }; + uint32_t vertex2[3] = { + (uint32_t)((b->x - stl->stats.min.x) / tolerance), + (uint32_t)((b->y - stl->stats.min.y) / tolerance), + (uint32_t)((b->z - stl->stats.min.z) / tolerance) + }; if( (vertex1[0] == vertex2[0]) && (vertex1[1] == vertex2[1]) @@ -337,33 +317,19 @@ stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge, return 0; } - if(diff_x == max_diff) { - if(a->x > b->x) { - memcpy(&edge->key[0], vertex1, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex2, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], vertex2, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex1, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } - } else if(diff_y == max_diff) { - if(a->y > b->y) { - memcpy(&edge->key[0], vertex1, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex2, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], vertex2, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex1, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } + // Ensure identical vertex ordering of edges, which vertices land into equal grid cells. + // This method is numerically robust. + if ((vertex1[0] != vertex2[0]) ? + (vertex1[0] < vertex2[0]) : + ((vertex1[1] != vertex2[1]) ? + (vertex1[1] < vertex2[1]) : + (vertex1[2] < vertex2[2]))) { + memcpy(&edge->key[0], vertex1, sizeof(stl_vertex)); + memcpy(&edge->key[3], vertex2, sizeof(stl_vertex)); } else { - if(a->z > b->z) { - memcpy(&edge->key[0], vertex1, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex2, sizeof(stl_vertex)); - } else { - memcpy(&edge->key[0], vertex2, sizeof(stl_vertex)); - memcpy(&edge->key[3], vertex1, sizeof(stl_vertex)); - edge->which_edge += 3; /* this edge is loaded backwards */ - } + memcpy(&edge->key[0], vertex2, sizeof(stl_vertex)); + memcpy(&edge->key[3], vertex1, sizeof(stl_vertex)); + edge->which_edge += 3; /* this edge is loaded backwards */ } return 1; } @@ -564,6 +530,33 @@ stl_change_vertices(stl_file *stl, int facet_num, int vnot, next_edge = pivot_vertex; } } +#if 0 + if (stl->facet_start[facet_num].vertex[pivot_vertex].x == new_vertex.x && + stl->facet_start[facet_num].vertex[pivot_vertex].y == new_vertex.y && + stl->facet_start[facet_num].vertex[pivot_vertex].z == new_vertex.z) + printf("Changing vertex %f,%f,%f: Same !!!\r\n", + new_vertex.x, new_vertex.y, new_vertex.z); + else { + if (stl->facet_start[facet_num].vertex[pivot_vertex].x != new_vertex.x) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex].x, + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].x), + new_vertex.x, + *reinterpret_cast<const int*>(&new_vertex.x)); + if (stl->facet_start[facet_num].vertex[pivot_vertex].y != new_vertex.y) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex].y, + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].y), + new_vertex.y, + *reinterpret_cast<const int*>(&new_vertex.y)); + if (stl->facet_start[facet_num].vertex[pivot_vertex].z != new_vertex.z) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex].z, + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].z), + new_vertex.z, + *reinterpret_cast<const int*>(&new_vertex.z)); + } +#endif stl->facet_start[facet_num].vertex[pivot_vertex] = new_vertex; vnot = stl->neighbors_start[facet_num].which_vertex_not[next_edge]; facet_num = stl->neighbors_start[facet_num].neighbor[next_edge]; diff --git a/xs/src/admesh/stl.h b/xs/src/admesh/stl.h index 3e498d929..e6e69c7f6 100644 --- a/xs/src/admesh/stl.h +++ b/xs/src/admesh/stl.h @@ -77,7 +77,10 @@ typedef struct { typedef struct stl_hash_edge { // Key of a hash edge: 2x binary copy of a floating point vertex. uint32_t key[6]; + // Index of a facet owning this edge. int facet_number; + // Index of this edge inside the facet with an index of facet_number. + // If this edge is stored backwards, which_edge is increased by 3. int which_edge; struct stl_hash_edge *next; } stl_hash_edge; diff --git a/xs/src/admesh/stlinit.c b/xs/src/admesh/stlinit.c index ff0e9227e..4e48ef48e 100644 --- a/xs/src/admesh/stlinit.c +++ b/xs/src/admesh/stlinit.c @@ -284,20 +284,6 @@ stl_read(stl_file *stl, int first_facet, int first) { stl->error = 1; return; } -#if 0 - // Report close to zero vertex coordinates. Due to the nature of the floating point numbers, - // close to zero values may be represented with singificantly higher precision than the rest of the vertices. - // It may be worth to round these numbers to zero during loading to reduce the number of errors reported - // during the STL import. - for (size_t j = 0; j < 3; ++ j) { - if (facet.vertex[j].x > -1e-12f && facet.vertex[j].x < 1e-12f) - printf("stl_read: facet %d.x = %e\r\n", j, facet.vertex[j].x); - if (facet.vertex[j].y > -1e-12f && facet.vertex[j].y < 1e-12f) - printf("stl_read: facet %d.y = %e\r\n", j, facet.vertex[j].y); - if (facet.vertex[j].z > -1e-12f && facet.vertex[j].z < 1e-12f) - printf("stl_read: facet %d.z = %e\r\n", j, facet.vertex[j].z); - } -#endif } else /* Read a single facet from an ASCII .STL file */ { @@ -319,6 +305,21 @@ stl_read(stl_file *stl, int first_facet, int first) { } } +#if 0 + // Report close to zero vertex coordinates. Due to the nature of the floating point numbers, + // close to zero values may be represented with singificantly higher precision than the rest of the vertices. + // It may be worth to round these numbers to zero during loading to reduce the number of errors reported + // during the STL import. + for (size_t j = 0; j < 3; ++ j) { + if (facet.vertex[j].x > -1e-12f && facet.vertex[j].x < 1e-12f) + printf("stl_read: facet %d.x = %e\r\n", j, facet.vertex[j].x); + if (facet.vertex[j].y > -1e-12f && facet.vertex[j].y < 1e-12f) + printf("stl_read: facet %d.y = %e\r\n", j, facet.vertex[j].y); + if (facet.vertex[j].z > -1e-12f && facet.vertex[j].z < 1e-12f) + printf("stl_read: facet %d.z = %e\r\n", j, facet.vertex[j].z); + } +#endif + #if 1 { // Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.