This commit is contained in:
bubnikv 2018-08-22 15:03:35 +02:00
parent ac72cd779f
commit 6829704475
16 changed files with 539 additions and 908 deletions

View file

@ -25,11 +25,11 @@
#include <string.h>
#include <math.h>
#include <boost/detail/endian.hpp>
#include "stl.h"
static void stl_match_neighbors_exact(stl_file *stl,
stl_hash_edge *edge_a, stl_hash_edge *edge_b);
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,
@ -43,7 +43,6 @@ static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge,
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_get_hash_for_edge(int M, stl_hash_edge *edge);
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);
@ -82,37 +81,20 @@ stl_check_facets_exact(stl_file *stl) {
for(i = 0; i < stl->stats.number_of_facets; i++) {
facet = stl->facet_start[i];
// Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
// When using a memcmp on raw floats, those numbers report to be different.
// Unify all +0 and -0 to +0 to make the floats equal under memcmp.
{
uint32_t *f = (uint32_t*)&facet;
for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
if (*f == 0x80000000)
// Negative zero, switch to positive zero.
*f = 0;
}
/* If any two of the three vertices are found to be exactally the same, call them degenerate and remove the facet. */
if( !memcmp(&facet.vertex[0], &facet.vertex[1],
sizeof(stl_vertex))
|| !memcmp(&facet.vertex[1], &facet.vertex[2],
sizeof(stl_vertex))
|| !memcmp(&facet.vertex[0], &facet.vertex[2],
sizeof(stl_vertex))) {
// 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--;
-- 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_match_neighbors_exact);
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);
@ -131,28 +113,33 @@ stl_load_edge_exact(stl_file *stl, stl_hash_edge *edge,
if (stl->error) return;
{
float diff_x = ABS(a->x - b->x);
float diff_y = ABS(a->y - b->y);
float diff_z = ABS(a->z - b->z);
float max_diff = STL_MAX(diff_x, diff_y);
max_diff = STL_MAX(diff_z, max_diff);
stl->stats.shortest_edge = STL_MIN(max_diff, stl->stats.shortest_edge);
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 ((a->x != b->x) ?
(a->x < b->x) :
((a->y != b->y) ?
(a->y < b->y) :
(a->z < b->z))) {
memcpy(&edge->key[0], a, sizeof(stl_vertex));
memcpy(&edge->key[3], b, sizeof(stl_vertex));
if (stl_vertex_lower(*a, *b)) {
} else {
memcpy(&edge->key[0], b, sizeof(stl_vertex));
memcpy(&edge->key[3], a, sizeof(stl_vertex));
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
@ -188,21 +175,17 @@ stl_initialize_facet_check_exact(stl_file *stl) {
}
}
static void
insert_hash_edge(stl_file *stl, stl_hash_edge edge,
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)) {
stl_hash_edge *link;
stl_hash_edge *new_edge;
stl_hash_edge *temp;
int chain_number;
stl_hash_edge *edge_a, stl_hash_edge *edge_b))
{
if (stl->error) return;
chain_number = stl_get_hash_for_edge(stl->M, &edge);
link = stl->heads[chain_number];
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));
@ -252,30 +235,17 @@ insert_hash_edge(stl_file *stl, stl_hash_edge edge,
}
}
static int
stl_get_hash_for_edge(int M, stl_hash_edge *edge) {
return ((edge->key[0] / 23 + edge->key[1] / 19 + edge->key[2] / 17
+ edge->key[3] /13 + edge->key[4] / 11 + edge->key[5] / 7 ) % M);
// 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);
}
static int
stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
if(edge_a->facet_number == edge_b->facet_number) {
return 1; /* Don't match edges of the same facet */
} else {
return memcmp(edge_a, edge_b, SIZEOF_EDGE_SORT);
}
}
void
stl_check_facets_nearby(stl_file *stl, float tolerance) {
stl_hash_edge edge[3];
stl_facet facet;
int i;
int j;
if (stl->error) return;
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)
@ -286,27 +256,19 @@ stl_check_facets_nearby(stl_file *stl, float tolerance) {
stl_initialize_facet_check_nearby(stl);
for(i = 0; i < stl->stats.number_of_facets; i++) {
facet = stl->facet_start[i];
// Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
// When using a memcmp on raw floats, those numbers report to be different.
// Unify all +0 and -0 to +0 to make the floats equal under memcmp.
{
uint32_t *f = (uint32_t*)&facet;
for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
if (*f == 0x80000000)
// Negative zero, switch to positive zero.
*f = 0;
}
for(j = 0; j < 3; j++) {
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) {
edge[j].facet_number = i;
edge[j].which_edge = j;
if(stl_load_edge_nearby(stl, &edge[j], &facet.vertex[j],
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[j], stl_match_neighbors_nearby);
insert_hash_edge(stl, edge, stl_match_neighbors_nearby);
}
}
}
@ -315,27 +277,17 @@ stl_check_facets_nearby(stl_file *stl, float tolerance) {
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) {
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.
uint32_t vertex1[3] = {
(uint32_t)((a->x - stl->stats.min.x) / tolerance),
(uint32_t)((a->y - stl->stats.min.y) / tolerance),
(uint32_t)((a->z - stl->stats.min.z) / tolerance)
};
uint32_t vertex2[3] = {
(uint32_t)((b->x - stl->stats.min.x) / tolerance),
(uint32_t)((b->y - stl->stats.min.y) / tolerance),
(uint32_t)((b->z - stl->stats.min.z) / tolerance)
};
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[0] == vertex2[0])
&& (vertex1[1] == vertex2[1])
&& (vertex1[2] == vertex2[2])) {
/* Both vertices hash to the same value */
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.
@ -344,30 +296,27 @@ stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge,
((vertex1[1] != vertex2[1]) ?
(vertex1[1] < vertex2[1]) :
(vertex1[2] < vertex2[2]))) {
memcpy(&edge->key[0], vertex1, sizeof(stl_vertex));
memcpy(&edge->key[3], vertex2, sizeof(stl_vertex));
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, sizeof(stl_vertex));
memcpy(&edge->key[3], vertex1, sizeof(stl_vertex));
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) {
int i;
stl_hash_edge *temp;
if (stl->error) return;
static void stl_free_edges(stl_file *stl)
{
if (stl->error)
return;
if(stl->stats.malloced != stl->stats.freed) {
for(i = 0; i < stl->M; i++) {
for(temp = stl->heads[i]; stl->heads[i] != stl->tail;
temp = stl->heads[i]) {
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++;
++ stl->stats.freed;
}
}
}
@ -375,8 +324,8 @@ stl_free_edges(stl_file *stl) {
free(stl->tail);
}
static void
stl_initialize_facet_check_nearby(stl_file *stl) {
static void stl_initialize_facet_check_nearby(stl_file *stl)
{
int i;
if (stl->error) return;
@ -467,16 +416,8 @@ stl_record_neighbors(stl_file *stl,
}
}
static void
stl_match_neighbors_exact(stl_file *stl,
stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
if (stl->error) return;
stl_record_neighbors(stl, edge_a, edge_b);
}
static void
stl_match_neighbors_nearby(stl_file *stl,
stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
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;
@ -517,9 +458,7 @@ stl_match_neighbors_nearby(stl_file *stl,
}
static void
stl_change_vertices(stl_file *stl, int facet_num, int vnot,
stl_vertex new_vertex) {
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;
@ -551,30 +490,30 @@ stl_change_vertices(stl_file *stl, int facet_num, int vnot,
}
}
#if 0
if (stl->facet_start[facet_num].vertex[pivot_vertex].x == new_vertex.x &&
stl->facet_start[facet_num].vertex[pivot_vertex].y == new_vertex.y &&
stl->facet_start[facet_num].vertex[pivot_vertex].z == new_vertex.z)
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.x, new_vertex.y, new_vertex.z);
new_vertex(0), new_vertex(1), new_vertex(2));
else {
if (stl->facet_start[facet_num].vertex[pivot_vertex].x != new_vertex.x)
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].x,
*reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].x),
new_vertex.x,
*reinterpret_cast<const int*>(&new_vertex.x));
if (stl->facet_start[facet_num].vertex[pivot_vertex].y != new_vertex.y)
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].y,
*reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].y),
new_vertex.y,
*reinterpret_cast<const int*>(&new_vertex.y));
if (stl->facet_start[facet_num].vertex[pivot_vertex].z != new_vertex.z)
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].z,
*reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].z),
new_vertex.z,
*reinterpret_cast<const int*>(&new_vertex.z));
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;
@ -595,7 +534,6 @@ Try using a smaller tolerance or don't do a nearby check\n");
}
}
static void
stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
stl_hash_edge *edge_b, int *facet1, int *vertex1,
@ -622,11 +560,10 @@ stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
v1b = (edge_b->which_edge + 1) % 3;
}
/* Of the first pair, which vertex, if any, should be changed */
if(!memcmp(&stl->facet_start[edge_a->facet_number].vertex[v1a],
&stl->facet_start[edge_b->facet_number].vertex[v1b],
sizeof(stl_vertex))) {
/* These facets are already equal. No need to change. */
// 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)
@ -644,10 +581,9 @@ stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
}
/* Of the second pair, which vertex, if any, should be changed */
if(!memcmp(&stl->facet_start[edge_a->facet_number].vertex[v2a],
&stl->facet_start[edge_b->facet_number].vertex[v2b],
sizeof(stl_vertex))) {
/* These facets are already equal. No need to change. */
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)
@ -718,40 +654,35 @@ in stl_remove_facet: neighbor = %d numfacets = %d this is wrong\n",
}
}
void
stl_remove_unconnected_facets(stl_file *stl) {
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;
int i;
if (stl->error) return;
/* remove degenerate facets */
for(i = 0; i < stl->stats.number_of_facets; i++) {
if( !memcmp(&stl->facet_start[i].vertex[0],
&stl->facet_start[i].vertex[1], sizeof(stl_vertex))
|| !memcmp(&stl->facet_start[i].vertex[1],
&stl->facet_start[i].vertex[2], sizeof(stl_vertex))
|| !memcmp(&stl->facet_start[i].vertex[0],
&stl->facet_start[i].vertex[2], sizeof(stl_vertex))) {
// 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(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. */
// 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--;
-- i;
}
}
}
@ -771,30 +702,24 @@ stl_remove_degenerate(stl_file *stl, int facet) {
if (stl->error) return;
if( !memcmp(&stl->facet_start[facet].vertex[0],
&stl->facet_start[facet].vertex[1], sizeof(stl_vertex))
&& !memcmp(&stl->facet_start[facet].vertex[1],
&stl->facet_start[facet].vertex[2], sizeof(stl_vertex))) {
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(!memcmp(&stl->facet_start[facet].vertex[0],
&stl->facet_start[facet].vertex[1], sizeof(stl_vertex))) {
if (stl->facet_start[facet].vertex[0] == stl->facet_start[facet].vertex[1]) {
edge1 = 1;
edge2 = 2;
edge3 = 0;
} else if(!memcmp(&stl->facet_start[facet].vertex[1],
&stl->facet_start[facet].vertex[2], sizeof(stl_vertex))) {
} else if (stl->facet_start[facet].vertex[1] == stl->facet_start[facet].vertex[2]) {
edge1 = 0;
edge2 = 2;
edge3 = 1;
} else if(!memcmp(&stl->facet_start[facet].vertex[2],
&stl->facet_start[facet].vertex[0], sizeof(stl_vertex))) {
} else if (stl->facet_start[facet].vertex[2] == stl->facet_start[facet].vertex[0]) {
edge1 = 0;
edge2 = 1;
edge3 = 2;
@ -883,7 +808,7 @@ stl_fill_holes(stl_file *stl) {
stl_load_edge_exact(stl, &edge, &facet.vertex[j],
&facet.vertex[(j + 1) % 3]);
insert_hash_edge(stl, edge, stl_match_neighbors_exact);
insert_hash_edge(stl, edge, stl_record_neighbors);
}
}
@ -939,7 +864,7 @@ stl_fill_holes(stl_file *stl) {
stl_load_edge_exact(stl, &edge, &new_facet.vertex[k],
&new_facet.vertex[(k + 1) % 3]);
insert_hash_edge(stl, edge, stl_match_neighbors_exact);
insert_hash_edge(stl, edge, stl_record_neighbors);
}
break;
} else {
@ -977,9 +902,7 @@ stl_add_facet(stl_file *stl, stl_facet *new_facet) {
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.x = 0.0;
stl->facet_start[stl->stats.number_of_facets].normal.y = 0.0;
stl->facet_start[stl->stats.number_of_facets].normal.z = 0.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;

View file

@ -27,12 +27,6 @@
#include "stl.h"
static void stl_reverse_vector(float v[]) {
v[0] *= -1;
v[1] *= -1;
v[2] *= -1;
}
static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag);
static void
@ -228,102 +222,52 @@ static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_
/* Returns 2 if the normal is not within tolerance and backwards */
/* Returns 4 if the status is unknown. */
float normal[3];
float test_norm[3];
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();
if( (ABS(normal[0] - facet->normal.x) < 0.001)
&& (ABS(normal[1] - facet->normal.y) < 0.001)
&& (ABS(normal[2] - facet->normal.z) < 0.001)) {
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.x = normal[0];
facet->normal.y = normal[1];
facet->normal.z = normal[2];
facet->normal = normal;
return 0;
}
test_norm[0] = facet->normal.x;
test_norm[1] = facet->normal.y;
test_norm[2] = facet->normal.z;
stl_normal test_norm = facet->normal;
stl_normalize_vector(test_norm);
if( (ABS(normal[0] - test_norm[0]) < 0.001)
&& (ABS(normal[1] - test_norm[1]) < 0.001)
&& (ABS(normal[2] - test_norm[2]) < 0.001)) {
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.x = normal[0];
facet->normal.y = normal[1];
facet->normal.z = normal[2];
facet->normal = normal;
stl->stats.normals_fixed += 1;
}
return 1;
}
stl_reverse_vector(test_norm);
if( (ABS(normal[0] - test_norm[0]) < 0.001)
&& (ABS(normal[1] - test_norm[1]) < 0.001)
&& (ABS(normal[2] - test_norm[2]) < 0.001)) {
/* Facet is backwards. */
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.x = normal[0];
facet->normal.y = normal[1];
facet->normal.z = normal[2];
facet->normal = normal;
stl->stats.normals_fixed += 1;
}
return 2;
}
if(normal_fix_flag) {
facet->normal.x = normal[0];
facet->normal.y = normal[1];
facet->normal.z = normal[2];
facet->normal = normal;
stl->stats.normals_fixed += 1;
}
return 4;
}
void stl_calculate_normal(float normal[], stl_facet *facet) {
float v1[3] = {
facet->vertex[1].x - facet->vertex[0].x,
facet->vertex[1].y - facet->vertex[0].y,
facet->vertex[1].z - facet->vertex[0].z
};
float v2[3] = {
facet->vertex[2].x - facet->vertex[0].x,
facet->vertex[2].y - facet->vertex[0].y,
facet->vertex[2].z - facet->vertex[0].z
};
normal[0] = (float)((double)v1[1] * (double)v2[2]) - ((double)v1[2] * (double)v2[1]);
normal[1] = (float)((double)v1[2] * (double)v2[0]) - ((double)v1[0] * (double)v2[2]);
normal[2] = (float)((double)v1[0] * (double)v2[1]) - ((double)v1[1] * (double)v2[0]);
}
void stl_normalize_vector(float v[]) {
double length;
double factor;
float min_normal_length;
length = sqrt((double)v[0] * (double)v[0] + (double)v[1] * (double)v[1] + (double)v[2] * (double)v[2]);
min_normal_length = 0.000000000001;
if(length < min_normal_length) {
v[0] = 0.0;
v[1] = 0.0;
v[2] = 0.0;
return;
}
factor = 1.0 / length;
v[0] *= factor;
v[1] *= factor;
v[2] *= factor;
}
void
stl_fix_normal_values(stl_file *stl) {
void stl_fix_normal_values(stl_file *stl) {
int i;
if (stl->error) return;
@ -333,20 +277,16 @@ stl_fix_normal_values(stl_file *stl) {
}
}
void
stl_reverse_all_facets(stl_file *stl) {
int i;
float normal[3];
void stl_reverse_all_facets(stl_file *stl)
{
if (stl->error)
return;
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
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.x = normal[0];
stl->facet_start[i].normal.y = normal[1];
stl->facet_start[i].normal.z = normal[2];
stl->facet_start[i].normal = normal;
}
}

View file

@ -169,7 +169,7 @@ stl_write_off(stl_file *stl, char *file) {
for(i = 0; i < stl->stats.shared_vertices; i++) {
fprintf(fp, "\t%f %f %f\n",
stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
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],
@ -216,10 +216,10 @@ stl_write_vrml(stl_file *stl, char *file) {
for(i = 0; i < (stl->stats.shared_vertices - 1); i++) {
fprintf(fp, "\t\t\t\t%f %f %f,\n",
stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
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].x, stl->v_shared[i].y, stl->v_shared[i].z);
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");
@ -254,7 +254,7 @@ void stl_write_obj (stl_file *stl, char *file) {
}
for (i = 0; i < stl->stats.shared_vertices; i++) {
fprintf(fp, "v %f %f %f\n", stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
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);

View file

@ -27,9 +27,7 @@
#include <stdint.h>
#include <stddef.h>
#define STL_MAX(A,B) ((A)>(B)? (A):(B))
#define STL_MIN(A,B) ((A)<(B)? (A):(B))
#define ABS(X) ((X) < 0 ? -(X) : (X))
#include <Eigen/Geometry>
// Size of the binary STL header, free form.
#define LABEL_SIZE 80
@ -39,31 +37,16 @@
#define HEADER_SIZE 84
#define STL_MIN_FILE_SIZE 284
#define ASCII_LINES_PER_FACET 7
// Comparing an edge by memcmp, 2x3x4 bytes = 24
#define SIZEOF_EDGE_SORT 24
typedef struct {
float x;
float y;
float z;
} stl_vertex;
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");
typedef struct {
float x;
float y;
float z;
} stl_normal;
static_assert(sizeof(stl_normal) == 12, "size of stl_normal incorrect");
typedef char stl_extra[2];
typedef struct {
stl_normal normal;
stl_vertex vertex[3];
stl_extra extra;
char extra[2];
} stl_facet;
#define SIZEOF_STL_FACET 50
@ -81,8 +64,12 @@ typedef struct {
} stl_edge;
typedef struct stl_hash_edge {
// Key of a hash edge: 2x binary copy of a floating point vertex.
uint32_t key[6];
// 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.
@ -91,8 +78,6 @@ typedef struct stl_hash_edge {
struct stl_hash_edge *next;
} stl_hash_edge;
static_assert(offsetof(stl_hash_edge, facet_number) == SIZEOF_EDGE_SORT, "size of stl_hash_edge.key incorrect");
typedef struct {
// Index of a neighbor facet.
int neighbor[3];
@ -179,8 +164,8 @@ 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, float versor[3]);
extern void stl_scale(stl_file *stl, float factor);
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);
@ -195,8 +180,20 @@ 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);
extern void stl_calculate_normal(float normal[], stl_facet *facet);
extern void stl_normalize_vector(float v[]);
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);
@ -204,8 +201,8 @@ extern void stl_repair(stl_file *stl, int fixall_flag, int exact_flag, int toler
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, int first);
extern void stl_facet_stats(stl_file *stl, stl_facet facet, int first);
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);

View file

@ -44,9 +44,9 @@ stl_print_edges(stl_file *stl, FILE *file) {
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.x, stl->edge_start[i].p1.y,
stl->edge_start[i].p1.z, stl->edge_start[i].p2.x,
stl->edge_start[i].p2.y, stl->edge_start[i].p2.z);
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));
}
}
@ -75,11 +75,11 @@ File type : ASCII STL file\n");
Header : %s\n", stl->stats.header);
fprintf(file, "============== Size ==============\n");
fprintf(file, "Min X = % f, Max X = % f\n",
stl->stats.min.x, stl->stats.max.x);
stl->stats.min(0), stl->stats.max(0));
fprintf(file, "Min Y = % f, Max Y = % f\n",
stl->stats.min.y, stl->stats.max.y);
stl->stats.min(1), stl->stats.max(1));
fprintf(file, "Min Z = % f, Max Z = % f\n",
stl->stats.min.z, stl->stats.max.z);
stl->stats.min(2), stl->stats.max(2));
fprintf(file, "\
========= Facet Status ========== Original ============ Final ====\n");
@ -149,18 +149,18 @@ stl_write_ascii(stl_file *stl, const char *file, const char *label) {
for(i = 0; i < stl->stats.number_of_facets; i++) {
fprintf(fp, " facet normal % .8E % .8E % .8E\n",
stl->facet_start[i].normal.x, stl->facet_start[i].normal.y,
stl->facet_start[i].normal.z);
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].x, stl->facet_start[i].vertex[0].y,
stl->facet_start[i].vertex[0].z);
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].x, stl->facet_start[i].vertex[1].y,
stl->facet_start[i].vertex[1].z);
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].x, stl->facet_start[i].vertex[2].y,
stl->facet_start[i].vertex[2].z);
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");
}
@ -264,9 +264,9 @@ 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].x,
stl->facet_start[facet].vertex[vertex].y,
stl->facet_start[facet].vertex[vertex].z);
stl->facet_start[facet].vertex[vertex](0),
stl->facet_start[facet].vertex[vertex](1),
stl->facet_start[facet].vertex[vertex](2));
}
void
@ -309,10 +309,10 @@ stl_write_quad_object(stl_file *stl, char *file) {
int i;
int j;
char *error_msg;
stl_vertex connect_color;
stl_vertex uncon_1_color;
stl_vertex uncon_2_color;
stl_vertex uncon_3_color;
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;
@ -330,19 +330,6 @@ stl_write_quad_object(stl_file *stl, char *file) {
return;
}
connect_color.x = 0.0;
connect_color.y = 0.0;
connect_color.z = 1.0;
uncon_1_color.x = 0.0;
uncon_1_color.y = 1.0;
uncon_1_color.z = 0.0;
uncon_2_color.x = 1.0;
uncon_2_color.y = 1.0;
uncon_2_color.z = 1.0;
uncon_3_color.x = 1.0;
uncon_3_color.y = 0.0;
uncon_3_color.z = 0.0;
fprintf(fp, "CQUAD\n");
for(i = 0; i < stl->stats.number_of_facets; i++) {
j = ((stl->neighbors_start[i].neighbor[0] == -1) +
@ -358,21 +345,21 @@ stl_write_quad_object(stl_file *stl, char *file) {
color = uncon_3_color;
}
fprintf(fp, "%f %f %f %1.1f %1.1f %1.1f 1\n",
stl->facet_start[i].vertex[0].x,
stl->facet_start[i].vertex[0].y,
stl->facet_start[i].vertex[0].z, color.x, color.y, color.z);
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].x,
stl->facet_start[i].vertex[1].y,
stl->facet_start[i].vertex[1].z, color.x, color.y, color.z);
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].x,
stl->facet_start[i].vertex[2].y,
stl->facet_start[i].vertex[2].z, color.x, color.y, color.z);
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].x,
stl->facet_start[i].vertex[2].y,
stl->facet_start[i].vertex[2].z, color.x, color.y, color.z);
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);
}
@ -409,17 +396,17 @@ stl_write_dxf(stl_file *stl, char *file, char *label) {
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].x, stl->facet_start[i].vertex[0].y,
stl->facet_start[i].vertex[0].z);
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].x, stl->facet_start[i].vertex[1].y,
stl->facet_start[i].vertex[1].z);
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].x, stl->facet_start[i].vertex[2].y,
stl->facet_start[i].vertex[2].z);
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].x, stl->facet_start[i].vertex[2].y,
stl->facet_start[i].vertex[2].z);
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");

