diff options
Diffstat (limited to 'src/admesh')
| -rw-r--r-- | src/admesh/CMakeLists.txt | 12 | ||||
| -rw-r--r-- | src/admesh/connect.cpp | 911 | ||||
| -rw-r--r-- | src/admesh/normals.cpp | 292 | ||||
| -rw-r--r-- | src/admesh/shared.cpp | 264 | ||||
| -rw-r--r-- | src/admesh/stl.h | 215 | ||||
| -rw-r--r-- | src/admesh/stl_io.cpp | 433 | ||||
| -rw-r--r-- | src/admesh/stlinit.cpp | 382 | ||||
| -rw-r--r-- | src/admesh/util.cpp | 518 |
8 files changed, 3027 insertions, 0 deletions
diff --git a/src/admesh/CMakeLists.txt b/src/admesh/CMakeLists.txt new file mode 100644 index 000000000..44c97c3f1 --- /dev/null +++ b/src/admesh/CMakeLists.txt @@ -0,0 +1,12 @@ +project(admesh) +cmake_minimum_required(VERSION 2.6) + +add_library(admesh STATIC + connect.cpp + normals.cpp + shared.cpp + stl.h + stl_io.cpp + stlinit.cpp + util.cpp +) diff --git a/src/admesh/connect.cpp b/src/admesh/connect.cpp new file mode 100644 index 000000000..da5b66720 --- /dev/null +++ b/src/admesh/connect.cpp @@ -0,0 +1,911 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include <boost/detail/endian.hpp> + +#include "stl.h" + + +static void stl_match_neighbors_nearby(stl_file *stl, + stl_hash_edge *edge_a, stl_hash_edge *edge_b); +static void stl_record_neighbors(stl_file *stl, + stl_hash_edge *edge_a, stl_hash_edge *edge_b); +static void stl_initialize_facet_check_exact(stl_file *stl); +static void stl_initialize_facet_check_nearby(stl_file *stl); +static void stl_load_edge_exact(stl_file *stl, stl_hash_edge *edge, + stl_vertex *a, stl_vertex *b); +static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge, + stl_vertex *a, stl_vertex *b, float tolerance); +static void insert_hash_edge(stl_file *stl, stl_hash_edge edge, + void (*match_neighbors)(stl_file *stl, + stl_hash_edge *edge_a, stl_hash_edge *edge_b)); +static int stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b); +static void stl_free_edges(stl_file *stl); +static void stl_remove_facet(stl_file *stl, int facet_number); +static void stl_change_vertices(stl_file *stl, int facet_num, int vnot, + stl_vertex new_vertex); +static void stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a, + stl_hash_edge *edge_b, int *facet1, int *vertex1, + int *facet2, int *vertex2, + stl_vertex *new_vertex1, stl_vertex *new_vertex2); +static void stl_remove_degenerate(stl_file *stl, int facet); +extern int stl_check_normal_vector(stl_file *stl, + int facet_num, int normal_fix_flag); +static void stl_update_connects_remove_1(stl_file *stl, int facet_num); + + +void +stl_check_facets_exact(stl_file *stl) { + /* This function builds the neighbors list. No modifications are made + * to any of the facets. The edges are said to match only if all six + * floats of the first edge matches all six floats of the second edge. + */ + + stl_hash_edge edge; + stl_facet facet; + int i; + int j; + + if (stl->error) return; + + stl->stats.connected_edges = 0; + stl->stats.connected_facets_1_edge = 0; + stl->stats.connected_facets_2_edge = 0; + stl->stats.connected_facets_3_edge = 0; + + stl_initialize_facet_check_exact(stl); + + for(i = 0; i < stl->stats.number_of_facets; i++) { + facet = stl->facet_start[i]; + // If any two of the three vertices are found to be exactally the same, call them degenerate and remove the facet. + if (facet.vertex[0] == facet.vertex[1] || + facet.vertex[1] == facet.vertex[2] || + facet.vertex[0] == facet.vertex[2]) { + stl->stats.degenerate_facets += 1; + stl_remove_facet(stl, i); + -- i; + continue; + } + for(j = 0; j < 3; j++) { + edge.facet_number = i; + edge.which_edge = j; + stl_load_edge_exact(stl, &edge, &facet.vertex[j], &facet.vertex[(j + 1) % 3]); + insert_hash_edge(stl, edge, stl_record_neighbors); + } + } + 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) { + + if (stl->error) return; + + { + stl_vertex diff = (*a - *b).cwiseAbs(); + float max_diff = std::max(diff(0), std::max(diff(1), diff(2))); + stl->stats.shortest_edge = std::min(max_diff, stl->stats.shortest_edge); + } + + // Ensure identical vertex ordering of equal edges. + // This method is numerically robust. + if (stl_vertex_lower(*a, *b)) { + } else { + std::swap(a, b); + edge->which_edge += 3; /* this edge is loaded backwards */ + } + memcpy(&edge->key[0], a->data(), sizeof(stl_vertex)); + memcpy(&edge->key[sizeof(stl_vertex)], b->data(), sizeof(stl_vertex)); + // Switch negative zeros to positive zeros, so memcmp will consider them to be equal. + for (size_t i = 0; i < 6; ++ i) { + unsigned char *p = edge->key + i * 4; +#ifdef BOOST_LITTLE_ENDIAN + if (p[0] == 0 && p[1] == 0 && p[2] == 0 && p[3] == 0x80) + // Negative zero, switch to positive zero. + p[3] = 0; +#else /* BOOST_LITTLE_ENDIAN */ + if (p[0] == 0x80 && p[1] == 0 && p[2] == 0 && p[3] == 0) + // Negative zero, switch to positive zero. + p[0] = 0; +#endif /* BOOST_LITTLE_ENDIAN */ + } +} + +static void +stl_initialize_facet_check_exact(stl_file *stl) { + int i; + + if (stl->error) return; + + stl->stats.malloced = 0; + stl->stats.freed = 0; + stl->stats.collisions = 0; + + + stl->M = 81397; + + for(i = 0; i < stl->stats.number_of_facets ; i++) { + /* initialize neighbors list to -1 to mark unconnected edges */ + stl->neighbors_start[i].neighbor[0] = -1; + stl->neighbors_start[i].neighbor[1] = -1; + stl->neighbors_start[i].neighbor[2] = -1; + } + + stl->heads = (stl_hash_edge**)calloc(stl->M, sizeof(*stl->heads)); + if(stl->heads == NULL) perror("stl_initialize_facet_check_exact"); + + stl->tail = (stl_hash_edge*)malloc(sizeof(stl_hash_edge)); + if(stl->tail == NULL) perror("stl_initialize_facet_check_exact"); + + stl->tail->next = stl->tail; + + for(i = 0; i < stl->M; i++) { + stl->heads[i] = stl->tail; + } +} + +static void insert_hash_edge(stl_file *stl, stl_hash_edge edge, + void (*match_neighbors)(stl_file *stl, + stl_hash_edge *edge_a, stl_hash_edge *edge_b)) +{ + if (stl->error) return; + + int chain_number = edge.hash(stl->M); + stl_hash_edge *link = stl->heads[chain_number]; + + stl_hash_edge *new_edge; + stl_hash_edge *temp; + if(link == stl->tail) { + /* This list doesn't have any edges currently in it. Add this one. */ + new_edge = (stl_hash_edge*)malloc(sizeof(stl_hash_edge)); + if(new_edge == NULL) perror("insert_hash_edge"); + stl->stats.malloced++; + *new_edge = edge; + new_edge->next = stl->tail; + stl->heads[chain_number] = new_edge; + return; + } else if(!stl_compare_function(&edge, link)) { + /* This is a match. Record result in neighbors list. */ + match_neighbors(stl, &edge, link); + /* Delete the matched edge from the list. */ + stl->heads[chain_number] = link->next; + free(link); + stl->stats.freed++; + return; + } else { + /* Continue through the rest of the list */ + for(;;) { + if(link->next == stl->tail) { + /* This is the last item in the list. Insert a new edge. */ + new_edge = (stl_hash_edge*)malloc(sizeof(stl_hash_edge)); + if(new_edge == NULL) perror("insert_hash_edge"); + stl->stats.malloced++; + *new_edge = edge; + new_edge->next = stl->tail; + link->next = new_edge; + stl->stats.collisions++; + return; + } else if(!stl_compare_function(&edge, link->next)) { + /* This is a match. Record result in neighbors list. */ + match_neighbors(stl, &edge, link->next); + + /* Delete the matched edge from the list. */ + temp = link->next; + link->next = link->next->next; + free(temp); + stl->stats.freed++; + return; + } else { + /* This is not a match. Go to the next link */ + link = link->next; + stl->stats.collisions++; + } + } + } +} + +// Return 1 if the edges are not matched. +static inline int stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b) +{ + // Don't match edges of the same facet + return (edge_a->facet_number == edge_b->facet_number) || (*edge_a != *edge_b); +} + +void stl_check_facets_nearby(stl_file *stl, float tolerance) +{ + if (stl->error) + return; + + if( (stl->stats.connected_facets_1_edge == stl->stats.number_of_facets) + && (stl->stats.connected_facets_2_edge == stl->stats.number_of_facets) + && (stl->stats.connected_facets_3_edge == stl->stats.number_of_facets)) { + /* No need to check any further. All facets are connected */ + return; + } + + stl_initialize_facet_check_nearby(stl); + + for (int i = 0; i < stl->stats.number_of_facets; ++ i) { + //FIXME is the copy necessary? + stl_facet facet = stl->facet_start[i]; + for (int j = 0; j < 3; j++) { + if(stl->neighbors_start[i].neighbor[j] == -1) { + stl_hash_edge edge; + edge.facet_number = i; + edge.which_edge = j; + if(stl_load_edge_nearby(stl, &edge, &facet.vertex[j], + &facet.vertex[(j + 1) % 3], + tolerance)) { + /* only insert edges that have different keys */ + insert_hash_edge(stl, edge, stl_match_neighbors_nearby); + } + } + } + } + + stl_free_edges(stl); +} + +static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge, stl_vertex *a, stl_vertex *b, float tolerance) +{ + // Index of a grid cell spaced by tolerance. + typedef Eigen::Matrix<int32_t, 3, 1, Eigen::DontAlign> Vec3i; + Vec3i vertex1 = (*a / tolerance).cast<int32_t>(); + Vec3i vertex2 = (*b / tolerance).cast<int32_t>(); + static_assert(sizeof(Vec3i) == 12, "size of Vec3i incorrect"); + + if (vertex1 == vertex2) + // Both vertices hash to the same value + return 0; + + // 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.data(), sizeof(stl_vertex)); + memcpy(&edge->key[sizeof(stl_vertex)], vertex2.data(), sizeof(stl_vertex)); + } else { + memcpy(&edge->key[0], vertex2.data(), sizeof(stl_vertex)); + memcpy(&edge->key[sizeof(stl_vertex)], vertex1.data(), sizeof(stl_vertex)); + edge->which_edge += 3; /* this edge is loaded backwards */ + } + return 1; +} + +static void stl_free_edges(stl_file *stl) +{ + if (stl->error) + return; + + if(stl->stats.malloced != stl->stats.freed) { + for (int i = 0; i < stl->M; i++) { + for (stl_hash_edge *temp = stl->heads[i]; stl->heads[i] != stl->tail; temp = stl->heads[i]) { + stl->heads[i] = stl->heads[i]->next; + free(temp); + ++ stl->stats.freed; + } + } + } + free(stl->heads); + free(stl->tail); +} + +static void stl_initialize_facet_check_nearby(stl_file *stl) +{ + int i; + + if (stl->error) return; + + stl->stats.malloced = 0; + stl->stats.freed = 0; + stl->stats.collisions = 0; + + /* tolerance = STL_MAX(stl->stats.shortest_edge, tolerance);*/ + /* tolerance = STL_MAX((stl->stats.bounding_diameter / 500000.0), tolerance);*/ + /* tolerance *= 0.5;*/ + + stl->M = 81397; + + stl->heads = (stl_hash_edge**)calloc(stl->M, sizeof(*stl->heads)); + if(stl->heads == NULL) perror("stl_initialize_facet_check_nearby"); + + stl->tail = (stl_hash_edge*)malloc(sizeof(stl_hash_edge)); + if(stl->tail == NULL) perror("stl_initialize_facet_check_nearby"); + + stl->tail->next = stl->tail; + + for(i = 0; i < stl->M; i++) { + stl->heads[i] = stl->tail; + } +} + + + +static void +stl_record_neighbors(stl_file *stl, + stl_hash_edge *edge_a, stl_hash_edge *edge_b) { + int i; + int j; + + if (stl->error) return; + + /* Facet a's neighbor is facet b */ + stl->neighbors_start[edge_a->facet_number].neighbor[edge_a->which_edge % 3] = + edge_b->facet_number; /* sets the .neighbor part */ + + stl->neighbors_start[edge_a->facet_number]. + which_vertex_not[edge_a->which_edge % 3] = + (edge_b->which_edge + 2) % 3; /* sets the .which_vertex_not part */ + + /* Facet b's neighbor is facet a */ + stl->neighbors_start[edge_b->facet_number].neighbor[edge_b->which_edge % 3] = + edge_a->facet_number; /* sets the .neighbor part */ + + stl->neighbors_start[edge_b->facet_number]. + which_vertex_not[edge_b->which_edge % 3] = + (edge_a->which_edge + 2) % 3; /* sets the .which_vertex_not part */ + + if( ((edge_a->which_edge < 3) && (edge_b->which_edge < 3)) + || ((edge_a->which_edge > 2) && (edge_b->which_edge > 2))) { + /* these facets are oriented in opposite directions. */ + /* their normals are probably messed up. */ + stl->neighbors_start[edge_a->facet_number]. + which_vertex_not[edge_a->which_edge % 3] += 3; + stl->neighbors_start[edge_b->facet_number]. + which_vertex_not[edge_b->which_edge % 3] += 3; + } + + + /* Count successful connects */ + /* Total connects */ + stl->stats.connected_edges += 2; + /* Count individual connects */ + i = ((stl->neighbors_start[edge_a->facet_number].neighbor[0] == -1) + + (stl->neighbors_start[edge_a->facet_number].neighbor[1] == -1) + + (stl->neighbors_start[edge_a->facet_number].neighbor[2] == -1)); + j = ((stl->neighbors_start[edge_b->facet_number].neighbor[0] == -1) + + (stl->neighbors_start[edge_b->facet_number].neighbor[1] == -1) + + (stl->neighbors_start[edge_b->facet_number].neighbor[2] == -1)); + if(i == 2) { + stl->stats.connected_facets_1_edge +=1; + } else if(i == 1) { + stl->stats.connected_facets_2_edge +=1; + } else { + stl->stats.connected_facets_3_edge +=1; + } + if(j == 2) { + stl->stats.connected_facets_1_edge +=1; + } else if(j == 1) { + stl->stats.connected_facets_2_edge +=1; + } else { + stl->stats.connected_facets_3_edge +=1; + } +} + +static void stl_match_neighbors_nearby(stl_file *stl, stl_hash_edge *edge_a, stl_hash_edge *edge_b) +{ + int facet1; + int facet2; + int vertex1; + int vertex2; + int vnot1; + int vnot2; + stl_vertex new_vertex1; + stl_vertex new_vertex2; + + if (stl->error) return; + + stl_record_neighbors(stl, edge_a, edge_b); + stl_which_vertices_to_change(stl, edge_a, edge_b, &facet1, &vertex1, + &facet2, &vertex2, &new_vertex1, &new_vertex2); + if(facet1 != -1) { + if(facet1 == edge_a->facet_number) { + vnot1 = (edge_a->which_edge + 2) % 3; + } else { + vnot1 = (edge_b->which_edge + 2) % 3; + } + if(((vnot1 + 2) % 3) == vertex1) { + vnot1 += 3; + } + stl_change_vertices(stl, facet1, vnot1, new_vertex1); + } + if(facet2 != -1) { + if(facet2 == edge_a->facet_number) { + vnot2 = (edge_a->which_edge + 2) % 3; + } else { + vnot2 = (edge_b->which_edge + 2) % 3; + } + if(((vnot2 + 2) % 3) == vertex2) { + vnot2 += 3; + } + stl_change_vertices(stl, facet2, vnot2, new_vertex2); + } + stl->stats.edges_fixed += 2; +} + + +static void stl_change_vertices(stl_file *stl, int facet_num, int vnot, stl_vertex new_vertex) { + int first_facet; + int direction; + int next_edge; + int pivot_vertex; + + if (stl->error) return; + + first_facet = facet_num; + direction = 0; + + for(;;) { + if(vnot > 2) { + if(direction == 0) { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + direction = 1; + } else { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot % 3; + direction = 0; + } + } else { + if(direction == 0) { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot; + } else { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + } + } +#if 0 + if (stl->facet_start[facet_num].vertex[pivot_vertex](0) == new_vertex(0) && + stl->facet_start[facet_num].vertex[pivot_vertex](1) == new_vertex(1) && + stl->facet_start[facet_num].vertex[pivot_vertex](2) == new_vertex(2)) + printf("Changing vertex %f,%f,%f: Same !!!\r\n", + new_vertex(0), new_vertex(1), new_vertex(2)); + else { + if (stl->facet_start[facet_num].vertex[pivot_vertex](0) != new_vertex(0)) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex](0), + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](0)), + new_vertex(0), + *reinterpret_cast<const int*>(&new_vertex(0))); + if (stl->facet_start[facet_num].vertex[pivot_vertex](1) != new_vertex(1)) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex](1), + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](1)), + new_vertex(1), + *reinterpret_cast<const int*>(&new_vertex(1))); + if (stl->facet_start[facet_num].vertex[pivot_vertex](2) != new_vertex(2)) + printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", + stl->facet_start[facet_num].vertex[pivot_vertex](2), + *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](2)), + new_vertex(2), + *reinterpret_cast<const int*>(&new_vertex(2))); + } +#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]; + + if(facet_num == -1) { + break; + } + + if(facet_num == first_facet) { + /* back to the beginning */ + printf("\ +Back to the first facet changing vertices: probably a mobius part.\n\ +Try using a smaller tolerance or don't do a nearby check\n"); + return; + } + } +} + +static void +stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a, + stl_hash_edge *edge_b, int *facet1, int *vertex1, + int *facet2, int *vertex2, + stl_vertex *new_vertex1, stl_vertex *new_vertex2) { + int v1a; /* pair 1, facet a */ + int v1b; /* pair 1, facet b */ + int v2a; /* pair 2, facet a */ + int v2b; /* pair 2, facet b */ + + /* Find first pair */ + if(edge_a->which_edge < 3) { + v1a = edge_a->which_edge; + v2a = (edge_a->which_edge + 1) % 3; + } else { + v2a = edge_a->which_edge % 3; + v1a = (edge_a->which_edge + 1) % 3; + } + if(edge_b->which_edge < 3) { + v1b = edge_b->which_edge; + v2b = (edge_b->which_edge + 1) % 3; + } else { + v2b = edge_b->which_edge % 3; + v1b = (edge_b->which_edge + 1) % 3; + } + + // Of the first pair, which vertex, if any, should be changed + if(stl->facet_start[edge_a->facet_number].vertex[v1a] == + stl->facet_start[edge_b->facet_number].vertex[v1b]) { + // These facets are already equal. No need to change. + *facet1 = -1; + } else { + if( (stl->neighbors_start[edge_a->facet_number].neighbor[v1a] == -1) + && (stl->neighbors_start[edge_a->facet_number]. + neighbor[(v1a + 2) % 3] == -1)) { + /* This vertex has no neighbors. This is a good one to change */ + *facet1 = edge_a->facet_number; + *vertex1 = v1a; + *new_vertex1 = stl->facet_start[edge_b->facet_number].vertex[v1b]; + } else { + *facet1 = edge_b->facet_number; + *vertex1 = v1b; + *new_vertex1 = stl->facet_start[edge_a->facet_number].vertex[v1a]; + } + } + + /* Of the second pair, which vertex, if any, should be changed */ + if(stl->facet_start[edge_a->facet_number].vertex[v2a] == + stl->facet_start[edge_b->facet_number].vertex[v2b]) { + // These facets are already equal. No need to change. + *facet2 = -1; + } else { + if( (stl->neighbors_start[edge_a->facet_number].neighbor[v2a] == -1) + && (stl->neighbors_start[edge_a->facet_number]. + neighbor[(v2a + 2) % 3] == -1)) { + /* This vertex has no neighbors. This is a good one to change */ + *facet2 = edge_a->facet_number; + *vertex2 = v2a; + *new_vertex2 = stl->facet_start[edge_b->facet_number].vertex[v2b]; + } else { + *facet2 = edge_b->facet_number; + *vertex2 = v2b; + *new_vertex2 = stl->facet_start[edge_a->facet_number].vertex[v2a]; + } + } +} + +static void +stl_remove_facet(stl_file *stl, int facet_number) { + int neighbor[3]; + int vnot[3]; + int i; + int j; + + if (stl->error) return; + + stl->stats.facets_removed += 1; + /* Update list of connected edges */ + j = ((stl->neighbors_start[facet_number].neighbor[0] == -1) + + (stl->neighbors_start[facet_number].neighbor[1] == -1) + + (stl->neighbors_start[facet_number].neighbor[2] == -1)); + if(j == 2) { + stl->stats.connected_facets_1_edge -= 1; + } else if(j == 1) { + stl->stats.connected_facets_2_edge -= 1; + stl->stats.connected_facets_1_edge -= 1; + } else if(j == 0) { + stl->stats.connected_facets_3_edge -= 1; + stl->stats.connected_facets_2_edge -= 1; + stl->stats.connected_facets_1_edge -= 1; + } + + stl->facet_start[facet_number] = + stl->facet_start[stl->stats.number_of_facets - 1]; + /* I could reallocate at this point, but it is not really necessary. */ + stl->neighbors_start[facet_number] = + stl->neighbors_start[stl->stats.number_of_facets - 1]; + stl->stats.number_of_facets -= 1; + + for(i = 0; i < 3; i++) { + neighbor[i] = stl->neighbors_start[facet_number].neighbor[i]; + vnot[i] = stl->neighbors_start[facet_number].which_vertex_not[i]; + } + + for(i = 0; i < 3; i++) { + if(neighbor[i] != -1) { + if(stl->neighbors_start[neighbor[i]].neighbor[(vnot[i] + 1)% 3] != + stl->stats.number_of_facets) { + printf("\ +in stl_remove_facet: neighbor = %d numfacets = %d this is wrong\n", + stl->neighbors_start[neighbor[i]].neighbor[(vnot[i] + 1)% 3], + stl->stats.number_of_facets); + return; + } + stl->neighbors_start[neighbor[i]].neighbor[(vnot[i] + 1)% 3] + = facet_number; + } + } +} + +void stl_remove_unconnected_facets(stl_file *stl) +{ + /* A couple of things need to be done here. One is to remove any */ + /* completely unconnected facets (0 edges connected) since these are */ + /* useless and could be completely wrong. The second thing that needs to */ + /* be done is to remove any degenerate facets that were created during */ + /* stl_check_facets_nearby(). */ + if (stl->error) + return; + + // remove degenerate facets + for (int i = 0; i < stl->stats.number_of_facets; ++ i) { + if(stl->facet_start[i].vertex[0] == stl->facet_start[i].vertex[1] || + stl->facet_start[i].vertex[0] == stl->facet_start[i].vertex[2] || + stl->facet_start[i].vertex[1] == stl->facet_start[i].vertex[2]) { + stl_remove_degenerate(stl, i); + i--; + } + } + + if(stl->stats.connected_facets_1_edge < stl->stats.number_of_facets) { + // remove completely unconnected facets + for (int i = 0; i < stl->stats.number_of_facets; i++) { + if (stl->neighbors_start[i].neighbor[0] == -1 && + stl->neighbors_start[i].neighbor[1] == -1 && + stl->neighbors_start[i].neighbor[2] == -1) { + // This facet is completely unconnected. Remove it. + stl_remove_facet(stl, i); + -- i; + } + } + } +} + +static void +stl_remove_degenerate(stl_file *stl, int facet) { + int edge1; + int edge2; + int edge3; + int neighbor1; + int neighbor2; + int neighbor3; + int vnot1; + int vnot2; + int vnot3; + + if (stl->error) return; + + if (stl->facet_start[facet].vertex[0] == stl->facet_start[facet].vertex[1] && + stl->facet_start[facet].vertex[1] == stl->facet_start[facet].vertex[2]) { + /* all 3 vertices are equal. Just remove the facet. I don't think*/ + /* this is really possible, but just in case... */ + printf("removing a facet in stl_remove_degenerate\n"); + stl_remove_facet(stl, facet); + return; + } + + if (stl->facet_start[facet].vertex[0] == stl->facet_start[facet].vertex[1]) { + edge1 = 1; + edge2 = 2; + edge3 = 0; + } else if (stl->facet_start[facet].vertex[1] == stl->facet_start[facet].vertex[2]) { + edge1 = 0; + edge2 = 2; + edge3 = 1; + } else if (stl->facet_start[facet].vertex[2] == stl->facet_start[facet].vertex[0]) { + edge1 = 0; + edge2 = 1; + edge3 = 2; + } else { + /* No degenerate. Function shouldn't have been called. */ + return; + } + neighbor1 = stl->neighbors_start[facet].neighbor[edge1]; + neighbor2 = stl->neighbors_start[facet].neighbor[edge2]; + + if(neighbor1 == -1) { + stl_update_connects_remove_1(stl, neighbor2); + } + if(neighbor2 == -1) { + stl_update_connects_remove_1(stl, neighbor1); + } + + + neighbor3 = stl->neighbors_start[facet].neighbor[edge3]; + vnot1 = stl->neighbors_start[facet].which_vertex_not[edge1]; + vnot2 = stl->neighbors_start[facet].which_vertex_not[edge2]; + vnot3 = stl->neighbors_start[facet].which_vertex_not[edge3]; + + if(neighbor1 >= 0){ + stl->neighbors_start[neighbor1].neighbor[(vnot1 + 1) % 3] = neighbor2; + stl->neighbors_start[neighbor1].which_vertex_not[(vnot1 + 1) % 3] = vnot2; + } + if(neighbor2 >= 0){ + stl->neighbors_start[neighbor2].neighbor[(vnot2 + 1) % 3] = neighbor1; + stl->neighbors_start[neighbor2].which_vertex_not[(vnot2 + 1) % 3] = vnot1; + } + + stl_remove_facet(stl, facet); + + if(neighbor3 >= 0) { + stl_update_connects_remove_1(stl, neighbor3); + stl->neighbors_start[neighbor3].neighbor[(vnot3 + 1) % 3] = -1; + } +} + +void +stl_update_connects_remove_1(stl_file *stl, int facet_num) { + int j; + + if (stl->error) return; + /* Update list of connected edges */ + j = ((stl->neighbors_start[facet_num].neighbor[0] == -1) + + (stl->neighbors_start[facet_num].neighbor[1] == -1) + + (stl->neighbors_start[facet_num].neighbor[2] == -1)); + if(j == 0) { /* Facet has 3 neighbors */ + stl->stats.connected_facets_3_edge -= 1; + } else if(j == 1) { /* Facet has 2 neighbors */ + stl->stats.connected_facets_2_edge -= 1; + } else if(j == 2) { /* Facet has 1 neighbor */ + stl->stats.connected_facets_1_edge -= 1; + } +} + +void +stl_fill_holes(stl_file *stl) { + stl_facet facet; + stl_facet new_facet; + int neighbors_initial[3]; + stl_hash_edge edge; + int first_facet; + int direction; + int facet_num; + int vnot; + int next_edge; + int pivot_vertex; + int next_facet; + int i; + int j; + int k; + + if (stl->error) return; + + /* Insert all unconnected edges into hash list */ + stl_initialize_facet_check_nearby(stl); + for(i = 0; i < stl->stats.number_of_facets; i++) { + facet = stl->facet_start[i]; + for(j = 0; j < 3; j++) { + if(stl->neighbors_start[i].neighbor[j] != -1) continue; + edge.facet_number = i; + edge.which_edge = j; + stl_load_edge_exact(stl, &edge, &facet.vertex[j], + &facet.vertex[(j + 1) % 3]); + + insert_hash_edge(stl, edge, stl_record_neighbors); + } + } + + for(i = 0; i < stl->stats.number_of_facets; i++) { + facet = stl->facet_start[i]; + neighbors_initial[0] = stl->neighbors_start[i].neighbor[0]; + neighbors_initial[1] = stl->neighbors_start[i].neighbor[1]; + neighbors_initial[2] = stl->neighbors_start[i].neighbor[2]; + first_facet = i; + for(j = 0; j < 3; j++) { + if(stl->neighbors_start[i].neighbor[j] != -1) continue; + + new_facet.vertex[0] = facet.vertex[j]; + new_facet.vertex[1] = facet.vertex[(j + 1) % 3]; + if(neighbors_initial[(j + 2) % 3] == -1) { + direction = 1; + } else { + direction = 0; + } + + facet_num = i; + vnot = (j + 2) % 3; + + for(;;) { + if(vnot > 2) { + if(direction == 0) { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + direction = 1; + } else { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot % 3; + direction = 0; + } + } else { + if(direction == 0) { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot; + } else { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + } + } + next_facet = stl->neighbors_start[facet_num].neighbor[next_edge]; + + if(next_facet == -1) { + new_facet.vertex[2] = stl->facet_start[facet_num]. + vertex[vnot % 3]; + stl_add_facet(stl, &new_facet); + for(k = 0; k < 3; k++) { + edge.facet_number = stl->stats.number_of_facets - 1; + edge.which_edge = k; + stl_load_edge_exact(stl, &edge, &new_facet.vertex[k], + &new_facet.vertex[(k + 1) % 3]); + + insert_hash_edge(stl, edge, stl_record_neighbors); + } + break; + } else { + vnot = stl->neighbors_start[facet_num]. + which_vertex_not[next_edge]; + facet_num = next_facet; + } + + if(facet_num == first_facet) { + /* back to the beginning */ + printf("\ +Back to the first facet filling holes: probably a mobius part.\n\ +Try using a smaller tolerance or don't do a nearby check\n"); + return; + } + } + } + } +} + +void +stl_add_facet(stl_file *stl, stl_facet *new_facet) { + if (stl->error) return; + + stl->stats.facets_added += 1; + if(stl->stats.facets_malloced < stl->stats.number_of_facets + 1) { + stl->facet_start = (stl_facet*)realloc(stl->facet_start, + (sizeof(stl_facet) * (stl->stats.facets_malloced + 256))); + if(stl->facet_start == NULL) perror("stl_add_facet"); + stl->neighbors_start = (stl_neighbors*)realloc(stl->neighbors_start, + (sizeof(stl_neighbors) * (stl->stats.facets_malloced + 256))); + if(stl->neighbors_start == NULL) perror("stl_add_facet"); + stl->stats.facets_malloced += 256; + } + stl->facet_start[stl->stats.number_of_facets] = *new_facet; + + /* note that the normal vector is not set here, just initialized to 0 */ + stl->facet_start[stl->stats.number_of_facets].normal = stl_normal::Zero(); + + stl->neighbors_start[stl->stats.number_of_facets].neighbor[0] = -1; + stl->neighbors_start[stl->stats.number_of_facets].neighbor[1] = -1; + stl->neighbors_start[stl->stats.number_of_facets].neighbor[2] = -1; + stl->stats.number_of_facets += 1; +} diff --git a/src/admesh/normals.cpp b/src/admesh/normals.cpp new file mode 100644 index 000000000..a8faa44bd --- /dev/null +++ b/src/admesh/normals.cpp @@ -0,0 +1,292 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "stl.h" + +static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag); + +static void +stl_reverse_facet(stl_file *stl, int facet_num) { + stl_vertex tmp_vertex; + /* int tmp_neighbor;*/ + int neighbor[3]; + int vnot[3]; + + stl->stats.facets_reversed += 1; + + neighbor[0] = stl->neighbors_start[facet_num].neighbor[0]; + neighbor[1] = stl->neighbors_start[facet_num].neighbor[1]; + neighbor[2] = stl->neighbors_start[facet_num].neighbor[2]; + vnot[0] = stl->neighbors_start[facet_num].which_vertex_not[0]; + vnot[1] = stl->neighbors_start[facet_num].which_vertex_not[1]; + vnot[2] = stl->neighbors_start[facet_num].which_vertex_not[2]; + + /* reverse the facet */ + tmp_vertex = stl->facet_start[facet_num].vertex[0]; + stl->facet_start[facet_num].vertex[0] = + stl->facet_start[facet_num].vertex[1]; + stl->facet_start[facet_num].vertex[1] = tmp_vertex; + + /* fix the vnots of the neighboring facets */ + if(neighbor[0] != -1) + stl->neighbors_start[neighbor[0]].which_vertex_not[(vnot[0] + 1) % 3] = + (stl->neighbors_start[neighbor[0]]. + which_vertex_not[(vnot[0] + 1) % 3] + 3) % 6; + if(neighbor[1] != -1) + stl->neighbors_start[neighbor[1]].which_vertex_not[(vnot[1] + 1) % 3] = + (stl->neighbors_start[neighbor[1]]. + which_vertex_not[(vnot[1] + 1) % 3] + 4) % 6; + if(neighbor[2] != -1) + stl->neighbors_start[neighbor[2]].which_vertex_not[(vnot[2] + 1) % 3] = + (stl->neighbors_start[neighbor[2]]. + which_vertex_not[(vnot[2] + 1) % 3] + 2) % 6; + + /* swap the neighbors of the facet that is being reversed */ + stl->neighbors_start[facet_num].neighbor[1] = neighbor[2]; + stl->neighbors_start[facet_num].neighbor[2] = neighbor[1]; + + /* swap the vnots of the facet that is being reversed */ + stl->neighbors_start[facet_num].which_vertex_not[1] = vnot[2]; + stl->neighbors_start[facet_num].which_vertex_not[2] = vnot[1]; + + /* reverse the values of the vnots of the facet that is being reversed */ + stl->neighbors_start[facet_num].which_vertex_not[0] = + (stl->neighbors_start[facet_num].which_vertex_not[0] + 3) % 6; + stl->neighbors_start[facet_num].which_vertex_not[1] = + (stl->neighbors_start[facet_num].which_vertex_not[1] + 3) % 6; + stl->neighbors_start[facet_num].which_vertex_not[2] = + (stl->neighbors_start[facet_num].which_vertex_not[2] + 3) % 6; +} + +void +stl_fix_normal_directions(stl_file *stl) { + char *norm_sw; + /* int edge_num;*/ + /* int vnot;*/ + int checked = 0; + int facet_num; + /* int next_facet;*/ + int i; + int j; + struct stl_normal { + int facet_num; + struct stl_normal *next; + }; + struct stl_normal *head; + struct stl_normal *tail; + struct stl_normal *newn; + struct stl_normal *temp; + + int* reversed_ids; + int reversed_count = 0; + int id; + int force_exit = 0; + + if (stl->error) return; + + /* Initialize linked list. */ + head = (struct stl_normal*)malloc(sizeof(struct stl_normal)); + if(head == NULL) perror("stl_fix_normal_directions"); + tail = (struct stl_normal*)malloc(sizeof(struct stl_normal)); + if(tail == NULL) perror("stl_fix_normal_directions"); + head->next = tail; + tail->next = tail; + + /* Initialize list that keeps track of already fixed facets. */ + norm_sw = (char*)calloc(stl->stats.number_of_facets, sizeof(char)); + if(norm_sw == NULL) perror("stl_fix_normal_directions"); + + /* Initialize list that keeps track of reversed facets. */ + reversed_ids = (int*)calloc(stl->stats.number_of_facets, sizeof(int)); + if (reversed_ids == NULL) perror("stl_fix_normal_directions reversed_ids"); + + facet_num = 0; + /* If normal vector is not within tolerance and backwards: + Arbitrarily starts at face 0. If this one is wrong, we're screwed. Thankfully, the chances + of it being wrong randomly are low if most of the triangles are right: */ + if (stl_check_normal_vector(stl, 0, 0) == 2) { + stl_reverse_facet(stl, 0); + reversed_ids[reversed_count++] = 0; + } + + /* Say that we've fixed this facet: */ + norm_sw[facet_num] = 1; + checked++; + + for(;;) { + /* Add neighbors_to_list. + Add unconnected neighbors to the list:a */ + for(j = 0; j < 3; j++) { + /* Reverse the neighboring facets if necessary. */ + if(stl->neighbors_start[facet_num].which_vertex_not[j] > 2) { + /* If the facet has a neighbor that is -1, it means that edge isn't shared by another facet */ + if(stl->neighbors_start[facet_num].neighbor[j] != -1) { + if (norm_sw[stl->neighbors_start[facet_num].neighbor[j]] == 1) { + /* trying to modify a facet already marked as fixed, revert all changes made until now and exit (fixes: #716, #574, #413, #269, #262, #259, #230, #228, #206) */ + for (id = reversed_count - 1; id >= 0; --id) { + stl_reverse_facet(stl, reversed_ids[id]); + } + force_exit = 1; + break; + } else { + stl_reverse_facet(stl, stl->neighbors_start[facet_num].neighbor[j]); + reversed_ids[reversed_count++] = stl->neighbors_start[facet_num].neighbor[j]; + } + } + } + /* If this edge of the facet is connected: */ + if(stl->neighbors_start[facet_num].neighbor[j] != -1) { + /* If we haven't fixed this facet yet, add it to the list: */ + if(norm_sw[stl->neighbors_start[facet_num].neighbor[j]] != 1) { + /* Add node to beginning of list. */ + newn = (struct stl_normal*)malloc(sizeof(struct stl_normal)); + if(newn == NULL) perror("stl_fix_normal_directions"); + newn->facet_num = stl->neighbors_start[facet_num].neighbor[j]; + newn->next = head->next; + head->next = newn; + } + } + } + + /* an error occourred, quit the for loop and exit */ + if (force_exit) break; + + /* Get next facet to fix from top of list. */ + if(head->next != tail) { + facet_num = head->next->facet_num; + if(norm_sw[facet_num] != 1) { /* If facet is in list mutiple times */ + norm_sw[facet_num] = 1; /* Record this one as being fixed. */ + checked++; + } + temp = head->next; /* Delete this facet from the list. */ + head->next = head->next->next; + free(temp); + } else { /* if we ran out of facets to fix: */ + /* All of the facets in this part have been fixed. */ + stl->stats.number_of_parts += 1; + if(checked >= stl->stats.number_of_facets) { + /* All of the facets have been checked. Bail out. */ + break; + } else { + /* There is another part here. Find it and continue. */ + for(i = 0; i < stl->stats.number_of_facets; i++) { + if(norm_sw[i] == 0) { + /* This is the first facet of the next part. */ + facet_num = i; + if(stl_check_normal_vector(stl, i, 0) == 2) { + stl_reverse_facet(stl, i); + reversed_ids[reversed_count++] = i; + } + + norm_sw[facet_num] = 1; + checked++; + break; + } + } + } + } + } + free(head); + free(tail); + free(reversed_ids); + free(norm_sw); +} + +static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag) { + /* Returns 0 if the normal is within tolerance */ + /* Returns 1 if the normal is not within tolerance, but direction is OK */ + /* Returns 2 if the normal is not within tolerance and backwards */ + /* Returns 4 if the status is unknown. */ + + stl_facet *facet; + + facet = &stl->facet_start[facet_num]; + + stl_normal normal; + stl_calculate_normal(normal, facet); + stl_normalize_vector(normal); + stl_normal normal_dif = (normal - facet->normal).cwiseAbs(); + + const float eps = 0.001f; + if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) { + /* It is not really necessary to change the values here */ + /* but just for consistency, I will. */ + facet->normal = normal; + return 0; + } + + stl_normal test_norm = facet->normal; + stl_normalize_vector(test_norm); + normal_dif = (normal - test_norm).cwiseAbs(); + if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) { + if(normal_fix_flag) { + facet->normal = normal; + stl->stats.normals_fixed += 1; + } + return 1; + } + + test_norm *= -1.f; + normal_dif = (normal - test_norm).cwiseAbs(); + if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) { + // Facet is backwards. + if(normal_fix_flag) { + facet->normal = normal; + stl->stats.normals_fixed += 1; + } + return 2; + } + if(normal_fix_flag) { + facet->normal = normal; + stl->stats.normals_fixed += 1; + } + return 4; +} + +void stl_fix_normal_values(stl_file *stl) { + int i; + + if (stl->error) return; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + stl_check_normal_vector(stl, i, 1); + } +} + +void stl_reverse_all_facets(stl_file *stl) +{ + if (stl->error) + return; + + stl_normal normal; + for(int i = 0; i < stl->stats.number_of_facets; i++) { + stl_reverse_facet(stl, i); + stl_calculate_normal(normal, &stl->facet_start[i]); + stl_normalize_vector(normal); + stl->facet_start[i].normal = normal; + } +} diff --git a/src/admesh/shared.cpp b/src/admesh/shared.cpp new file mode 100644 index 000000000..91bb82e00 --- /dev/null +++ b/src/admesh/shared.cpp @@ -0,0 +1,264 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdlib.h> +#include <string.h> + +#include <boost/nowide/cstdio.hpp> + +#include "stl.h" + +void +stl_invalidate_shared_vertices(stl_file *stl) { + if (stl->error) return; + + if (stl->v_indices != NULL) { + free(stl->v_indices); + stl->v_indices = NULL; + } + if (stl->v_shared != NULL) { + free(stl->v_shared); + stl->v_shared = NULL; + } +} + +void +stl_generate_shared_vertices(stl_file *stl) { + int i; + int j; + int first_facet; + int direction; + int facet_num; + int vnot; + int next_edge; + int pivot_vertex; + int next_facet; + int reversed; + + if (stl->error) return; + + /* make sure this function is idempotent and does not leak memory */ + stl_invalidate_shared_vertices(stl); + + stl->v_indices = (v_indices_struct*) + calloc(stl->stats.number_of_facets, sizeof(v_indices_struct)); + if(stl->v_indices == NULL) perror("stl_generate_shared_vertices"); + stl->v_shared = (stl_vertex*) + calloc((stl->stats.number_of_facets / 2), sizeof(stl_vertex)); + if(stl->v_shared == NULL) perror("stl_generate_shared_vertices"); + stl->stats.shared_malloced = stl->stats.number_of_facets / 2; + stl->stats.shared_vertices = 0; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + stl->v_indices[i].vertex[0] = -1; + stl->v_indices[i].vertex[1] = -1; + stl->v_indices[i].vertex[2] = -1; + } + + + for(i = 0; i < stl->stats.number_of_facets; i++) { + first_facet = i; + for(j = 0; j < 3; j++) { + if(stl->v_indices[i].vertex[j] != -1) { + continue; + } + if(stl->stats.shared_vertices == stl->stats.shared_malloced) { + stl->stats.shared_malloced += 1024; + stl->v_shared = (stl_vertex*)realloc(stl->v_shared, + stl->stats.shared_malloced * sizeof(stl_vertex)); + if(stl->v_shared == NULL) perror("stl_generate_shared_vertices"); + } + + stl->v_shared[stl->stats.shared_vertices] = + stl->facet_start[i].vertex[j]; + + direction = 0; + reversed = 0; + facet_num = i; + vnot = (j + 2) % 3; + + for(;;) { + if(vnot > 2) { + if(direction == 0) { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + direction = 1; + } else { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot % 3; + direction = 0; + } + } else { + if(direction == 0) { + pivot_vertex = (vnot + 1) % 3; + next_edge = vnot; + } else { + pivot_vertex = (vnot + 2) % 3; + next_edge = pivot_vertex; + } + } + stl->v_indices[facet_num].vertex[pivot_vertex] = + stl->stats.shared_vertices; + + next_facet = stl->neighbors_start[facet_num].neighbor[next_edge]; + if(next_facet == -1) { + if(reversed) { + break; + } else { + direction = 1; + vnot = (j + 1) % 3; + reversed = 1; + facet_num = first_facet; + } + } else if(next_facet != first_facet) { + vnot = stl->neighbors_start[facet_num]. + which_vertex_not[next_edge]; + facet_num = next_facet; + } else { + break; + } + } + stl->stats.shared_vertices += 1; + } + } +} + +void +stl_write_off(stl_file *stl, char *file) { + int i; + FILE *fp; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_ascii: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "OFF\n"); + fprintf(fp, "%d %d 0\n", + stl->stats.shared_vertices, stl->stats.number_of_facets); + + for(i = 0; i < stl->stats.shared_vertices; i++) { + fprintf(fp, "\t%f %f %f\n", + stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2)); + } + for(i = 0; i < stl->stats.number_of_facets; i++) { + fprintf(fp, "\t3 %d %d %d\n", stl->v_indices[i].vertex[0], + stl->v_indices[i].vertex[1], stl->v_indices[i].vertex[2]); + } + fclose(fp); +} + +void +stl_write_vrml(stl_file *stl, char *file) { + int i; + FILE *fp; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_ascii: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "#VRML V1.0 ascii\n\n"); + fprintf(fp, "Separator {\n"); + fprintf(fp, "\tDEF STLShape ShapeHints {\n"); + fprintf(fp, "\t\tvertexOrdering COUNTERCLOCKWISE\n"); + fprintf(fp, "\t\tfaceType CONVEX\n"); + fprintf(fp, "\t\tshapeType SOLID\n"); + fprintf(fp, "\t\tcreaseAngle 0.0\n"); + fprintf(fp, "\t}\n"); + fprintf(fp, "\tDEF STLModel Separator {\n"); + fprintf(fp, "\t\tDEF STLColor Material {\n"); + fprintf(fp, "\t\t\temissiveColor 0.700000 0.700000 0.000000\n"); + fprintf(fp, "\t\t}\n"); + fprintf(fp, "\t\tDEF STLVertices Coordinate3 {\n"); + fprintf(fp, "\t\t\tpoint [\n"); + + for(i = 0; i < (stl->stats.shared_vertices - 1); i++) { + fprintf(fp, "\t\t\t\t%f %f %f,\n", + stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2)); + } + fprintf(fp, "\t\t\t\t%f %f %f]\n", + stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2)); + fprintf(fp, "\t\t}\n"); + fprintf(fp, "\t\tDEF STLTriangles IndexedFaceSet {\n"); + fprintf(fp, "\t\t\tcoordIndex [\n"); + + for(i = 0; i < (stl->stats.number_of_facets - 1); i++) { + fprintf(fp, "\t\t\t\t%d, %d, %d, -1,\n", stl->v_indices[i].vertex[0], + stl->v_indices[i].vertex[1], stl->v_indices[i].vertex[2]); + } + fprintf(fp, "\t\t\t\t%d, %d, %d, -1]\n", stl->v_indices[i].vertex[0], + stl->v_indices[i].vertex[1], stl->v_indices[i].vertex[2]); + fprintf(fp, "\t\t}\n"); + fprintf(fp, "\t}\n"); + fprintf(fp, "}\n"); + fclose(fp); +} + +void stl_write_obj (stl_file *stl, char *file) { + int i; + FILE* fp; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if (fp == NULL) { + char* error_msg = (char*)malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_ascii: Couldn't open %s for writing", file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + for (i = 0; i < stl->stats.shared_vertices; i++) { + fprintf(fp, "v %f %f %f\n", stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2)); + } + for (i = 0; i < stl->stats.number_of_facets; i++) { + fprintf(fp, "f %d %d %d\n", stl->v_indices[i].vertex[0]+1, stl->v_indices[i].vertex[1]+1, stl->v_indices[i].vertex[2]+1); + } + + fclose(fp); +} diff --git a/src/admesh/stl.h b/src/admesh/stl.h new file mode 100644 index 000000000..096430d15 --- /dev/null +++ b/src/admesh/stl.h @@ -0,0 +1,215 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#ifndef __admesh_stl__ +#define __admesh_stl__ + +#include <stdio.h> +#include <stdint.h> +#include <stddef.h> + +#include <Eigen/Geometry> + +// Size of the binary STL header, free form. +#define LABEL_SIZE 80 +// Binary STL, length of the "number of faces" counter. +#define NUM_FACET_SIZE 4 +// Binary STL, sizeof header + number of faces. +#define HEADER_SIZE 84 +#define STL_MIN_FILE_SIZE 284 +#define ASCII_LINES_PER_FACET 7 + +typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_vertex; +typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_normal; +static_assert(sizeof(stl_vertex) == 12, "size of stl_vertex incorrect"); +static_assert(sizeof(stl_normal) == 12, "size of stl_normal incorrect"); + +typedef struct { + stl_normal normal; + stl_vertex vertex[3]; + char extra[2]; +} stl_facet; +#define SIZEOF_STL_FACET 50 + +static_assert(offsetof(stl_facet, normal) == 0, "stl_facet.normal has correct offset"); +static_assert(offsetof(stl_facet, vertex) == 12, "stl_facet.vertex has correct offset"); +static_assert(offsetof(stl_facet, extra ) == 48, "stl_facet.extra has correct offset"); +static_assert(sizeof(stl_facet) >= SIZEOF_STL_FACET, "size of stl_facet incorrect"); + +typedef enum {binary, ascii, inmemory} stl_type; + +typedef struct { + stl_vertex p1; + stl_vertex p2; + int facet_number; +} stl_edge; + +typedef struct stl_hash_edge { + // Key of a hash edge: sorted vertices of the edge. + unsigned char key[2 * sizeof(stl_vertex)]; + // Compare two keys. + bool operator==(const stl_hash_edge &rhs) { return memcmp(key, rhs.key, sizeof(key)) == 0; } + bool operator!=(const stl_hash_edge &rhs) { return ! (*this == rhs); } + int hash(int M) const { return ((key[0] / 23 + key[1] / 19 + key[2] / 17 + key[3] /13 + key[4] / 11 + key[5] / 7 ) % M); } + // 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; + +typedef struct { + // Index of a neighbor facet. + int neighbor[3]; + // Index of an opposite vertex at the neighbor face. + char which_vertex_not[3]; +} stl_neighbors; + +typedef struct { + int vertex[3]; +} v_indices_struct; + +typedef struct { + char header[81]; + stl_type type; + uint32_t number_of_facets; + stl_vertex max; + stl_vertex min; + stl_vertex size; + float bounding_diameter; + float shortest_edge; + float volume; + unsigned number_of_blocks; + int connected_edges; + int connected_facets_1_edge; + int connected_facets_2_edge; + int connected_facets_3_edge; + int facets_w_1_bad_edge; + int facets_w_2_bad_edge; + int facets_w_3_bad_edge; + int original_num_facets; + int edges_fixed; + int degenerate_facets; + int facets_removed; + int facets_added; + int facets_reversed; + int backwards_edges; + int normals_fixed; + int number_of_parts; + int malloced; + int freed; + int facets_malloced; + int collisions; + int shared_vertices; + int shared_malloced; +} stl_stats; + +typedef struct { + FILE *fp; + stl_facet *facet_start; + stl_edge *edge_start; + stl_hash_edge **heads; + stl_hash_edge *tail; + int M; + stl_neighbors *neighbors_start; + v_indices_struct *v_indices; + stl_vertex *v_shared; + stl_stats stats; + char error; +} stl_file; + + +extern void stl_open(stl_file *stl, const char *file); +extern void stl_close(stl_file *stl); +extern void stl_stats_out(stl_file *stl, FILE *file, char *input_file); +extern void stl_print_edges(stl_file *stl, FILE *file); +extern void stl_print_neighbors(stl_file *stl, char *file); +extern void stl_put_little_int(FILE *fp, int value_in); +extern void stl_put_little_float(FILE *fp, float value_in); +extern void stl_write_ascii(stl_file *stl, const char *file, const char *label); +extern void stl_write_binary(stl_file *stl, const char *file, const char *label); +extern void stl_write_binary_block(stl_file *stl, FILE *fp); +extern void stl_check_facets_exact(stl_file *stl); +extern void stl_check_facets_nearby(stl_file *stl, float tolerance); +extern void stl_remove_unconnected_facets(stl_file *stl); +extern void stl_write_vertex(stl_file *stl, int facet, int vertex); +extern void stl_write_facet(stl_file *stl, char *label, int facet); +extern void stl_write_edge(stl_file *stl, char *label, stl_hash_edge edge); +extern void stl_write_neighbor(stl_file *stl, int facet); +extern void stl_write_quad_object(stl_file *stl, char *file); +extern void stl_verify_neighbors(stl_file *stl); +extern void stl_fill_holes(stl_file *stl); +extern void stl_fix_normal_directions(stl_file *stl); +extern void stl_fix_normal_values(stl_file *stl); +extern void stl_reverse_all_facets(stl_file *stl); +extern void stl_translate(stl_file *stl, float x, float y, float z); +extern void stl_translate_relative(stl_file *stl, float x, float y, float z); +extern void stl_scale_versor(stl_file *stl, const stl_vertex &versor); +inline void stl_scale(stl_file *stl, float factor) { stl_scale_versor(stl, stl_vertex(factor, factor, factor)); } +extern void stl_rotate_x(stl_file *stl, float angle); +extern void stl_rotate_y(stl_file *stl, float angle); +extern void stl_rotate_z(stl_file *stl, float angle); +extern void stl_mirror_xy(stl_file *stl); +extern void stl_mirror_yz(stl_file *stl); +extern void stl_mirror_xz(stl_file *stl); +extern void stl_transform(stl_file *stl, float *trafo3x4); +extern void stl_transform(stl_file *stl, const Eigen::Transform<float, 3, Eigen::Affine, Eigen::DontAlign>& t); +extern void stl_open_merge(stl_file *stl, char *file); +extern void stl_invalidate_shared_vertices(stl_file *stl); +extern void stl_generate_shared_vertices(stl_file *stl); +extern void stl_write_obj(stl_file *stl, char *file); +extern void stl_write_off(stl_file *stl, char *file); +extern void stl_write_dxf(stl_file *stl, char *file, char *label); +extern void stl_write_vrml(stl_file *stl, char *file); +inline void stl_calculate_normal(stl_normal &normal, stl_facet *facet) { + normal = (facet->vertex[1] - facet->vertex[0]).cross(facet->vertex[2] - facet->vertex[0]); +} +inline void stl_normalize_vector(stl_normal &normal) { + double length = normal.cast<double>().norm(); + if (length < 0.000000000001) + normal = stl_normal::Zero(); + else + normal *= (1.0 / length); +} +inline bool stl_vertex_lower(const stl_vertex &a, const stl_vertex &b) { + return (a(0) != b(0)) ? (a(0) < b(0)) : + ((a(1) != b(1)) ? (a(1) < b(1)) : (a(2) < b(2))); +} +extern void stl_calculate_volume(stl_file *stl); + +extern void stl_repair(stl_file *stl, int fixall_flag, int exact_flag, int tolerance_flag, float tolerance, int increment_flag, float increment, int nearby_flag, int iterations, int remove_unconnected_flag, int fill_holes_flag, int normal_directions_flag, int normal_values_flag, int reverse_all_flag, int verbose_flag); + +extern void stl_initialize(stl_file *stl); +extern void stl_count_facets(stl_file *stl, const char *file); +extern void stl_allocate(stl_file *stl); +extern void stl_read(stl_file *stl, int first_facet, bool first); +extern void stl_facet_stats(stl_file *stl, stl_facet facet, bool &first); +extern void stl_reallocate(stl_file *stl); +extern void stl_add_facet(stl_file *stl, stl_facet *new_facet); +extern void stl_get_size(stl_file *stl); + +extern void stl_clear_error(stl_file *stl); +extern int stl_get_error(stl_file *stl); +extern void stl_exit_on_error(stl_file *stl); + +#endif diff --git a/src/admesh/stl_io.cpp b/src/admesh/stl_io.cpp new file mode 100644 index 000000000..81e29d286 --- /dev/null +++ b/src/admesh/stl_io.cpp @@ -0,0 +1,433 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdlib.h> +#include <string.h> +#include "stl.h" + +#include <boost/nowide/cstdio.hpp> +#include <boost/detail/endian.hpp> + +#if !defined(SEEK_SET) +#define SEEK_SET 0 +#define SEEK_CUR 1 +#define SEEK_END 2 +#endif + +void +stl_print_edges(stl_file *stl, FILE *file) { + int i; + int edges_allocated; + + if (stl->error) return; + + edges_allocated = stl->stats.number_of_facets * 3; + for(i = 0; i < edges_allocated; i++) { + fprintf(file, "%d, %f, %f, %f, %f, %f, %f\n", + stl->edge_start[i].facet_number, + stl->edge_start[i].p1(0), stl->edge_start[i].p1(1), + stl->edge_start[i].p1(2), stl->edge_start[i].p2(0), + stl->edge_start[i].p2(1), stl->edge_start[i].p2(2)); + } +} + + +void +stl_stats_out(stl_file *stl, FILE *file, char *input_file) { + if (stl->error) return; + + /* this is here for Slic3r, without our config.h + it won't use this part of the code anyway */ +#ifndef VERSION +#define VERSION "unknown" +#endif + fprintf(file, "\n\ +================= Results produced by ADMesh version " VERSION " ================\n"); + fprintf(file, "\ +Input file : %s\n", input_file); + if(stl->stats.type == binary) { + fprintf(file, "\ +File type : Binary STL file\n"); + } else { + fprintf(file, "\ +File type : ASCII STL file\n"); + } + fprintf(file, "\ +Header : %s\n", stl->stats.header); + fprintf(file, "============== Size ==============\n"); + fprintf(file, "Min X = % f, Max X = % f\n", + stl->stats.min(0), stl->stats.max(0)); + fprintf(file, "Min Y = % f, Max Y = % f\n", + stl->stats.min(1), stl->stats.max(1)); + fprintf(file, "Min Z = % f, Max Z = % f\n", + stl->stats.min(2), stl->stats.max(2)); + + fprintf(file, "\ +========= Facet Status ========== Original ============ Final ====\n"); + fprintf(file, "\ +Number of facets : %5d %5d\n", + stl->stats.original_num_facets, stl->stats.number_of_facets); + fprintf(file, "\ +Facets with 1 disconnected edge : %5d %5d\n", + stl->stats.facets_w_1_bad_edge, stl->stats.connected_facets_2_edge - + stl->stats.connected_facets_3_edge); + fprintf(file, "\ +Facets with 2 disconnected edges : %5d %5d\n", + stl->stats.facets_w_2_bad_edge, stl->stats.connected_facets_1_edge - + stl->stats.connected_facets_2_edge); + fprintf(file, "\ +Facets with 3 disconnected edges : %5d %5d\n", + stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets - + stl->stats.connected_facets_1_edge); + fprintf(file, "\ +Total disconnected facets : %5d %5d\n", + stl->stats.facets_w_1_bad_edge + stl->stats.facets_w_2_bad_edge + + stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets - + stl->stats.connected_facets_3_edge); + + fprintf(file, + "=== Processing Statistics === ===== Other Statistics =====\n"); + fprintf(file, "\ +Number of parts : %5d Volume : % f\n", + stl->stats.number_of_parts, stl->stats.volume); + fprintf(file, "\ +Degenerate facets : %5d\n", stl->stats.degenerate_facets); + fprintf(file, "\ +Edges fixed : %5d\n", stl->stats.edges_fixed); + fprintf(file, "\ +Facets removed : %5d\n", stl->stats.facets_removed); + fprintf(file, "\ +Facets added : %5d\n", stl->stats.facets_added); + fprintf(file, "\ +Facets reversed : %5d\n", stl->stats.facets_reversed); + fprintf(file, "\ +Backwards edges : %5d\n", stl->stats.backwards_edges); + fprintf(file, "\ +Normals fixed : %5d\n", stl->stats.normals_fixed); +} + +void +stl_write_ascii(stl_file *stl, const char *file, const char *label) { + int i; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + FILE *fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_ascii: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "solid %s\n", label); + + for(i = 0; i < stl->stats.number_of_facets; i++) { + fprintf(fp, " facet normal % .8E % .8E % .8E\n", + stl->facet_start[i].normal(0), stl->facet_start[i].normal(1), + stl->facet_start[i].normal(2)); + fprintf(fp, " outer loop\n"); + fprintf(fp, " vertex % .8E % .8E % .8E\n", + stl->facet_start[i].vertex[0](0), stl->facet_start[i].vertex[0](1), + stl->facet_start[i].vertex[0](2)); + fprintf(fp, " vertex % .8E % .8E % .8E\n", + stl->facet_start[i].vertex[1](0), stl->facet_start[i].vertex[1](1), + stl->facet_start[i].vertex[1](2)); + fprintf(fp, " vertex % .8E % .8E % .8E\n", + stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1), + stl->facet_start[i].vertex[2](2)); + fprintf(fp, " endloop\n"); + fprintf(fp, " endfacet\n"); + } + + fprintf(fp, "endsolid %s\n", label); + + fclose(fp); +} + +void +stl_print_neighbors(stl_file *stl, char *file) { + int i; + FILE *fp; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_print_neighbors: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + for(i = 0; i < stl->stats.number_of_facets; i++) { + fprintf(fp, "%d, %d,%d, %d,%d, %d,%d\n", + i, + stl->neighbors_start[i].neighbor[0], + (int)stl->neighbors_start[i].which_vertex_not[0], + stl->neighbors_start[i].neighbor[1], + (int)stl->neighbors_start[i].which_vertex_not[1], + stl->neighbors_start[i].neighbor[2], + (int)stl->neighbors_start[i].which_vertex_not[2]); + } + fclose(fp); +} + +#ifndef BOOST_LITTLE_ENDIAN +// Swap a buffer of 32bit data from little endian to big endian and vice versa. +void stl_internal_reverse_quads(char *buf, size_t cnt) +{ + for (size_t i = 0; i < cnt; i += 4) { + std::swap(buf[i], buf[i+3]); + std::swap(buf[i+1], buf[i+2]); + } +} +#endif + +void +stl_write_binary(stl_file *stl, const char *file, const char *label) { + FILE *fp; + int i; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "wb"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_binary: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "%s", label); + for(i = strlen(label); i < LABEL_SIZE; i++) putc(0, fp); + + fseek(fp, LABEL_SIZE, SEEK_SET); +#ifdef BOOST_LITTLE_ENDIAN + fwrite(&stl->stats.number_of_facets, 4, 1, fp); + for (i = 0; i < stl->stats.number_of_facets; ++ i) + fwrite(stl->facet_start + i, SIZEOF_STL_FACET, 1, fp); +#else /* BOOST_LITTLE_ENDIAN */ + char buffer[50]; + // Convert the number of facets to little endian. + memcpy(buffer, &stl->stats.number_of_facets, 4); + stl_internal_reverse_quads(buffer, 4); + fwrite(buffer, 4, 1, fp); + for (i = 0; i < stl->stats.number_of_facets; ++ i) { + memcpy(buffer, stl->facet_start + i, 50); + // Convert to little endian. + stl_internal_reverse_quads(buffer, 48); + fwrite(buffer, SIZEOF_STL_FACET, 1, fp); + } +#endif /* BOOST_LITTLE_ENDIAN */ + fclose(fp); +} + +void +stl_write_vertex(stl_file *stl, int facet, int vertex) { + if (stl->error) return; + printf(" vertex %d/%d % .8E % .8E % .8E\n", vertex, facet, + stl->facet_start[facet].vertex[vertex](0), + stl->facet_start[facet].vertex[vertex](1), + stl->facet_start[facet].vertex[vertex](2)); +} + +void +stl_write_facet(stl_file *stl, char *label, int facet) { + if (stl->error) return; + printf("facet (%d)/ %s\n", facet, label); + stl_write_vertex(stl, facet, 0); + stl_write_vertex(stl, facet, 1); + stl_write_vertex(stl, facet, 2); +} + +void +stl_write_edge(stl_file *stl, char *label, stl_hash_edge edge) { + if (stl->error) return; + printf("edge (%d)/(%d) %s\n", edge.facet_number, edge.which_edge, label); + if(edge.which_edge < 3) { + stl_write_vertex(stl, edge.facet_number, edge.which_edge % 3); + stl_write_vertex(stl, edge.facet_number, (edge.which_edge + 1) % 3); + } else { + stl_write_vertex(stl, edge.facet_number, (edge.which_edge + 1) % 3); + stl_write_vertex(stl, edge.facet_number, edge.which_edge % 3); + } +} + +void +stl_write_neighbor(stl_file *stl, int facet) { + if (stl->error) return; + printf("Neighbors %d: %d, %d, %d ; %d, %d, %d\n", facet, + stl->neighbors_start[facet].neighbor[0], + stl->neighbors_start[facet].neighbor[1], + stl->neighbors_start[facet].neighbor[2], + stl->neighbors_start[facet].which_vertex_not[0], + stl->neighbors_start[facet].which_vertex_not[1], + stl->neighbors_start[facet].which_vertex_not[2]); +} + +void +stl_write_quad_object(stl_file *stl, char *file) { + FILE *fp; + int i; + int j; + char *error_msg; + stl_vertex connect_color = stl_vertex::Zero(); + stl_vertex uncon_1_color = stl_vertex::Zero(); + stl_vertex uncon_2_color = stl_vertex::Zero(); + stl_vertex uncon_3_color = stl_vertex::Zero(); + stl_vertex color; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_quad_object: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "CQUAD\n"); + for(i = 0; i < stl->stats.number_of_facets; i++) { + j = ((stl->neighbors_start[i].neighbor[0] == -1) + + (stl->neighbors_start[i].neighbor[1] == -1) + + (stl->neighbors_start[i].neighbor[2] == -1)); + if(j == 0) { + color = connect_color; + } else if(j == 1) { + color = uncon_1_color; + } else if(j == 2) { + color = uncon_2_color; + } else { + color = uncon_3_color; + } + fprintf(fp, "%f %f %f %1.1f %1.1f %1.1f 1\n", + stl->facet_start[i].vertex[0](0), + stl->facet_start[i].vertex[0](1), + stl->facet_start[i].vertex[0](2), color(0), color(1), color(2)); + fprintf(fp, "%f %f %f %1.1f %1.1f %1.1f 1\n", + stl->facet_start[i].vertex[1](0), + stl->facet_start[i].vertex[1](1), + stl->facet_start[i].vertex[1](2), color(0), color(1), color(2)); + fprintf(fp, "%f %f %f %1.1f %1.1f %1.1f 1\n", + stl->facet_start[i].vertex[2](0), + stl->facet_start[i].vertex[2](1), + stl->facet_start[i].vertex[2](2), color(0), color(1), color(2)); + fprintf(fp, "%f %f %f %1.1f %1.1f %1.1f 1\n", + stl->facet_start[i].vertex[2](0), + stl->facet_start[i].vertex[2](1), + stl->facet_start[i].vertex[2](2), color(0), color(1), color(2)); + } + fclose(fp); +} + +void +stl_write_dxf(stl_file *stl, char *file, char *label) { + int i; + FILE *fp; + char *error_msg; + + if (stl->error) return; + + /* Open the file */ + fp = boost::nowide::fopen(file, "w"); + if(fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_write_ascii: Couldn't open %s for writing", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + fprintf(fp, "999\n%s\n", label); + fprintf(fp, "0\nSECTION\n2\nHEADER\n0\nENDSEC\n"); + fprintf(fp, "0\nSECTION\n2\nTABLES\n0\nTABLE\n2\nLAYER\n70\n1\n\ +0\nLAYER\n2\n0\n70\n0\n62\n7\n6\nCONTINUOUS\n0\nENDTAB\n0\nENDSEC\n"); + fprintf(fp, "0\nSECTION\n2\nBLOCKS\n0\nENDSEC\n"); + + fprintf(fp, "0\nSECTION\n2\nENTITIES\n"); + + for(i = 0; i < stl->stats.number_of_facets; i++) { + fprintf(fp, "0\n3DFACE\n8\n0\n"); + fprintf(fp, "10\n%f\n20\n%f\n30\n%f\n", + stl->facet_start[i].vertex[0](0), stl->facet_start[i].vertex[0](1), + stl->facet_start[i].vertex[0](2)); + fprintf(fp, "11\n%f\n21\n%f\n31\n%f\n", + stl->facet_start[i].vertex[1](0), stl->facet_start[i].vertex[1](1), + stl->facet_start[i].vertex[1](2)); + fprintf(fp, "12\n%f\n22\n%f\n32\n%f\n", + stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1), + stl->facet_start[i].vertex[2](2)); + fprintf(fp, "13\n%f\n23\n%f\n33\n%f\n", + stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1), + stl->facet_start[i].vertex[2](2)); + } + + fprintf(fp, "0\nENDSEC\n0\nEOF\n"); + + fclose(fp); +} + +void +stl_clear_error(stl_file *stl) { + stl->error = 0; +} + +void +stl_exit_on_error(stl_file *stl) { + if (!stl->error) return; + stl->error = 0; + stl_close(stl); + exit(1); +} + +int +stl_get_error(stl_file *stl) { + return stl->error; +} diff --git a/src/admesh/stlinit.cpp b/src/admesh/stlinit.cpp new file mode 100644 index 000000000..e2939b8af --- /dev/null +++ b/src/admesh/stlinit.cpp @@ -0,0 +1,382 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> +#include <assert.h> + +#include <boost/nowide/cstdio.hpp> +#include <boost/detail/endian.hpp> + +#include "stl.h" + +#ifndef SEEK_SET +#error "SEEK_SET not defined" +#endif + +void +stl_open(stl_file *stl, const char *file) { + stl_initialize(stl); + stl_count_facets(stl, file); + stl_allocate(stl); + stl_read(stl, 0, true); + if (!stl->error) fclose(stl->fp); +} + + +void +stl_initialize(stl_file *stl) { + memset(stl, 0, sizeof(stl_file)); + stl->stats.volume = -1.0; +} + +#ifndef BOOST_LITTLE_ENDIAN +extern void stl_internal_reverse_quads(char *buf, size_t cnt); +#endif /* BOOST_LITTLE_ENDIAN */ + +void +stl_count_facets(stl_file *stl, const char *file) { + long file_size; + uint32_t header_num_facets; + uint32_t num_facets; + int i; + size_t s; + unsigned char chtest[128]; + int num_lines = 1; + char *error_msg; + + if (stl->error) return; + + /* Open the file in binary mode first */ + stl->fp = boost::nowide::fopen(file, "rb"); + if(stl->fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_initialize: Couldn't open %s for reading", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + /* Find size of file */ + fseek(stl->fp, 0, SEEK_END); + file_size = ftell(stl->fp); + + /* Check for binary or ASCII file */ + fseek(stl->fp, HEADER_SIZE, SEEK_SET); + if (!fread(chtest, sizeof(chtest), 1, stl->fp)) { + perror("The input is an empty file"); + stl->error = 1; + return; + } + stl->stats.type = ascii; + for(s = 0; s < sizeof(chtest); s++) { + if(chtest[s] > 127) { + stl->stats.type = binary; + break; + } + } + rewind(stl->fp); + + /* Get the header and the number of facets in the .STL file */ + /* If the .STL file is binary, then do the following */ + if(stl->stats.type == binary) { + /* Test if the STL file has the right size */ + if(((file_size - HEADER_SIZE) % SIZEOF_STL_FACET != 0) + || (file_size < STL_MIN_FILE_SIZE)) { + fprintf(stderr, "The file %s has the wrong size.\n", file); + stl->error = 1; + return; + } + num_facets = (file_size - HEADER_SIZE) / SIZEOF_STL_FACET; + + /* Read the header */ + if (fread(stl->stats.header, LABEL_SIZE, 1, stl->fp) > 79) { + stl->stats.header[80] = '\0'; + } + + /* Read the int following the header. This should contain # of facets */ + bool header_num_faces_read = fread(&header_num_facets, sizeof(uint32_t), 1, stl->fp); +#ifndef BOOST_LITTLE_ENDIAN + // Convert from little endian to big endian. + stl_internal_reverse_quads((char*)&header_num_facets, 4); +#endif /* BOOST_LITTLE_ENDIAN */ + if (! header_num_faces_read || num_facets != header_num_facets) { + fprintf(stderr, + "Warning: File size doesn't match number of facets in the header\n"); + } + } + /* Otherwise, if the .STL file is ASCII, then do the following */ + else { + /* Reopen the file in text mode (for getting correct newlines on Windows) */ + // fix to silence a warning about unused return value. + // obviously if it fails we have problems.... + stl->fp = boost::nowide::freopen(file, "r", stl->fp); + + // do another null check to be safe + if(stl->fp == NULL) { + error_msg = (char*) + malloc(81 + strlen(file)); /* Allow 80 chars+file size for message */ + sprintf(error_msg, "stl_initialize: Couldn't open %s for reading", + file); + perror(error_msg); + free(error_msg); + stl->error = 1; + return; + } + + /* Find the number of facets */ + char linebuf[100]; + while (fgets(linebuf, 100, stl->fp) != NULL) { + /* don't count short lines */ + if (strlen(linebuf) <= 4) continue; + + /* skip solid/endsolid lines as broken STL file generators may put several of them */ + if (strncmp(linebuf, "solid", 5) == 0 || strncmp(linebuf, "endsolid", 8) == 0) continue; + + ++num_lines; + } + + rewind(stl->fp); + + /* Get the header */ + for(i = 0; + (i < 80) && (stl->stats.header[i] = getc(stl->fp)) != '\n'; i++); + stl->stats.header[i] = '\0'; /* Lose the '\n' */ + stl->stats.header[80] = '\0'; + + num_facets = num_lines / ASCII_LINES_PER_FACET; + } + stl->stats.number_of_facets += num_facets; + stl->stats.original_num_facets = stl->stats.number_of_facets; +} + +void +stl_allocate(stl_file *stl) { + if (stl->error) return; + + /* Allocate memory for the entire .STL file */ + stl->facet_start = (stl_facet*)calloc(stl->stats.number_of_facets, + sizeof(stl_facet)); + if(stl->facet_start == NULL) perror("stl_initialize"); + stl->stats.facets_malloced = stl->stats.number_of_facets; + + /* Allocate memory for the neighbors list */ + stl->neighbors_start = (stl_neighbors*) + calloc(stl->stats.number_of_facets, sizeof(stl_neighbors)); + if(stl->facet_start == NULL) perror("stl_initialize"); +} + +void +stl_open_merge(stl_file *stl, char *file_to_merge) { + int num_facets_so_far; + stl_type origStlType; + FILE *origFp; + stl_file stl_to_merge; + + if (stl->error) return; + + /* Record how many facets we have so far from the first file. We will start putting + facets in the next position. Since we're 0-indexed, it'l be the same position. */ + num_facets_so_far = stl->stats.number_of_facets; + + /* Record the file type we started with: */ + origStlType=stl->stats.type; + /* Record the file pointer too: */ + origFp=stl->fp; + + /* Initialize the sturucture with zero stats, header info and sizes: */ + stl_initialize(&stl_to_merge); + stl_count_facets(&stl_to_merge, file_to_merge); + + /* Copy what we need to into stl so that we can read the file_to_merge directly into it + using stl_read: Save the rest of the valuable info: */ + stl->stats.type=stl_to_merge.stats.type; + stl->fp=stl_to_merge.fp; + + /* Add the number of facets we already have in stl with what we we found in stl_to_merge but + haven't read yet. */ + stl->stats.number_of_facets=num_facets_so_far+stl_to_merge.stats.number_of_facets; + + /* Allocate enough room for stl->stats.number_of_facets facets and neighbors: */ + stl_reallocate(stl); + + /* Read the file to merge directly into stl, adding it to what we have already. + Start at num_facets_so_far, the index to the first unused facet. Also say + that this isn't our first time so we should augment stats like min and max + instead of erasing them. */ + stl_read(stl, num_facets_so_far, false); + + /* Restore the stl information we overwrote (for stl_read) so that it still accurately + reflects the subject part: */ + stl->stats.type=origStlType; + stl->fp=origFp; +} + +extern void +stl_reallocate(stl_file *stl) { + if (stl->error) return; + /* Reallocate more memory for the .STL file(s) */ + stl->facet_start = (stl_facet*)realloc(stl->facet_start, stl->stats.number_of_facets * + sizeof(stl_facet)); + if(stl->facet_start == NULL) perror("stl_initialize"); + stl->stats.facets_malloced = stl->stats.number_of_facets; + + /* Reallocate more memory for the neighbors list */ + stl->neighbors_start = (stl_neighbors*) + realloc(stl->neighbors_start, stl->stats.number_of_facets * + sizeof(stl_neighbors)); + if(stl->facet_start == NULL) perror("stl_initialize"); +} + + +/* Reads the contents of the file pointed to by stl->fp into the stl structure, + starting at facet first_facet. The second argument says if it's our first + time running this for the stl and therefore we should reset our max and min stats. */ +void stl_read(stl_file *stl, int first_facet, bool first) { + stl_facet facet; + int i; + + if (stl->error) return; + + if(stl->stats.type == binary) { + fseek(stl->fp, HEADER_SIZE, SEEK_SET); + } else { + rewind(stl->fp); + } + + char normal_buf[3][32]; + for(i = first_facet; i < stl->stats.number_of_facets; i++) { + if(stl->stats.type == binary) + /* Read a single facet from a binary .STL file */ + { + /* we assume little-endian architecture! */ + if (fread(&facet, 1, SIZEOF_STL_FACET, stl->fp) != SIZEOF_STL_FACET) { + stl->error = 1; + return; + } +#ifndef BOOST_LITTLE_ENDIAN + // Convert the loaded little endian data to big endian. + stl_internal_reverse_quads((char*)&facet, 48); +#endif /* BOOST_LITTLE_ENDIAN */ + } else + /* Read a single facet from an ASCII .STL file */ + { + // skip solid/endsolid + // (in this order, otherwise it won't work when they are paired in the middle of a file) + fscanf(stl->fp, "endsolid%*[^\n]\n"); + fscanf(stl->fp, "solid%*[^\n]\n"); // name might contain spaces so %*s doesn't work and it also can be empty (just "solid") + // Leading space in the fscanf format skips all leading white spaces including numerous new lines and tabs. + int res_normal = fscanf(stl->fp, " facet normal %31s %31s %31s", normal_buf[0], normal_buf[1], normal_buf[2]); + assert(res_normal == 3); + int res_outer_loop = fscanf(stl->fp, " outer loop"); + assert(res_outer_loop == 0); + int res_vertex1 = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[0](0), &facet.vertex[0](1), &facet.vertex[0](2)); + assert(res_vertex1 == 3); + int res_vertex2 = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[1](0), &facet.vertex[1](1), &facet.vertex[1](2)); + assert(res_vertex2 == 3); + int res_vertex3 = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[2](0), &facet.vertex[2](1), &facet.vertex[2](2)); + assert(res_vertex3 == 3); + int res_endloop = fscanf(stl->fp, " endloop"); + assert(res_endloop == 0); + // There is a leading and trailing white space around endfacet to eat up all leading and trailing white spaces including numerous tabs and new lines. + int res_endfacet = fscanf(stl->fp, " endfacet "); + if (res_normal != 3 || res_outer_loop != 0 || res_vertex1 != 3 || res_vertex2 != 3 || res_vertex3 != 3 || res_endloop != 0 || res_endfacet != 0) { + perror("Something is syntactically very wrong with this ASCII STL!"); + stl->error = 1; + return; + } + + // The facet normal has been parsed as a single string as to workaround for not a numbers in the normal definition. + if (sscanf(normal_buf[0], "%f", &facet.normal(0)) != 1 || + sscanf(normal_buf[1], "%f", &facet.normal(1)) != 1 || + sscanf(normal_buf[2], "%f", &facet.normal(2)) != 1) { + // Normal was mangled. Maybe denormals or "not a number" were stored? + // Just reset the normal and silently ignore it. + memset(&facet.normal, 0, sizeof(facet.normal)); + } + } + +#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](0) > -1e-12f && facet.vertex[j](0) < 1e-12f) + printf("stl_read: facet %d(0) = %e\r\n", j, facet.vertex[j](0)); + if (facet.vertex[j](1) > -1e-12f && facet.vertex[j](1) < 1e-12f) + printf("stl_read: facet %d(1) = %e\r\n", j, facet.vertex[j](1)); + if (facet.vertex[j](2) > -1e-12f && facet.vertex[j](2) < 1e-12f) + printf("stl_read: facet %d(2) = %e\r\n", j, facet.vertex[j](2)); + } +#endif + + /* Write the facet into memory. */ + stl->facet_start[i] = facet; + stl_facet_stats(stl, facet, first); + } + stl->stats.size = stl->stats.max - stl->stats.min; + stl->stats.bounding_diameter = stl->stats.size.norm(); +} + +void stl_facet_stats(stl_file *stl, stl_facet facet, bool &first) +{ + if (stl->error) + return; + + // While we are going through all of the facets, let's find the + // maximum and minimum values for x, y, and z + + if (first) { + // Initialize the max and min values the first time through + stl->stats.min = facet.vertex[0]; + stl->stats.max = facet.vertex[0]; + stl_vertex diff = (facet.vertex[1] - facet.vertex[0]).cwiseAbs(); + stl->stats.shortest_edge = std::max(diff(0), std::max(diff(1), diff(2))); + first = false; + } + + // Now find the max and min values. + for (size_t i = 0; i < 3; ++ i) { + stl->stats.min = stl->stats.min.cwiseMin(facet.vertex[i]); + stl->stats.max = stl->stats.max.cwiseMax(facet.vertex[i]); + } +} + +void +stl_close(stl_file *stl) { + if (stl->error) return; + + if(stl->neighbors_start != NULL) + free(stl->neighbors_start); + if(stl->facet_start != NULL) + free(stl->facet_start); + if(stl->v_indices != NULL) + free(stl->v_indices); + if(stl->v_shared != NULL) + free(stl->v_shared); +} + diff --git a/src/admesh/util.cpp b/src/admesh/util.cpp new file mode 100644 index 000000000..cc104fdd1 --- /dev/null +++ b/src/admesh/util.cpp @@ -0,0 +1,518 @@ +/* ADMesh -- process triangulated solid meshes + * Copyright (C) 1995, 1996 Anthony D. Martin <amartin@engr.csulb.edu> + * Copyright (C) 2013, 2014 several contributors, see AUTHORS + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Questions, comments, suggestions, etc to + * https://github.com/admesh/admesh/issues + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "stl.h" + +static void stl_rotate(float *x, float *y, const double c, const double s); +static float get_area(stl_facet *facet); +static float get_volume(stl_file *stl); + + +void +stl_verify_neighbors(stl_file *stl) { + int i; + int j; + stl_edge edge_a; + stl_edge edge_b; + int neighbor; + int vnot; + + if (stl->error) return; + + stl->stats.backwards_edges = 0; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + for(j = 0; j < 3; j++) { + edge_a.p1 = stl->facet_start[i].vertex[j]; + edge_a.p2 = stl->facet_start[i].vertex[(j + 1) % 3]; + neighbor = stl->neighbors_start[i].neighbor[j]; + vnot = stl->neighbors_start[i].which_vertex_not[j]; + + if(neighbor == -1) + continue; /* this edge has no neighbor... Continue. */ + if(vnot < 3) { + edge_b.p1 = stl->facet_start[neighbor].vertex[(vnot + 2) % 3]; + edge_b.p2 = stl->facet_start[neighbor].vertex[(vnot + 1) % 3]; + } else { + stl->stats.backwards_edges += 1; + edge_b.p1 = stl->facet_start[neighbor].vertex[(vnot + 1) % 3]; + edge_b.p2 = stl->facet_start[neighbor].vertex[(vnot + 2) % 3]; + } + if (edge_a.p1 != edge_b.p1 || edge_a.p2 != edge_b.p2) { + /* These edges should match but they don't. Print results. */ + printf("edge %d of facet %d doesn't match edge %d of facet %d\n", + j, i, vnot + 1, neighbor); + stl_write_facet(stl, (char*)"first facet", i); + stl_write_facet(stl, (char*)"second facet", neighbor); + } + } + } +} + +void stl_translate(stl_file *stl, float x, float y, float z) +{ + if (stl->error) + return; + + stl_vertex new_min(x, y, z); + stl_vertex shift = new_min - stl->stats.min; + for (int i = 0; i < stl->stats.number_of_facets; ++ i) + for (int j = 0; j < 3; ++ j) + stl->facet_start[i].vertex[j] += shift; + stl->stats.min = new_min; + stl->stats.max += shift; + stl_invalidate_shared_vertices(stl); +} + +/* Translates the stl by x,y,z, relatively from wherever it is currently */ +void stl_translate_relative(stl_file *stl, float x, float y, float z) +{ + if (stl->error) + return; + + stl_vertex shift(x, y, z); + for (int i = 0; i < stl->stats.number_of_facets; ++ i) + for (int j = 0; j < 3; ++ j) + stl->facet_start[i].vertex[j] += shift; + stl->stats.min += shift; + stl->stats.max += shift; + stl_invalidate_shared_vertices(stl); +} + +void stl_scale_versor(stl_file *stl, const stl_vertex &versor) +{ + if (stl->error) + return; + + // Scale extents. + auto s = versor.array(); + stl->stats.min.array() *= s; + stl->stats.max.array() *= s; + // Scale size. + stl->stats.size.array() *= s; + // Scale volume. + if (stl->stats.volume > 0.0) + stl->stats.volume *= versor(0) * versor(1) * versor(2); + // Scale the mesh. + for (int i = 0; i < stl->stats.number_of_facets; ++ i) + for (int j = 0; j < 3; ++ j) + stl->facet_start[i].vertex[j].array() *= s; + stl_invalidate_shared_vertices(stl); +} + +static void calculate_normals(stl_file *stl) +{ + if (stl->error) + return; + + stl_normal normal; + for(uint32_t i = 0; i < stl->stats.number_of_facets; i++) { + stl_calculate_normal(normal, &stl->facet_start[i]); + stl_normalize_vector(normal); + stl->facet_start[i].normal = normal; + } +} + +void stl_transform(stl_file *stl, float *trafo3x4) { + int i_face, i_vertex; + if (stl->error) + return; + for (i_face = 0; i_face < stl->stats.number_of_facets; ++ i_face) { + stl_vertex *vertices = stl->facet_start[i_face].vertex; + for (i_vertex = 0; i_vertex < 3; ++ i_vertex) { + stl_vertex &v_dst = vertices[i_vertex]; + stl_vertex v_src = v_dst; + v_dst(0) = trafo3x4[0] * v_src(0) + trafo3x4[1] * v_src(1) + trafo3x4[2] * v_src(2) + trafo3x4[3]; + v_dst(1) = trafo3x4[4] * v_src(0) + trafo3x4[5] * v_src(1) + trafo3x4[6] * v_src(2) + trafo3x4[7]; + v_dst(2) = trafo3x4[8] * v_src(0) + trafo3x4[9] * v_src(1) + trafo3x4[10] * v_src(2) + trafo3x4[11]; + } + } + stl_get_size(stl); + calculate_normals(stl); +} + +void stl_transform(stl_file *stl, const Eigen::Transform<float, 3, Eigen::Affine, Eigen::DontAlign>& t) +{ + if (stl->error) + return; + + unsigned int vertices_count = 3 * (unsigned int)stl->stats.number_of_facets; + if (vertices_count == 0) + return; + + Eigen::MatrixXf src_vertices(3, vertices_count); + stl_facet* facet_ptr = stl->facet_start; + unsigned int v_id = 0; + while (facet_ptr < stl->facet_start + stl->stats.number_of_facets) + { + for (int i = 0; i < 3; ++i) + { + ::memcpy((void*)src_vertices.col(v_id).data(), (const void*)&facet_ptr->vertex[i], 3 * sizeof(float)); + ++v_id; + } + facet_ptr += 1; + } + + Eigen::MatrixXf dst_vertices(3, vertices_count); + dst_vertices = t * src_vertices.colwise().homogeneous(); + + facet_ptr = stl->facet_start; + v_id = 0; + while (facet_ptr < stl->facet_start + stl->stats.number_of_facets) + { + for (int i = 0; i < 3; ++i) + { + ::memcpy((void*)&facet_ptr->vertex[i], (const void*)dst_vertices.col(v_id).data(), 3 * sizeof(float)); + ++v_id; + } + facet_ptr += 1; + } + + stl_get_size(stl); + calculate_normals(stl); +} + +void +stl_rotate_x(stl_file *stl, float angle) { + int i; + int j; + double radian_angle = (angle / 180.0) * M_PI; + double c = cos(radian_angle); + double s = sin(radian_angle); + + if (stl->error) return; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + for(j = 0; j < 3; j++) { + stl_rotate(&stl->facet_start[i].vertex[j](1), + &stl->facet_start[i].vertex[j](2), c, s); + } + } + stl_get_size(stl); + calculate_normals(stl); +} + +void +stl_rotate_y(stl_file *stl, float angle) { + int i; + int j; + double radian_angle = (angle / 180.0) * M_PI; + double c = cos(radian_angle); + double s = sin(radian_angle); + + if (stl->error) return; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + for(j = 0; j < 3; j++) { + stl_rotate(&stl->facet_start[i].vertex[j](2), + &stl->facet_start[i].vertex[j](0), c, s); + } + } + stl_get_size(stl); + calculate_normals(stl); +} + +void +stl_rotate_z(stl_file *stl, float angle) { + int i; + int j; + double radian_angle = (angle / 180.0) * M_PI; + double c = cos(radian_angle); + double s = sin(radian_angle); + + if (stl->error) return; + + for(i = 0; i < stl->stats.number_of_facets; i++) { + for(j = 0; j < 3; j++) { + stl_rotate(&stl->facet_start[i].vertex[j](0), + &stl->facet_start[i].vertex[j](1), c, s); + } + } + stl_get_size(stl); + calculate_normals(stl); +} + + + +static void +stl_rotate(float *x, float *y, const double c, const double s) { + double xold = *x; + double yold = *y; + *x = float(c * xold - s * yold); + *y = float(s * xold + c * yold); +} + +void stl_get_size(stl_file *stl) +{ + if (stl->error || stl->stats.number_of_facets == 0) + return; + stl->stats.min = stl->facet_start[0].vertex[0]; + stl->stats.max = stl->stats.min; + for (int i = 0; i < stl->stats.number_of_facets; ++ i) { + const stl_facet &face = stl->facet_start[i]; + for (int j = 0; j < 3; ++ j) { + stl->stats.min = stl->stats.min.cwiseMin(face.vertex[j]); + stl->stats.max = stl->stats.max.cwiseMax(face.vertex[j]); + } + } + stl->stats.size = stl->stats.max - stl->stats.min; + stl->stats.bounding_diameter = stl->stats.size.norm(); +} + +void stl_mirror_xy(stl_file *stl) +{ + if (stl->error) + return; + + for(int i = 0; i < stl->stats.number_of_facets; i++) { + for(int j = 0; j < 3; j++) { + stl->facet_start[i].vertex[j](2) *= -1.0; + } + } + float temp_size = stl->stats.min(2); + stl->stats.min(2) = stl->stats.max(2); + stl->stats.max(2) = temp_size; + stl->stats.min(2) *= -1.0; + stl->stats.max(2) *= -1.0; + stl_reverse_all_facets(stl); + stl->stats.facets_reversed -= stl->stats.number_of_facets; /* for not altering stats */ +} + +void stl_mirror_yz(stl_file *stl) +{ + if (stl->error) return; + + for (int i = 0; i < stl->stats.number_of_facets; i++) { + for (int j = 0; j < 3; j++) { + stl->facet_start[i].vertex[j](0) *= -1.0; + } + } + float temp_size = stl->stats.min(0); + stl->stats.min(0) = stl->stats.max(0); + stl->stats.max(0) = temp_size; + stl->stats.min(0) *= -1.0; + stl->stats.max(0) *= -1.0; + stl_reverse_all_facets(stl); + stl->stats.facets_reversed -= stl->stats.number_of_facets; /* for not altering stats */ +} + +void stl_mirror_xz(stl_file *stl) +{ + if (stl->error) + return; + + for (int i = 0; i < stl->stats.number_of_facets; i++) { + for (int j = 0; j < 3; j++) { + stl->facet_start[i].vertex[j](1) *= -1.0; + } + } + float temp_size = stl->stats.min(1); + stl->stats.min(1) = stl->stats.max(1); + stl->stats.max(1) = temp_size; + stl->stats.min(1) *= -1.0; + stl->stats.max(1) *= -1.0; + stl_reverse_all_facets(stl); + stl->stats.facets_reversed -= stl->stats.number_of_facets; /* for not altering stats */ +} + +static float get_volume(stl_file *stl) +{ + if (stl->error) + return 0; + + // Choose a point, any point as the reference. + stl_vertex p0 = stl->facet_start[0].vertex[0]; + float volume = 0.f; + for(uint32_t i = 0; i < stl->stats.number_of_facets; ++ i) { + // Do dot product to get distance from point to plane. + float height = stl->facet_start[i].normal.dot(stl->facet_start[i].vertex[0] - p0); + float area = get_area(&stl->facet_start[i]); + volume += (area * height) / 3.0f; + } + return volume; +} + +void stl_calculate_volume(stl_file *stl) +{ + if (stl->error) return; + stl->stats.volume = get_volume(stl); + if(stl->stats.volume < 0.0) { + stl_reverse_all_facets(stl); + stl->stats.volume = -stl->stats.volume; + } +} + +static float get_area(stl_facet *facet) +{ + /* cast to double before calculating cross product because large coordinates + can result in overflowing product + (bad area is responsible for bad volume and bad facets reversal) */ + double cross[3][3]; + for (int i = 0; i < 3; i++) { + cross[i][0]=(((double)facet->vertex[i](1) * (double)facet->vertex[(i + 1) % 3](2)) - + ((double)facet->vertex[i](2) * (double)facet->vertex[(i + 1) % 3](1))); + cross[i][1]=(((double)facet->vertex[i](2) * (double)facet->vertex[(i + 1) % 3](0)) - + ((double)facet->vertex[i](0) * (double)facet->vertex[(i + 1) % 3](2))); + cross[i][2]=(((double)facet->vertex[i](0) * (double)facet->vertex[(i + 1) % 3](1)) - + ((double)facet->vertex[i](1) * (double)facet->vertex[(i + 1) % 3](0))); + } + + stl_normal sum; + sum(0) = cross[0][0] + cross[1][0] + cross[2][0]; + sum(1) = cross[0][1] + cross[1][1] + cross[2][1]; + sum(2) = cross[0][2] + cross[1][2] + cross[2][2]; + + // This should already be done. But just in case, let's do it again. + //FIXME this is questionable. the "sum" normal should be accurate, while the normal "n" may be calculated with a low accuracy. + stl_normal n; + stl_calculate_normal(n, facet); + stl_normalize_vector(n); + return 0.5f * n.dot(sum); +} + +void stl_repair(stl_file *stl, + int fixall_flag, + int exact_flag, + int tolerance_flag, + float tolerance, + int increment_flag, + float increment, + int nearby_flag, + int iterations, + int remove_unconnected_flag, + int fill_holes_flag, + int normal_directions_flag, + int normal_values_flag, + int reverse_all_flag, + int verbose_flag) { + + int i; + int last_edges_fixed = 0; + + if (stl->error) return; + + if(exact_flag || fixall_flag || nearby_flag || remove_unconnected_flag + || fill_holes_flag || normal_directions_flag) { + if (verbose_flag) + printf("Checking exact...\n"); + exact_flag = 1; + stl_check_facets_exact(stl); + stl->stats.facets_w_1_bad_edge = + (stl->stats.connected_facets_2_edge - + stl->stats.connected_facets_3_edge); + stl->stats.facets_w_2_bad_edge = + (stl->stats.connected_facets_1_edge - + stl->stats.connected_facets_2_edge); + stl->stats.facets_w_3_bad_edge = + (stl->stats.number_of_facets - + stl->stats.connected_facets_1_edge); + } + + if(nearby_flag || fixall_flag) { + if(!tolerance_flag) { + tolerance = stl->stats.shortest_edge; + } + if(!increment_flag) { + increment = stl->stats.bounding_diameter / 10000.0; + } + + if(stl->stats.connected_facets_3_edge < stl->stats.number_of_facets) { + for(i = 0; i < iterations; i++) { + if(stl->stats.connected_facets_3_edge < + stl->stats.number_of_facets) { + if (verbose_flag) + printf("\ +Checking nearby. Tolerance= %f Iteration=%d of %d...", + tolerance, i + 1, iterations); + stl_check_facets_nearby(stl, tolerance); + if (verbose_flag) + printf(" Fixed %d edges.\n", + stl->stats.edges_fixed - last_edges_fixed); + last_edges_fixed = stl->stats.edges_fixed; + tolerance += increment; + } else { + if (verbose_flag) + printf("\ +All facets connected. No further nearby check necessary.\n"); + break; + } + } + } else { + if (verbose_flag) + printf("All facets connected. No nearby check necessary.\n"); + } + } + + if(remove_unconnected_flag || fixall_flag || fill_holes_flag) { + if(stl->stats.connected_facets_3_edge < stl->stats.number_of_facets) { + if (verbose_flag) + printf("Removing unconnected facets...\n"); + stl_remove_unconnected_facets(stl); + } else + if (verbose_flag) + printf("No unconnected need to be removed.\n"); + } + + if(fill_holes_flag || fixall_flag) { + if(stl->stats.connected_facets_3_edge < stl->stats.number_of_facets) { + if (verbose_flag) + printf("Filling holes...\n"); + stl_fill_holes(stl); + } else + if (verbose_flag) + printf("No holes need to be filled.\n"); + } + + if(reverse_all_flag) { + if (verbose_flag) + printf("Reversing all facets...\n"); + stl_reverse_all_facets(stl); + } + + if(normal_directions_flag || fixall_flag) { + if (verbose_flag) + printf("Checking normal directions...\n"); + stl_fix_normal_directions(stl); + } + + if(normal_values_flag || fixall_flag) { + if (verbose_flag) + printf("Checking normal values...\n"); + stl_fix_normal_values(stl); + } + + /* Always calculate the volume. It shouldn't take too long */ + if (verbose_flag) + printf("Calculating volume...\n"); + stl_calculate_volume(stl); + + if(exact_flag) { + if (verbose_flag) + printf("Verifying neighbors...\n"); + stl_verify_neighbors(stl); + } +} |