From 519ed91c68b78ad6c336669fe9e5f84dabbbfe60 Mon Sep 17 00:00:00 2001 From: Alessandro Ranellucci Date: Wed, 15 Jan 2014 20:31:38 +0100 Subject: [PATCH] Refactored mesh slicing code into a new TriangleMeshSlicer class --- xs/src/TriangleMesh.cpp | 902 ++++++++++++++++++++-------------------- xs/src/TriangleMesh.hpp | 20 +- xs/t/01_trianglemesh.t | 2 +- xs/xsp/TriangleMesh.xsp | 3 +- 4 files changed, 479 insertions(+), 448 deletions(-) diff --git a/xs/src/TriangleMesh.cpp b/xs/src/TriangleMesh.cpp index 58a4a588af..40578baa15 100644 --- a/xs/src/TriangleMesh.cpp +++ b/xs/src/TriangleMesh.cpp @@ -164,450 +164,6 @@ void TriangleMesh::rotate(double angle, Point* center) this->translate(+center->x, +center->y, 0); } -void -TriangleMesh::slice(const std::vector &z, std::vector* layers) -{ - /* - This method gets called with a list of unscaled Z coordinates and outputs - a vector pointer having the same number of items as the original list. - Each item is a vector of polygons created by slicing our mesh at the - given heights. - - This method should basically combine the behavior of the existing - Perl methods defined in lib/Slic3r/TriangleMesh.pm: - - - analyze(): this creates the 'facets_edges' and the 'edges_facets' - tables (we don't need the 'edges' table) - - - slice_facet(): this has to be done for each facet. It generates - intersection lines with each plane identified by the Z list. - The get_layer_range() binary search used to identify the Z range - of the facet is already ported to C++ (see Object.xsp) - - - make_loops(): this has to be done for each layer. It creates polygons - from the lines generated by the previous step. - - At the end, we free the tables generated by analyze() as we don't - need them anymore. - FUTURE: parallelize slice_facet() and make_loops() - - NOTE: this method accepts a vector of floats because the mesh coordinate - type is float. - */ - - // build a table to map a facet_idx to its three edge indices - this->require_shared_vertices(); - typedef std::pair t_edge; - typedef std::vector t_edges; // edge_idx => a_id,b_id - typedef std::map t_edges_map; // a_id,b_id => edge_idx - typedef std::vector< std::vector > t_facets_edges; - t_facets_edges facets_edges; - - facets_edges.resize(this->stl.stats.number_of_facets); - - { - t_edges edges; - // reserve() instad of resize() because otherwise we couldn't read .size() below to assign edge_idx - edges.reserve(this->stl.stats.number_of_facets * 3); // number of edges = number of facets * 3 - t_edges_map edges_map; - for (int facet_idx = 0; facet_idx < this->stl.stats.number_of_facets; facet_idx++) { - facets_edges[facet_idx].resize(3); - for (int i = 0; i <= 2; i++) { - int a_id = this->stl.v_indices[facet_idx].vertex[i]; - int b_id = this->stl.v_indices[facet_idx].vertex[(i+1) % 3]; - - int edge_idx; - t_edges_map::const_iterator my_edge = edges_map.find(std::make_pair(b_id,a_id)); - if (my_edge != edges_map.end()) { - edge_idx = my_edge->second; - } else { - /* admesh can assign the same edge ID to more than two facets (which is - still topologically correct), so we have to search for a duplicate of - this edge too in case it was already seen in this orientation */ - my_edge = edges_map.find(std::make_pair(a_id,b_id)); - - if (my_edge != edges_map.end()) { - edge_idx = my_edge->second; - } else { - // edge isn't listed in table, so we insert it - edge_idx = edges.size(); - edges.push_back(std::make_pair(a_id,b_id)); - edges_map[ edges[edge_idx] ] = edge_idx; - } - } - facets_edges[facet_idx][i] = edge_idx; - - #ifdef SLIC3R_DEBUG - printf(" [facet %d, edge %d] a_id = %d, b_id = %d --> edge %d\n", facet_idx, i, a_id, b_id, edge_idx); - #endif - } - } - } - - std::vector lines(z.size()); - - // clone shared vertices coordinates and scale them - stl_vertex* v_scaled_shared = (stl_vertex*)calloc(this->stl.stats.shared_vertices, sizeof(stl_vertex)); - std::copy(this->stl.v_shared, this->stl.v_shared + this->stl.stats.shared_vertices, v_scaled_shared); - for (int i = 0; i < this->stl.stats.shared_vertices; i++) { - v_scaled_shared[i].x /= SCALING_FACTOR; - v_scaled_shared[i].y /= SCALING_FACTOR; - v_scaled_shared[i].z /= SCALING_FACTOR; - } - - for (int facet_idx = 0; facet_idx < this->stl.stats.number_of_facets; facet_idx++) { - stl_facet* facet = &this->stl.facet_start[facet_idx]; - - // find facet extents - float min_z = fminf(facet->vertex[0].z, fminf(facet->vertex[1].z, facet->vertex[2].z)); - float max_z = fmaxf(facet->vertex[0].z, fmaxf(facet->vertex[1].z, facet->vertex[2].z)); - - #ifdef SLIC3R_DEBUG - printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx, - facet->vertex[0].x, facet->vertex[0].y, facet->vertex[0].z, - facet->vertex[1].x, facet->vertex[1].y, facet->vertex[1].z, - facet->vertex[2].x, facet->vertex[2].y, facet->vertex[2].z); - printf("z: min = %.2f, max = %.2f\n", min_z, max_z); - #endif - - // find layer extents - std::vector::const_iterator min_layer, max_layer; - min_layer = std::lower_bound(z.begin(), z.end(), min_z); // first layer whose slice_z is >= min_z - max_layer = std::upper_bound(z.begin() + (min_layer - z.begin()), z.end(), max_z) - 1; // last layer whose slice_z is <= max_z - #ifdef SLIC3R_DEBUG - printf("layers: min = %d, max = %d\n", (int)(min_layer - z.begin()), (int)(max_layer - z.begin())); - #endif - - for (std::vector::const_iterator it = min_layer; it != max_layer + 1; ++it) { - std::vector::size_type layer_idx = it - z.begin(); - float slice_z = *it / SCALING_FACTOR; - std::vector points; - std::vector< std::vector::size_type > points_on_layer; - bool found_horizontal_edge = false; - - /* reorder vertices so that the first one is the one with lowest Z - this is needed to get all intersection lines in a consistent order - (external on the right of the line) */ - int i = 0; - if (facet->vertex[1].z == min_z) { - // vertex 1 has lowest Z - i = 1; - } else if (facet->vertex[2].z == min_z) { - // vertex 2 has lowest Z - i = 2; - } - for (int j = i; (j-i) < 3; j++) { // loop through facet edges - int edge_id = facets_edges[facet_idx][j % 3]; - int a_id = this->stl.v_indices[facet_idx].vertex[j % 3]; - int b_id = this->stl.v_indices[facet_idx].vertex[(j+1) % 3]; - stl_vertex* a = &v_scaled_shared[a_id]; - stl_vertex* b = &v_scaled_shared[b_id]; - - if (a->z == b->z && a->z == slice_z) { - // edge is horizontal and belongs to the current layer - - /* We assume that this method is never being called for horizontal - facets, so no other edge is going to be on this layer. */ - stl_vertex* v0 = &v_scaled_shared[ this->stl.v_indices[facet_idx].vertex[0] ]; - stl_vertex* v1 = &v_scaled_shared[ this->stl.v_indices[facet_idx].vertex[1] ]; - stl_vertex* v2 = &v_scaled_shared[ this->stl.v_indices[facet_idx].vertex[2] ]; - IntersectionLine line; - if (min_z == max_z) { - line.edge_type = feHorizontal; - } else if (v0->z < slice_z || v1->z < slice_z || v2->z < slice_z) { - line.edge_type = feTop; - std::swap(a, b); - std::swap(a_id, b_id); - } else { - line.edge_type = feBottom; - } - line.a.x = a->x; - line.a.y = a->y; - line.b.x = b->x; - line.b.y = b->y; - line.a_id = a_id; - line.b_id = b_id; - - lines[layer_idx].push_back(line); - found_horizontal_edge = true; - - // if this is a top or bottom edge, we can stop looping through edges - // because we won't find anything interesting - if (line.edge_type != feHorizontal) break; - } else if (a->z == slice_z) { - IntersectionPoint point; - point.x = a->x; - point.y = a->y; - point.point_id = a_id; - points.push_back(point); - points_on_layer.push_back(points.size()-1); - } else if (b->z == slice_z) { - IntersectionPoint point; - point.x = b->x; - point.y = b->y; - point.point_id = b_id; - points.push_back(point); - points_on_layer.push_back(points.size()-1); - } else if ((a->z < slice_z && b->z > slice_z) || (b->z < slice_z && a->z > slice_z)) { - // edge intersects the current layer; calculate intersection - - IntersectionPoint point; - point.x = b->x + (a->x - b->x) * (slice_z - b->z) / (a->z - b->z); - point.y = b->y + (a->y - b->y) * (slice_z - b->z) / (a->z - b->z); - point.edge_id = edge_id; - points.push_back(point); - } - } - if (found_horizontal_edge) continue; - - if (!points_on_layer.empty()) { - // we can't have only one point on layer because each vertex gets detected - // twice (once for each edge), and we can't have three points on layer because - // we assume this code is not getting called for horizontal facets - assert(points_on_layer.size() == 2); - assert( points[ points_on_layer[0] ].point_id == points[ points_on_layer[1] ].point_id ); - if (points.size() < 3) continue; // no intersection point, this is a V-shaped facet tangent to plane - points.erase( points.begin() + points_on_layer[1] ); - } - - if (!points.empty()) { - assert(points.size() == 2); // facets must intersect each plane 0 or 2 times - IntersectionLine line; - line.a.x = points[1].x; - line.a.y = points[1].y; - line.b.x = points[0].x; - line.b.y = points[0].y; - line.a_id = points[1].point_id; - line.b_id = points[0].point_id; - line.edge_a_id = points[1].edge_id; - line.edge_b_id = points[0].edge_id; - lines[layer_idx].push_back(line); - } - } - } - - free(v_scaled_shared); - - // build loops - layers->resize(z.size()); - for (std::vector::iterator it = lines.begin(); it != lines.end(); ++it) { - int layer_idx = it - lines.begin(); - #ifdef SLIC3R_DEBUG - printf("Layer %d:\n", layer_idx); - #endif - - /* - SVG svg("lines.svg"); - for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) - svg.AddLine(*line); - svg.Close(); - */ - - // remove tangent edges - for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { - if (line->skip || line->edge_type == feNone) continue; - - /* if the line is a facet edge, find another facet edge - having the same endpoints but in reverse order */ - for (IntersectionLines::iterator line2 = line + 1; line2 != it->end(); ++line2) { - if (line2->skip || line2->edge_type == feNone) continue; - - // are these facets adjacent? (sharing a common edge on this layer) - if (line->a_id == line2->a_id && line->b_id == line2->b_id) { - line2->skip = true; - - /* if they are both oriented upwards or downwards (like a 'V') - then we can remove both edges from this layer since it won't - affect the sliced shape */ - /* if one of them is oriented upwards and the other is oriented - downwards, let's only keep one of them (it doesn't matter which - one since all 'top' lines were reversed at slicing) */ - if (line->edge_type == line2->edge_type) { - line->skip = true; - break; - } - } else if (line->a_id == line2->b_id && line->b_id == line2->a_id) { - /* if this edge joins two horizontal facets, remove both of them */ - if (line->edge_type == feHorizontal && line2->edge_type == feHorizontal) { - line->skip = true; - line2->skip = true; - break; - } - } - } - } - - // build a map of lines by edge_a_id and a_id - std::vector by_edge_a_id, by_a_id; - by_edge_a_id.resize(this->stl.stats.number_of_facets * 3); - by_a_id.resize(this->stl.stats.shared_vertices); - for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { - if (line->skip) continue; - if (line->edge_a_id != -1) by_edge_a_id[line->edge_a_id].push_back(&(*line)); - if (line->a_id != -1) by_a_id[line->a_id].push_back(&(*line)); - } - - CYCLE: while (1) { - // take first spare line and start a new loop - IntersectionLine* first_line = NULL; - for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { - if (line->skip) continue; - first_line = &(*line); - break; - } - if (first_line == NULL) break; - first_line->skip = true; - IntersectionLinePtrs loop; - loop.push_back(first_line); - - /* - printf("first_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", - first_line->edge_a_id, first_line->edge_b_id, first_line->a_id, first_line->b_id, - first_line->a.x, first_line->a.y, first_line->b.x, first_line->b.y); - */ - - while (1) { - // find a line starting where last one finishes - IntersectionLine* next_line = NULL; - if (loop.back()->edge_b_id != -1) { - IntersectionLinePtrs* candidates = &(by_edge_a_id[loop.back()->edge_b_id]); - for (IntersectionLinePtrs::iterator lineptr = candidates->begin(); lineptr != candidates->end(); ++lineptr) { - if ((*lineptr)->skip) continue; - next_line = *lineptr; - break; - } - } - if (next_line == NULL && loop.back()->b_id != -1) { - IntersectionLinePtrs* candidates = &(by_a_id[loop.back()->b_id]); - for (IntersectionLinePtrs::iterator lineptr = candidates->begin(); lineptr != candidates->end(); ++lineptr) { - if ((*lineptr)->skip) continue; - next_line = *lineptr; - break; - } - } - - if (next_line == NULL) { - // check whether we closed this loop - if ((loop.front()->edge_a_id != -1 && loop.front()->edge_a_id == loop.back()->edge_b_id) - || (loop.front()->a_id != -1 && loop.front()->a_id == loop.back()->b_id)) { - // loop is complete - Polygon p; - p.points.reserve(loop.size()); - for (IntersectionLinePtrs::iterator lineptr = loop.begin(); lineptr != loop.end(); ++lineptr) { - p.points.push_back((*lineptr)->a); - } - (*layers)[layer_idx].push_back(p); - - #ifdef SLIC3R_DEBUG - printf(" Discovered %s polygon of %d points\n", (p.is_counter_clockwise() ? "ccw" : "cw"), (int)p.points.size()); - #endif - - goto CYCLE; - } - - // we can't close this loop! - //// push @failed_loops, [@loop]; - #ifdef SLIC3R_DEBUG - printf(" Unable to close this loop having %d points\n", (int)loop.size()); - #endif - goto CYCLE; - } - /* - printf("next_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", - next_line->edge_a_id, next_line->edge_b_id, next_line->a_id, next_line->b_id, - next_line->a.x, next_line->a.y, next_line->b.x, next_line->b.y); - */ - loop.push_back(next_line); - next_line->skip = true; - } - } - } -} - -class _area_comp { - public: - _area_comp(std::vector* _aa) : abs_area(_aa) {}; - bool operator() (const size_t &a, const size_t &b) { - return (*this->abs_area)[a] > (*this->abs_area)[b]; - } - - private: - std::vector* abs_area; -}; - -void -TriangleMesh::slice(const std::vector &z, std::vector* layers) -{ - std::vector layers_p; - this->slice(z, &layers_p); - - /* - Input loops are not suitable for evenodd nor nonzero fill types, as we might get - two consecutive concentric loops having the same winding order - and we have to - respect such order. In that case, evenodd would create wrong inversions, and nonzero - would ignore holes inside two concentric contours. - So we're ordering loops and collapse consecutive concentric loops having the same - winding order. - TODO: find a faster algorithm for this, maybe with some sort of binary search. - If we computed a "nesting tree" we could also just remove the consecutive loops - having the same winding order, and remove the extra one(s) so that we could just - supply everything to offset_ex() instead of performing several union/diff calls. - - we sort by area assuming that the outermost loops have larger area; - the previous sorting method, based on $b->contains_point($a->[0]), failed to nest - loops correctly in some edge cases when original model had overlapping facets - */ - - layers->resize(z.size()); - - for (std::vector::const_iterator loops = layers_p.begin(); loops != layers_p.end(); ++loops) { - size_t layer_id = loops - layers_p.begin(); - - std::vector area; - std::vector abs_area; - std::vector sorted_area; // vector of indices - for (Polygons::const_iterator loop = loops->begin(); loop != loops->end(); ++loop) { - double a = loop->area(); - area.push_back(a); - abs_area.push_back(std::fabs(a)); - sorted_area.push_back(loop - loops->begin()); - } - - std::sort(sorted_area.begin(), sorted_area.end(), _area_comp(&abs_area)); // outer first - - // we don't perform a safety offset now because it might reverse cw loops - Polygons slices; - for (std::vector::const_iterator loop_idx = sorted_area.begin(); loop_idx != sorted_area.end(); ++loop_idx) { - /* we rely on the already computed area to determine the winding order - of the loops, since the Orientation() function provided by Clipper - would do the same, thus repeating the calculation */ - Polygons::const_iterator loop = loops->begin() + *loop_idx; - if (area[*loop_idx] >= 0) { - slices.push_back(*loop); - } else { - diff(slices, *loop, slices); - } - } - - // perform a safety offset to merge very close facets (TODO: find test case for this) - double safety_offset = scale_(0.0499); - ExPolygons ex_slices; - offset2_ex(slices, ex_slices, +safety_offset, -safety_offset); - - #ifdef SLIC3R_DEBUG - size_t holes_count = 0; - for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++e) { - holes_count += e->holes.size(); - } - printf("Layer %zu (slice_z = %.2f): %zu surface(s) having %zu holes detected from %zu polylines\n", - layer_id, z[layer_id], ex_slices.size(), holes_count, loops->size()); - #endif - - ExPolygons* layer = &(*layers)[layer_id]; - layer->insert(layer->end(), ex_slices.begin(), ex_slices.end()); - } -} - TriangleMeshPtrs TriangleMesh::split() const { @@ -777,4 +333,462 @@ void TriangleMesh::ReadFromPerl(SV* vertices, SV* facets) } #endif +void +TriangleMeshSlicer::slice(const std::vector &z, std::vector* layers) +{ + /* + This method gets called with a list of unscaled Z coordinates and outputs + a vector pointer having the same number of items as the original list. + Each item is a vector of polygons created by slicing our mesh at the + given heights. + + This method should basically combine the behavior of the existing + Perl methods defined in lib/Slic3r/TriangleMesh.pm: + + - analyze(): this creates the 'facets_edges' and the 'edges_facets' + tables (we don't need the 'edges' table) + + - slice_facet(): this has to be done for each facet. It generates + intersection lines with each plane identified by the Z list. + The get_layer_range() binary search used to identify the Z range + of the facet is already ported to C++ (see Object.xsp) + + - make_loops(): this has to be done for each layer. It creates polygons + from the lines generated by the previous step. + + At the end, we free the tables generated by analyze() as we don't + need them anymore. + FUTURE: parallelize slice_facet() and make_loops() + + NOTE: this method accepts a vector of floats because the mesh coordinate + type is float. + */ + + std::vector lines(z.size()); + + for (int facet_idx = 0; facet_idx < this->mesh->stl.stats.number_of_facets; facet_idx++) { + stl_facet* facet = &this->mesh->stl.facet_start[facet_idx]; + + // find facet extents + float min_z = fminf(facet->vertex[0].z, fminf(facet->vertex[1].z, facet->vertex[2].z)); + float max_z = fmaxf(facet->vertex[0].z, fmaxf(facet->vertex[1].z, facet->vertex[2].z)); + + #ifdef SLIC3R_DEBUG + printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx, + facet->vertex[0].x, facet->vertex[0].y, facet->vertex[0].z, + facet->vertex[1].x, facet->vertex[1].y, facet->vertex[1].z, + facet->vertex[2].x, facet->vertex[2].y, facet->vertex[2].z); + printf("z: min = %.2f, max = %.2f\n", min_z, max_z); + #endif + + // find layer extents + std::vector::const_iterator min_layer, max_layer; + min_layer = std::lower_bound(z.begin(), z.end(), min_z); // first layer whose slice_z is >= min_z + max_layer = std::upper_bound(z.begin() + (min_layer - z.begin()), z.end(), max_z) - 1; // last layer whose slice_z is <= max_z + #ifdef SLIC3R_DEBUG + printf("layers: min = %d, max = %d\n", (int)(min_layer - z.begin()), (int)(max_layer - z.begin())); + #endif + + for (std::vector::const_iterator it = min_layer; it != max_layer + 1; ++it) { + std::vector::size_type layer_idx = it - z.begin(); + this->slice_facet(*it / SCALING_FACTOR, *facet, facet_idx, min_z, max_z, &lines[layer_idx]); + } + } + + // v_scaled_shared could be freed here + + // build loops + layers->resize(z.size()); + for (std::vector::iterator it = lines.begin(); it != lines.end(); ++it) { + int layer_idx = it - lines.begin(); + #ifdef SLIC3R_DEBUG + printf("Layer %d:\n", layer_idx); + #endif + + /* + SVG svg("lines.svg"); + for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { + svg.AddLine(*line); + } + svg.Close(); + */ + + // remove tangent edges + for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { + if (line->skip || line->edge_type == feNone) continue; + + /* if the line is a facet edge, find another facet edge + having the same endpoints but in reverse order */ + for (IntersectionLines::iterator line2 = line + 1; line2 != it->end(); ++line2) { + if (line2->skip || line2->edge_type == feNone) continue; + + // are these facets adjacent? (sharing a common edge on this layer) + if (line->a_id == line2->a_id && line->b_id == line2->b_id) { + line2->skip = true; + + /* if they are both oriented upwards or downwards (like a 'V') + then we can remove both edges from this layer since it won't + affect the sliced shape */ + /* if one of them is oriented upwards and the other is oriented + downwards, let's only keep one of them (it doesn't matter which + one since all 'top' lines were reversed at slicing) */ + if (line->edge_type == line2->edge_type) { + line->skip = true; + break; + } + } else if (line->a_id == line2->b_id && line->b_id == line2->a_id) { + /* if this edge joins two horizontal facets, remove both of them */ + if (line->edge_type == feHorizontal && line2->edge_type == feHorizontal) { + line->skip = true; + line2->skip = true; + break; + } + } + } + } + + // build a map of lines by edge_a_id and a_id + std::vector by_edge_a_id, by_a_id; + by_edge_a_id.resize(this->mesh->stl.stats.number_of_facets * 3); + by_a_id.resize(this->mesh->stl.stats.shared_vertices); + for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { + if (line->skip) continue; + if (line->edge_a_id != -1) by_edge_a_id[line->edge_a_id].push_back(&(*line)); + if (line->a_id != -1) by_a_id[line->a_id].push_back(&(*line)); + } + + CYCLE: while (1) { + // take first spare line and start a new loop + IntersectionLine* first_line = NULL; + for (IntersectionLines::iterator line = it->begin(); line != it->end(); ++line) { + if (line->skip) continue; + first_line = &(*line); + break; + } + if (first_line == NULL) break; + first_line->skip = true; + IntersectionLinePtrs loop; + loop.push_back(first_line); + + /* + printf("first_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", + first_line->edge_a_id, first_line->edge_b_id, first_line->a_id, first_line->b_id, + first_line->a.x, first_line->a.y, first_line->b.x, first_line->b.y); + */ + + while (1) { + // find a line starting where last one finishes + IntersectionLine* next_line = NULL; + if (loop.back()->edge_b_id != -1) { + IntersectionLinePtrs* candidates = &(by_edge_a_id[loop.back()->edge_b_id]); + for (IntersectionLinePtrs::iterator lineptr = candidates->begin(); lineptr != candidates->end(); ++lineptr) { + if ((*lineptr)->skip) continue; + next_line = *lineptr; + break; + } + } + if (next_line == NULL && loop.back()->b_id != -1) { + IntersectionLinePtrs* candidates = &(by_a_id[loop.back()->b_id]); + for (IntersectionLinePtrs::iterator lineptr = candidates->begin(); lineptr != candidates->end(); ++lineptr) { + if ((*lineptr)->skip) continue; + next_line = *lineptr; + break; + } + } + + if (next_line == NULL) { + // check whether we closed this loop + if ((loop.front()->edge_a_id != -1 && loop.front()->edge_a_id == loop.back()->edge_b_id) + || (loop.front()->a_id != -1 && loop.front()->a_id == loop.back()->b_id)) { + // loop is complete + Polygon p; + p.points.reserve(loop.size()); + for (IntersectionLinePtrs::iterator lineptr = loop.begin(); lineptr != loop.end(); ++lineptr) { + p.points.push_back((*lineptr)->a); + } + (*layers)[layer_idx].push_back(p); + + #ifdef SLIC3R_DEBUG + printf(" Discovered %s polygon of %d points\n", (p.is_counter_clockwise() ? "ccw" : "cw"), (int)p.points.size()); + #endif + + goto CYCLE; + } + + // we can't close this loop! + //// push @failed_loops, [@loop]; + //#ifdef SLIC3R_DEBUG + printf(" Unable to close this loop having %d points\n", (int)loop.size()); + //#endif + goto CYCLE; + } + /* + printf("next_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", + next_line->edge_a_id, next_line->edge_b_id, next_line->a_id, next_line->b_id, + next_line->a.x, next_line->a.y, next_line->b.x, next_line->b.y); + */ + loop.push_back(next_line); + next_line->skip = true; + } + } + } +} + +void +TriangleMeshSlicer::slice_facet(float slice_z, const stl_facet &facet, const int &facet_idx, const float &min_z, const float &max_z, std::vector* lines) const +{ + std::vector points; + std::vector< std::vector::size_type > points_on_layer; + bool found_horizontal_edge = false; + + /* reorder vertices so that the first one is the one with lowest Z + this is needed to get all intersection lines in a consistent order + (external on the right of the line) */ + int i = 0; + if (facet.vertex[1].z == min_z) { + // vertex 1 has lowest Z + i = 1; + } else if (facet.vertex[2].z == min_z) { + // vertex 2 has lowest Z + i = 2; + } + for (int j = i; (j-i) < 3; j++) { // loop through facet edges + int edge_id = this->facets_edges[facet_idx][j % 3]; + int a_id = this->mesh->stl.v_indices[facet_idx].vertex[j % 3]; + int b_id = this->mesh->stl.v_indices[facet_idx].vertex[(j+1) % 3]; + stl_vertex* a = &this->v_scaled_shared[a_id]; + stl_vertex* b = &this->v_scaled_shared[b_id]; + + if (a->z == b->z && a->z == slice_z) { + // edge is horizontal and belongs to the current layer + + /* We assume that this method is never being called for horizontal + facets, so no other edge is going to be on this layer. */ + stl_vertex* v0 = &this->v_scaled_shared[ this->mesh->stl.v_indices[facet_idx].vertex[0] ]; + stl_vertex* v1 = &this->v_scaled_shared[ this->mesh->stl.v_indices[facet_idx].vertex[1] ]; + stl_vertex* v2 = &this->v_scaled_shared[ this->mesh->stl.v_indices[facet_idx].vertex[2] ]; + IntersectionLine line; + if (min_z == max_z) { + line.edge_type = feHorizontal; + } else if (v0->z < slice_z || v1->z < slice_z || v2->z < slice_z) { + line.edge_type = feTop; + std::swap(a, b); + std::swap(a_id, b_id); + } else { + line.edge_type = feBottom; + } + line.a.x = a->x; + line.a.y = a->y; + line.b.x = b->x; + line.b.y = b->y; + line.a_id = a_id; + line.b_id = b_id; + lines->push_back(line); + + found_horizontal_edge = true; + + // if this is a top or bottom edge, we can stop looping through edges + // because we won't find anything interesting + + if (line.edge_type != feHorizontal) return; + } else if (a->z == slice_z) { + IntersectionPoint point; + point.x = a->x; + point.y = a->y; + point.point_id = a_id; + points.push_back(point); + points_on_layer.push_back(points.size()-1); + } else if (b->z == slice_z) { + IntersectionPoint point; + point.x = b->x; + point.y = b->y; + point.point_id = b_id; + points.push_back(point); + points_on_layer.push_back(points.size()-1); + } else if ((a->z < slice_z && b->z > slice_z) || (b->z < slice_z && a->z > slice_z)) { + // edge intersects the current layer; calculate intersection + + IntersectionPoint point; + point.x = b->x + (a->x - b->x) * (slice_z - b->z) / (a->z - b->z); + point.y = b->y + (a->y - b->y) * (slice_z - b->z) / (a->z - b->z); + point.edge_id = edge_id; + points.push_back(point); + } + } + if (found_horizontal_edge) return; + + if (!points_on_layer.empty()) { + // we can't have only one point on layer because each vertex gets detected + // twice (once for each edge), and we can't have three points on layer because + // we assume this code is not getting called for horizontal facets + assert(points_on_layer.size() == 2); + assert( points[ points_on_layer[0] ].point_id == points[ points_on_layer[1] ].point_id ); + if (points.size() < 3) return; // no intersection point, this is a V-shaped facet tangent to plane + points.erase( points.begin() + points_on_layer[1] ); + } + + if (!points.empty()) { + assert(points.size() == 2); // facets must intersect each plane 0 or 2 times + IntersectionLine line; + line.a.x = points[1].x; + line.a.y = points[1].y; + line.b.x = points[0].x; + line.b.y = points[0].y; + line.a_id = points[1].point_id; + line.b_id = points[0].point_id; + line.edge_a_id = points[1].edge_id; + line.edge_b_id = points[0].edge_id; + lines->push_back(line); + return; + } +} + +class _area_comp { + public: + _area_comp(std::vector* _aa) : abs_area(_aa) {}; + bool operator() (const size_t &a, const size_t &b) { + return (*this->abs_area)[a] > (*this->abs_area)[b]; + } + + private: + std::vector* abs_area; +}; + +void +TriangleMeshSlicer::slice(const std::vector &z, std::vector* layers) +{ + std::vector layers_p; + this->slice(z, &layers_p); + + /* + Input loops are not suitable for evenodd nor nonzero fill types, as we might get + two consecutive concentric loops having the same winding order - and we have to + respect such order. In that case, evenodd would create wrong inversions, and nonzero + would ignore holes inside two concentric contours. + So we're ordering loops and collapse consecutive concentric loops having the same + winding order. + TODO: find a faster algorithm for this, maybe with some sort of binary search. + If we computed a "nesting tree" we could also just remove the consecutive loops + having the same winding order, and remove the extra one(s) so that we could just + supply everything to offset_ex() instead of performing several union/diff calls. + + we sort by area assuming that the outermost loops have larger area; + the previous sorting method, based on $b->contains_point($a->[0]), failed to nest + loops correctly in some edge cases when original model had overlapping facets + */ + + layers->resize(z.size()); + + for (std::vector::const_iterator loops = layers_p.begin(); loops != layers_p.end(); ++loops) { + size_t layer_id = loops - layers_p.begin(); + + std::vector area; + std::vector abs_area; + std::vector sorted_area; // vector of indices + for (Polygons::const_iterator loop = loops->begin(); loop != loops->end(); ++loop) { + double a = loop->area(); + area.push_back(a); + abs_area.push_back(std::fabs(a)); + sorted_area.push_back(loop - loops->begin()); + } + + std::sort(sorted_area.begin(), sorted_area.end(), _area_comp(&abs_area)); // outer first + + // we don't perform a safety offset now because it might reverse cw loops + Polygons slices; + for (std::vector::const_iterator loop_idx = sorted_area.begin(); loop_idx != sorted_area.end(); ++loop_idx) { + /* we rely on the already computed area to determine the winding order + of the loops, since the Orientation() function provided by Clipper + would do the same, thus repeating the calculation */ + Polygons::const_iterator loop = loops->begin() + *loop_idx; + if (area[*loop_idx] >= 0) { + slices.push_back(*loop); + } else { + diff(slices, *loop, slices); + } + } + + // perform a safety offset to merge very close facets (TODO: find test case for this) + double safety_offset = scale_(0.0499); + ExPolygons ex_slices; + offset2_ex(slices, ex_slices, +safety_offset, -safety_offset); + + #ifdef SLIC3R_DEBUG + size_t holes_count = 0; + for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++e) { + holes_count += e->holes.size(); + } + printf("Layer %zu (slice_z = %.2f): %zu surface(s) having %zu holes detected from %zu polylines\n", + layer_id, z[layer_id], ex_slices.size(), holes_count, loops->size()); + #endif + + ExPolygons* layer = &(*layers)[layer_id]; + layer->insert(layer->end(), ex_slices.begin(), ex_slices.end()); + } +} + +TriangleMeshSlicer::TriangleMeshSlicer(TriangleMesh* _mesh) : mesh(_mesh), v_scaled_shared(NULL) +{ + // build a table to map a facet_idx to its three edge indices + this->mesh->require_shared_vertices(); + typedef std::pair t_edge; + typedef std::vector t_edges; // edge_idx => a_id,b_id + typedef std::map t_edges_map; // a_id,b_id => edge_idx + + this->facets_edges.resize(this->mesh->stl.stats.number_of_facets); + + { + t_edges edges; + // reserve() instad of resize() because otherwise we couldn't read .size() below to assign edge_idx + edges.reserve(this->mesh->stl.stats.number_of_facets * 3); // number of edges = number of facets * 3 + t_edges_map edges_map; + for (int facet_idx = 0; facet_idx < this->mesh->stl.stats.number_of_facets; facet_idx++) { + this->facets_edges[facet_idx].resize(3); + for (int i = 0; i <= 2; i++) { + int a_id = this->mesh->stl.v_indices[facet_idx].vertex[i]; + int b_id = this->mesh->stl.v_indices[facet_idx].vertex[(i+1) % 3]; + + int edge_idx; + t_edges_map::const_iterator my_edge = edges_map.find(std::make_pair(b_id,a_id)); + if (my_edge != edges_map.end()) { + edge_idx = my_edge->second; + } else { + /* admesh can assign the same edge ID to more than two facets (which is + still topologically correct), so we have to search for a duplicate of + this edge too in case it was already seen in this orientation */ + my_edge = edges_map.find(std::make_pair(a_id,b_id)); + + if (my_edge != edges_map.end()) { + edge_idx = my_edge->second; + } else { + // edge isn't listed in table, so we insert it + edge_idx = edges.size(); + edges.push_back(std::make_pair(a_id,b_id)); + edges_map[ edges[edge_idx] ] = edge_idx; + } + } + this->facets_edges[facet_idx][i] = edge_idx; + + #ifdef SLIC3R_DEBUG + printf(" [facet %d, edge %d] a_id = %d, b_id = %d --> edge %d\n", facet_idx, i, a_id, b_id, edge_idx); + #endif + } + } + } + + // clone shared vertices coordinates and scale them + this->v_scaled_shared = (stl_vertex*)calloc(this->mesh->stl.stats.shared_vertices, sizeof(stl_vertex)); + std::copy(this->mesh->stl.v_shared, this->mesh->stl.v_shared + this->mesh->stl.stats.shared_vertices, this->v_scaled_shared); + for (int i = 0; i < this->mesh->stl.stats.shared_vertices; i++) { + this->v_scaled_shared[i].x /= SCALING_FACTOR; + this->v_scaled_shared[i].y /= SCALING_FACTOR; + this->v_scaled_shared[i].z /= SCALING_FACTOR; + } +} + +TriangleMeshSlicer::~TriangleMeshSlicer() +{ + if (this->v_scaled_shared != NULL) free(this->v_scaled_shared); +} + } diff --git a/xs/src/TriangleMesh.hpp b/xs/src/TriangleMesh.hpp index 7efb0945e2..b98357d317 100644 --- a/xs/src/TriangleMesh.hpp +++ b/xs/src/TriangleMesh.hpp @@ -12,6 +12,7 @@ namespace Slic3r { class TriangleMesh; +class TriangleMeshSlicer; typedef std::vector TriangleMeshPtrs; class TriangleMesh @@ -30,8 +31,6 @@ class TriangleMesh void translate(float x, float y, float z); void align_to_origin(); void rotate(double angle, Point* center); - void slice(const std::vector &z, std::vector* layers); - void slice(const std::vector &z, std::vector* layers); TriangleMeshPtrs split() const; void merge(const TriangleMesh* mesh); void horizontal_projection(ExPolygons &retval) const; @@ -47,6 +46,7 @@ class TriangleMesh private: void require_shared_vertices(); + friend class TriangleMeshSlicer; }; enum FacetEdgeType { feNone, feTop, feBottom, feHorizontal }; @@ -75,6 +75,22 @@ class IntersectionLine typedef std::vector IntersectionLines; typedef std::vector IntersectionLinePtrs; +class TriangleMeshSlicer +{ + public: + TriangleMesh* mesh; + TriangleMeshSlicer(TriangleMesh* _mesh); + ~TriangleMeshSlicer(); + void slice(const std::vector &z, std::vector* layers); + void slice(const std::vector &z, std::vector* layers); + void slice_facet(float slice_z, const stl_facet &facet, const int &facet_idx, const float &min_z, const float &max_z, std::vector* lines) const; + + private: + typedef std::vector< std::vector > t_facets_edges; + t_facets_edges facets_edges; + stl_vertex* v_scaled_shared; +}; + } #endif diff --git a/xs/t/01_trianglemesh.t b/xs/t/01_trianglemesh.t index 86dcdaa4a5..500729a588 100644 --- a/xs/t/01_trianglemesh.t +++ b/xs/t/01_trianglemesh.t @@ -83,7 +83,7 @@ my $cube = { my $result = $m->slice(\@z); my $SCALING_FACTOR = 0.000001; for my $i (0..$#z) { - is scalar(@{$result->[$i]}), 1, 'number of returned polygons per layer'; + is scalar(@{$result->[$i]}), 1, "number of returned polygons per layer (z = " . $z[$i] . ")"; is $result->[$i][0]->area, 20*20/($SCALING_FACTOR**2), 'size of returned polygon'; } } diff --git a/xs/xsp/TriangleMesh.xsp b/xs/xsp/TriangleMesh.xsp index 26670a1db7..6b587d5fc4 100644 --- a/xs/xsp/TriangleMesh.xsp +++ b/xs/xsp/TriangleMesh.xsp @@ -142,7 +142,8 @@ TriangleMesh::slice(z) delete z; std::vector layers; - THIS->slice(z_f, &layers); + TriangleMeshSlicer mslicer(THIS); + mslicer.slice(z_f, &layers); AV* layers_av = newAV(); av_extend(layers_av, layers.size()-1);