View file

@ -40,7 +40,7 @@ stl_open(stl_file *stl, const char *file) {
stl_initialize(stl);
stl_count_facets(stl, file);
stl_allocate(stl);
stl_read(stl, 0, 1);
stl_read(stl, 0, true);
if (!stl->error) fclose(stl->fp);
}
@ -227,7 +227,7 @@ stl_open_merge(stl_file *stl, char *file_to_merge) {
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, 0);
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: */
@ -255,8 +255,7 @@ stl_reallocate(stl_file *stl) {
/* 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, int first) {
void stl_read(stl_file *stl, int first_facet, bool first) {
stl_facet facet;
int i;
@ -294,11 +293,11 @@ stl_read(stl_file *stl, int first_facet, int first) {
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].x, &facet.vertex[0].y, &facet.vertex[0].z);
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].x, &facet.vertex[1].y, &facet.vertex[1].z);
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].x, &facet.vertex[2].y, &facet.vertex[2].z);
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);
@ -311,9 +310,9 @@ stl_read(stl_file *stl, int first_facet, int first) {
}
// 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.x) != 1 ||
sscanf(normal_buf[1], "%f", &facet.normal.y) != 1 ||
sscanf(normal_buf[2], "%f", &facet.normal.z) != 1) {
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));
@ -326,104 +325,45 @@ stl_read(stl_file *stl, int first_facet, int first) {
// It may be worth to round these numbers to zero during loading to reduce the number of errors reported
// during the STL import.
for (size_t j = 0; j < 3; ++ j) {
if (facet.vertex[j].x > -1e-12f && facet.vertex[j].x < 1e-12f)
printf("stl_read: facet %d.x = %e\r\n", j, facet.vertex[j].x);
if (facet.vertex[j].y > -1e-12f && facet.vertex[j].y < 1e-12f)
printf("stl_read: facet %d.y = %e\r\n", j, facet.vertex[j].y);
if (facet.vertex[j].z > -1e-12f && facet.vertex[j].z < 1e-12f)
printf("stl_read: facet %d.z = %e\r\n", j, facet.vertex[j].z);
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
#if 1
{
// Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
// When using a memcmp on raw floats, those numbers report to be different.
// Unify all +0 and -0 to +0 to make the floats equal under memcmp.
uint32_t *f = (uint32_t*)&facet;
for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
if (*f == 0x80000000)
// Negative zero, switch to positive zero.
*f = 0;
}
#else
{
// 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. Round them to zero.
float *f = (float*)&facet;
for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
if (*f > -1e-12f && *f < 1e-12f)
// Negative zero, switch to positive zero.
*f = 0;
}
#endif
/* Write the facet into memory. */
memcpy(stl->facet_start+i, &facet, SIZEOF_STL_FACET);
stl->facet_start[i] = facet;
stl_facet_stats(stl, facet, first);
first = 0;
}
stl->stats.size.x = stl->stats.max.x - stl->stats.min.x;
stl->stats.size.y = stl->stats.max.y - stl->stats.min.y;
stl->stats.size.z = stl->stats.max.z - stl->stats.min.z;
stl->stats.bounding_diameter = sqrt(
stl->stats.size.x * stl->stats.size.x +
stl->stats.size.y * stl->stats.size.y +
stl->stats.size.z * stl->stats.size.z
);
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, int first) {
float diff_x;
float diff_y;
float diff_z;
float max_diff;
void stl_facet_stats(stl_file *stl, stl_facet facet, bool &first)
{
if (stl->error)
return;
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
/* while we are going through all of the facets, let's find the */
/* maximum and minimum values for x, y, and z */
/* Initialize the max and min values the first time through*/
if (first) {
stl->stats.max.x = facet.vertex[0].x;
stl->stats.min.x = facet.vertex[0].x;
stl->stats.max.y = facet.vertex[0].y;
stl->stats.min.y = facet.vertex[0].y;
stl->stats.max.z = facet.vertex[0].z;
stl->stats.min.z = facet.vertex[0].z;
diff_x = ABS(facet.vertex[0].x - facet.vertex[1].x);
diff_y = ABS(facet.vertex[0].y - facet.vertex[1].y);
diff_z = ABS(facet.vertex[0].z - facet.vertex[1].z);
max_diff = STL_MAX(diff_x, diff_y);
max_diff = STL_MAX(diff_z, max_diff);
stl->stats.shortest_edge = max_diff;
first = 0;
// 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 */
stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[0].x);
stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[0].x);
stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[0].y);
stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[0].y);
stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[0].z);
stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[0].z);
stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[1].x);
stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[1].x);
stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[1].y);
stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[1].y);
stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[1].z);
stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[1].z);
stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[2].x);
stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[2].x);
stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[2].y);
stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[2].y);
stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[2].z);
stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[2].z);
// 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.cwiseMin(facet.vertex[i]);
}
}
void

View file

@ -62,7 +62,7 @@ stl_verify_neighbors(stl_file *stl) {
edge_b.p1 = stl->facet_start[neighbor].vertex[(vnot + 1) % 3];
edge_b.p2 = stl->facet_start[neighbor].vertex[(vnot + 2) % 3];
}
if(memcmp(&edge_a, &edge_b, SIZEOF_EDGE_SORT) != 0) {
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);
@ -73,114 +73,67 @@ stl_verify_neighbors(stl_file *stl) {
}
}
void
stl_translate(stl_file *stl, float x, float y, float z) {
int i;
int j;
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].x -= (stl->stats.min.x - x);
stl->facet_start[i].vertex[j].y -= (stl->stats.min.y - y);
stl->facet_start[i].vertex[j].z -= (stl->stats.min.z - z);
}
}
stl->stats.max.x -= (stl->stats.min.x - x);
stl->stats.max.y -= (stl->stats.min.y - y);
stl->stats.max.z -= (stl->stats.min.z - z);
stl->stats.min.x = x;
stl->stats.min.y = y;
stl->stats.min.z = z;
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) {
int i;
int j;
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].x += x;
stl->facet_start[i].vertex[j].y += y;
stl->facet_start[i].vertex[j].z += z;
}
}
stl->stats.min.x += x;
stl->stats.min.y += y;
stl->stats.min.z += z;
stl->stats.max.x += x;
stl->stats.max.y += y;
stl->stats.max.z += z;
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, float versor[3]) {
int i;
int j;
if (stl->error) return;
/* scale extents */
stl->stats.min.x *= versor[0];
stl->stats.min.y *= versor[1];
stl->stats.min.z *= versor[2];
stl->stats.max.x *= versor[0];
stl->stats.max.y *= versor[1];
stl->stats.max.z *= versor[2];
/* scale size */
stl->stats.size.x *= versor[0];
stl->stats.size.y *= versor[1];
stl->stats.size.z *= versor[2];
/* scale volume */
if (stl->stats.volume > 0.0) {
stl->stats.volume *= (versor[0] * versor[1] * versor[2]);
}
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].x *= versor[0];
stl->facet_start[i].vertex[j].y *= versor[1];
stl->facet_start[i].vertex[j].z *= versor[2];
}
}
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);
}
void
stl_scale(stl_file *stl, float factor) {
float versor[3];
if (stl->error) return;
versor[0] = factor;
versor[1] = factor;
versor[2] = factor;
stl_scale_versor(stl, versor);
}
static void calculate_normals(stl_file *stl) {
float normal[3];
if (stl->error) return;
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.x = normal[0];
stl->facet_start[i].normal.y = normal[1];
stl->facet_start[i].normal.z = normal[2];
stl->facet_start[i].normal = normal;
}
}
@ -193,9 +146,9 @@ void stl_transform(stl_file *stl, float *trafo3x4) {
for (i_vertex = 0; i_vertex < 3; ++ i_vertex) {
stl_vertex &v_dst = vertices[i_vertex];
stl_vertex v_src = v_dst;
v_dst.x = trafo3x4[0] * v_src.x + trafo3x4[1] * v_src.y + trafo3x4[2] * v_src.z + trafo3x4[3];
v_dst.y = trafo3x4[4] * v_src.x + trafo3x4[5] * v_src.y + trafo3x4[6] * v_src.z + trafo3x4[7];
v_dst.z = trafo3x4[8] * v_src.x + trafo3x4[9] * v_src.y + trafo3x4[10] * v_src.z + trafo3x4[11];
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);
@ -214,8 +167,8 @@ stl_rotate_x(stl_file *stl, float angle) {
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl_rotate(&stl->facet_start[i].vertex[j].y,
&stl->facet_start[i].vertex[j].z, c, s);
stl_rotate(&stl->facet_start[i].vertex[j](1),
&stl->facet_start[i].vertex[j](2), c, s);
}
}
stl_get_size(stl);
@ -234,8 +187,8 @@ stl_rotate_y(stl_file *stl, float angle) {
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl_rotate(&stl->facet_start[i].vertex[j].z,
&stl->facet_start[i].vertex[j].x, c, s);
stl_rotate(&stl->facet_start[i].vertex[j](2),
&stl->facet_start[i].vertex[j](0), c, s);
}
}
stl_get_size(stl);
@ -254,8 +207,8 @@ stl_rotate_z(stl_file *stl, float angle) {
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl_rotate(&stl->facet_start[i].vertex[j].x,
&stl->facet_start[i].vertex[j].y, c, s);
stl_rotate(&stl->facet_start[i].vertex[j](0),
&stl->facet_start[i].vertex[j](1), c, s);
}
}
stl_get_size(stl);
@ -272,142 +225,98 @@ stl_rotate(float *x, float *y, const double c, const double s) {
*y = float(s * xold + c * yold);
}
extern void
stl_get_size(stl_file *stl) {
int i;
int j;
if (stl->error) return;
if (stl->stats.number_of_facets == 0) return;
stl->stats.min.x = stl->facet_start[0].vertex[0].x;
stl->stats.min.y = stl->facet_start[0].vertex[0].y;
stl->stats.min.z = stl->facet_start[0].vertex[0].z;
stl->stats.max.x = stl->facet_start[0].vertex[0].x;
stl->stats.max.y = stl->facet_start[0].vertex[0].y;
stl->stats.max.z = stl->facet_start[0].vertex[0].z;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->stats.min.x = STL_MIN(stl->stats.min.x,
stl->facet_start[i].vertex[j].x);
stl->stats.min.y = STL_MIN(stl->stats.min.y,
stl->facet_start[i].vertex[j].y);
stl->stats.min.z = STL_MIN(stl->stats.min.z,
stl->facet_start[i].vertex[j].z);
stl->stats.max.x = STL_MAX(stl->stats.max.x,
stl->facet_start[i].vertex[j].x);
stl->stats.max.y = STL_MAX(stl->stats.max.y,
stl->facet_start[i].vertex[j].y);
stl->stats.max.z = STL_MAX(stl->stats.max.z,
stl->facet_start[i].vertex[j].z);
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.x = stl->stats.max.x - stl->stats.min.x;
stl->stats.size.y = stl->stats.max.y - stl->stats.min.y;
stl->stats.size.z = stl->stats.max.z - stl->stats.min.z;
stl->stats.bounding_diameter = sqrt(
stl->stats.size.x * stl->stats.size.x +
stl->stats.size.y * stl->stats.size.y +
stl->stats.size.z * stl->stats.size.z
);
stl->stats.size = stl->stats.max - stl->stats.min;
stl->stats.bounding_diameter = stl->stats.size.norm();
}
void
stl_mirror_xy(stl_file *stl) {
int i;
int j;
float temp_size;
void stl_mirror_xy(stl_file *stl)
{
if (stl->error)
return;
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].z *= -1.0;
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;
}
}
temp_size = stl->stats.min.z;
stl->stats.min.z = stl->stats.max.z;
stl->stats.max.z = temp_size;
stl->stats.min.z *= -1.0;
stl->stats.max.z *= -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) {
int i;
int j;
float temp_size;
void stl_mirror_yz(stl_file *stl)
{
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].x *= -1.0;
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;
}
}
temp_size = stl->stats.min.x;
stl->stats.min.x = stl->stats.max.x;
stl->stats.max.x = temp_size;
stl->stats.min.x *= -1.0;
stl->stats.max.x *= -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) {
int i;
int j;
float temp_size;
void stl_mirror_xz(stl_file *stl)
{
if (stl->error)
return;
if (stl->error) return;
for(i = 0; i < stl->stats.number_of_facets; i++) {
for(j = 0; j < 3; j++) {
stl->facet_start[i].vertex[j].y *= -1.0;
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;
}
}
temp_size = stl->stats.min.y;
stl->stats.min.y = stl->stats.max.y;
stl->stats.max.y = temp_size;
stl->stats.min.y *= -1.0;
stl->stats.max.y *= -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) {
stl_vertex p0;
stl_vertex p;
stl_normal n;
float height;
float area;
float volume = 0.0;
static float get_volume(stl_file *stl)
{
if (stl->error)
return 0;
if (stl->error) return 0;
/* Choose a point, any point as the reference */
p0.x = stl->facet_start[0].vertex[0].x;
p0.y = stl->facet_start[0].vertex[0].y;
p0.z = stl->facet_start[0].vertex[0].z;
for(uint32_t i = 0; i < stl->stats.number_of_facets; i++) {
p.x = stl->facet_start[i].vertex[0].x - p0.x;
p.y = stl->facet_start[i].vertex[0].y - p0.y;
p.z = stl->facet_start[i].vertex[0].z - p0.z;
/* Do dot product to get distance from point to plane */
n = stl->facet_start[i].normal;
height = (n.x * p.x) + (n.y * p.y) + (n.z * p.z);
area = get_area(&stl->facet_start[i]);
// 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) {
void stl_calculate_volume(stl_file *stl)
{
if (stl->error) return;
stl->stats.volume = get_volume(stl);
if(stl->stats.volume < 0.0) {
@ -416,35 +325,32 @@ void stl_calculate_volume(stl_file *stl) {
}
}
static float get_area(stl_facet *facet) {
double cross[3][3];
float sum[3];
float n[3];
float area;
int i;
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) */
for(i = 0; i < 3; i++) {
cross[i][0]=(((double)facet->vertex[i].y * (double)facet->vertex[(i + 1) % 3].z) -
((double)facet->vertex[i].z * (double)facet->vertex[(i + 1) % 3].y));
cross[i][1]=(((double)facet->vertex[i].z * (double)facet->vertex[(i + 1) % 3].x) -
((double)facet->vertex[i].x * (double)facet->vertex[(i + 1) % 3].z));
cross[i][2]=(((double)facet->vertex[i].x * (double)facet->vertex[(i + 1) % 3].y) -
((double)facet->vertex[i].y * (double)facet->vertex[(i + 1) % 3].x));
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)));
}
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];
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 */
// 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);
area = 0.5 * (n[0] * sum[0] + n[1] * sum[1] + n[2] * sum[2]);
return area;
return 0.5f * n.dot(sum);
}
void stl_repair(stl_file *stl,