From 48bc166d6d1ee257fc9f0d540ab2c349091fc900 Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Fri, 2 Nov 2018 11:57:57 +0100 Subject: [PATCH] Importing the SLA computing module into the native source tree. --- CMakeLists.txt | 12 + sandboxes/CMakeLists.txt | 2 + sandboxes/slabasebed/CMakeLists.txt | 2 + {src => sandboxes/slabasebed}/slabasebed.cpp | 0 sandboxes/slasupporttree/CMakeLists.txt | 2 + sandboxes/slasupporttree/slasupporttree.cpp | 48 + src/igl/AABB.cpp | 3 +- src/libnest2d/tools/benchmark.h | 58 + src/libslic3r/CMakeLists.txt | 10 + src/libslic3r/SLA/SLABasePool.cpp | 534 ++++++ src/libslic3r/SLA/SLABasePool.hpp | 37 + src/libslic3r/SLA/SLABoilerPlate.hpp | 92 + src/libslic3r/SLA/SLABoostAdapter.hpp | 132 ++ src/libslic3r/SLA/SLARotfinder.cpp | 156 ++ src/libslic3r/SLA/SLARotfinder.hpp | 38 + src/libslic3r/SLA/SLASpatIndex.hpp | 112 ++ src/libslic3r/SLA/SLASupportTree.cpp | 1665 ++++++++++++++++++ src/libslic3r/SLA/SLASupportTree.hpp | 132 ++ src/libslic3r/SLA/SLASupportTreeIGL.cpp | 237 +++ tests/CMakeLists.txt | 3 + 20 files changed, 3274 insertions(+), 1 deletion(-) create mode 100644 sandboxes/CMakeLists.txt create mode 100644 sandboxes/slabasebed/CMakeLists.txt rename {src => sandboxes/slabasebed}/slabasebed.cpp (100%) create mode 100644 sandboxes/slasupporttree/CMakeLists.txt create mode 100644 sandboxes/slasupporttree/slasupporttree.cpp create mode 100644 src/libnest2d/tools/benchmark.h create mode 100644 src/libslic3r/SLA/SLABasePool.cpp create mode 100644 src/libslic3r/SLA/SLABasePool.hpp create mode 100644 src/libslic3r/SLA/SLABoilerPlate.hpp create mode 100644 src/libslic3r/SLA/SLABoostAdapter.hpp create mode 100644 src/libslic3r/SLA/SLARotfinder.cpp create mode 100644 src/libslic3r/SLA/SLARotfinder.hpp create mode 100644 src/libslic3r/SLA/SLASpatIndex.hpp create mode 100644 src/libslic3r/SLA/SLASupportTree.cpp create mode 100644 src/libslic3r/SLA/SLASupportTree.hpp create mode 100644 src/libslic3r/SLA/SLASupportTreeIGL.cpp create mode 100644 tests/CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 0656b824e..2e1b6331e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,10 @@ option(SLIC3R_MSVC_COMPILE_PARALLEL "Compile on Visual Studio in parallel" 1) option(SLIC3R_MSVC_PDB "Generate PDB files on MSVC in Release mode" 1) option(SLIC3R_PERL_XS "Compile XS Perl module and enable Perl unit and integration tests" 0) +# Proposal for C++ unit tests and sandboxes +option(SLIC3R_BUILD_SANDBOXES "Build development sandboxes" OFF) +option(SLIC3R_BUILD_TESTS "Build unit tests" OFF) + if (MSVC) if (SLIC3R_MSVC_COMPILE_PARALLEL) add_compile_options(/MP) @@ -244,5 +248,13 @@ if (SLIC3R_PERL_XS) add_subdirectory(xs) endif () +if(SLIC3R_BUILD_SANDBOXES) + add_subdirectory(sandboxes) +endif() + +if(SLIC3R_BUILD_TESTS) + add_subdirectory(tests) +endif() + file(GLOB MyVar var/*.png) install(FILES ${MyVar} DESTINATION share/slic3r-prusa3d) diff --git a/sandboxes/CMakeLists.txt b/sandboxes/CMakeLists.txt new file mode 100644 index 000000000..5905c438e --- /dev/null +++ b/sandboxes/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(slabasebed) +add_subdirectory(slasupporttree) diff --git a/sandboxes/slabasebed/CMakeLists.txt b/sandboxes/slabasebed/CMakeLists.txt new file mode 100644 index 000000000..bff5ca588 --- /dev/null +++ b/sandboxes/slabasebed/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(slabasebed EXCLUDE_FROM_ALL slabasebed.cpp) +target_link_libraries(slabasebed libslic3r) \ No newline at end of file diff --git a/src/slabasebed.cpp b/sandboxes/slabasebed/slabasebed.cpp similarity index 100% rename from src/slabasebed.cpp rename to sandboxes/slabasebed/slabasebed.cpp diff --git a/sandboxes/slasupporttree/CMakeLists.txt b/sandboxes/slasupporttree/CMakeLists.txt new file mode 100644 index 000000000..a3323a83d --- /dev/null +++ b/sandboxes/slasupporttree/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(slasupporttree EXCLUDE_FROM_ALL slasupporttree.cpp) +target_link_libraries(slasupporttree libslic3r) \ No newline at end of file diff --git a/sandboxes/slasupporttree/slasupporttree.cpp b/sandboxes/slasupporttree/slasupporttree.cpp new file mode 100644 index 000000000..0cfc60f7a --- /dev/null +++ b/sandboxes/slasupporttree/slasupporttree.cpp @@ -0,0 +1,48 @@ +#include +#include + +#include +#include "TriangleMesh.hpp" +#include "Model.hpp" +#include "callback.hpp" +#include "SLA/SLASupportTree.hpp" +#include "benchmark.h" + +const std::string USAGE_STR = { + "Usage: slasupporttree stlfilename.stl" +}; + +void confess_at(const char * /*file*/, + int /*line*/, + const char * /*func*/, + const char * /*pat*/, + ...) {} + +namespace Slic3r { + +void PerlCallback::deregister_callback() {} +} + +int main(const int argc, const char *argv[]) { + using namespace Slic3r; + using std::cout; using std::endl; + + if(argc < 2) { + cout << USAGE_STR << endl; + return EXIT_SUCCESS; + } + + Benchmark bench; + TriangleMesh result; + + bench.start(); + sla::create_head(result, 3, 1, 4); + bench.stop(); + + cout << "Support tree creation time: " << std::setprecision(10) + << bench.getElapsedSec() << " seconds." << endl; + + result.write_ascii("out.stl"); + + return EXIT_SUCCESS; +} diff --git a/src/igl/AABB.cpp b/src/igl/AABB.cpp index 095373354..1b88e9faa 100644 --- a/src/igl/AABB.cpp +++ b/src/igl/AABB.cpp @@ -885,7 +885,8 @@ igl::AABB::intersect_ray( const RowVectorDIMS & dir, igl::Hit & hit) const { -#if false + // FIXME: Needs a proper path +#if /*false*/ 0 // BFS std::queue Q; // Or DFS diff --git a/src/libnest2d/tools/benchmark.h b/src/libnest2d/tools/benchmark.h new file mode 100644 index 000000000..19870b37b --- /dev/null +++ b/src/libnest2d/tools/benchmark.h @@ -0,0 +1,58 @@ +/* + * Copyright (C) Tamás Mészáros + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ +#ifndef INCLUDE_BENCHMARK_H_ +#define INCLUDE_BENCHMARK_H_ + +#include +#include + +/** + * A class for doing benchmarks. + */ +class Benchmark { + typedef std::chrono::high_resolution_clock Clock; + typedef Clock::duration Duration; + typedef Clock::time_point TimePoint; + + TimePoint t1, t2; + Duration d; + + inline double to_sec(Duration d) { + return d.count() * double(Duration::period::num) / Duration::period::den; + } + +public: + + /** + * Measure time from the moment of this call. + */ + void start() { t1 = Clock::now(); } + + /** + * Measure time to the moment of this call. + */ + void stop() { t2 = Clock::now(); } + + /** + * Get the time elapsed between a start() end a stop() call. + * @return Returns the elapsed time in seconds. + */ + double getElapsedSec() { d = t2 - t1; return to_sec(d); } +}; + + +#endif /* INCLUDE_BENCHMARK_H_ */ diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index c9953e38a..ab6098708 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -153,6 +153,16 @@ add_library(libslic3r STATIC SLABasePool.cpp utils.cpp Utils.hpp + SLA/SLABoilerPlate.hpp + SLA/SLABasePool.hpp + SLA/SLABasePool.cpp + SLA/SLASupportTree.hpp + SLA/SLASupportTree.cpp + SLA/SLASupportTreeIGL.cpp + SLA/SLARotfinder.hpp + SLA/SLARotfinder.cpp + SLA/SLABoostAdapter.hpp + SLA/SLASpatIndex.hpp ) add_precompiled_header(libslic3r pchheader.hpp FORCEINCLUDE) diff --git a/src/libslic3r/SLA/SLABasePool.cpp b/src/libslic3r/SLA/SLABasePool.cpp new file mode 100644 index 000000000..fe31ed57d --- /dev/null +++ b/src/libslic3r/SLA/SLABasePool.cpp @@ -0,0 +1,534 @@ +#include "SLABasePool.hpp" +#include "SLABoilerPlate.hpp" + +#include "boost/log/trivial.hpp" +#include "SLABoostAdapter.hpp" +#include "ClipperUtils.hpp" + +//#include "SVG.hpp" +//#include "benchmark.h" + +namespace Slic3r { namespace sla { + +/// Convert the triangulation output to an intermediate mesh. +Contour3D convert(const Polygons& triangles, coord_t z, bool dir) { + + Pointf3s points; + points.reserve(3*triangles.size()); + Indices indices; + indices.reserve(points.size()); + + for(auto& tr : triangles) { + auto c = coord_t(points.size()), b = c++, a = c++; + if(dir) indices.emplace_back(a, b, c); + else indices.emplace_back(c, b, a); + for(auto& p : tr.points) { + points.emplace_back(unscale(x(p), y(p), z)); + } + } + + return {points, indices}; +} + +Contour3D walls(const ExPolygon& floor_plate, const ExPolygon& ceiling, + double floor_z_mm, double ceiling_z_mm) { + using std::transform; using std::back_inserter; + + ExPolygon poly; + poly.contour.points = floor_plate.contour.points; + poly.holes.emplace_back(ceiling.contour); + auto& h = poly.holes.front(); + std::reverse(h.points.begin(), h.points.end()); + Polygons tri = triangulate(poly); + + Contour3D ret; + ret.points.reserve(tri.size() * 3); + + double fz = floor_z_mm; + double cz = ceiling_z_mm; + auto& rp = ret.points; + auto& rpi = ret.indices; + ret.indices.reserve(tri.size() * 3); + + coord_t idx = 0; + + auto hlines = h.lines(); + auto is_upper = [&hlines](const Point& p) { + return std::any_of(hlines.begin(), hlines.end(), + [&p](const Line& l) { + return l.distance_to(p) < mm(1e-6); + }); + }; + + std::for_each(tri.begin(), tri.end(), + [&rp, &rpi, &poly, &idx, is_upper, fz, cz](const Polygon& pp) + { + for(auto& p : pp.points) + if(is_upper(p)) + rp.emplace_back(unscale(x(p), y(p), mm(cz))); + else rp.emplace_back(unscale(x(p), y(p), mm(fz))); + + coord_t a = idx++, b = idx++, c = idx++; + if(fz > cz) rpi.emplace_back(c, b, a); + else rpi.emplace_back(a, b, c); + }); + + return ret; +} + +/// Offsetting with clipper and smoothing the edges into a curvature. +void offset(ExPolygon& sh, coord_t distance) { + using ClipperLib::ClipperOffset; + using ClipperLib::jtRound; + using ClipperLib::etClosedPolygon; + using ClipperLib::Paths; + using ClipperLib::Path; + + auto&& ctour = Slic3rMultiPoint_to_ClipperPath(sh.contour); + auto&& holes = Slic3rMultiPoints_to_ClipperPaths(sh.holes); + + // If the input is not at least a triangle, we can not do this algorithm + if(ctour.size() < 3 || + std::any_of(holes.begin(), holes.end(), + [](const Path& p) { return p.size() < 3; }) + ) { + BOOST_LOG_TRIVIAL(error) << "Invalid geometry for offsetting!"; + return; + } + + ClipperOffset offs; + offs.ArcTolerance = 0.01*mm(1); + Paths result; + offs.AddPath(ctour, jtRound, etClosedPolygon); + offs.AddPaths(holes, jtRound, etClosedPolygon); + offs.Execute(result, static_cast(distance)); + + // Offsetting reverts the orientation and also removes the last vertex + // so boost will not have a closed polygon. + + bool found_the_contour = false; + sh.holes.clear(); + for(auto& r : result) { + if(ClipperLib::Orientation(r)) { + // We don't like if the offsetting generates more than one contour + // but throwing would be an overkill. Instead, we should warn the + // caller about the inability to create correct geometries + if(!found_the_contour) { + auto rr = ClipperPath_to_Slic3rPolygon(r); + sh.contour.points.swap(rr.points); + found_the_contour = true; + } else { + BOOST_LOG_TRIVIAL(warning) + << "Warning: offsetting result is invalid!"; + } + } else { + // TODO If there are multiple contours we can't be sure which hole + // belongs to the first contour. (But in this case the situation is + // bad enough to let it go...) + sh.holes.emplace_back(ClipperPath_to_Slic3rPolygon(r)); + } + } +} + +/// Unification of polygons (with clipper) preserving holes as well. +ExPolygons unify(const ExPolygons& shapes) { + using ClipperLib::ptSubject; + + ExPolygons retv; + + bool closed = true; + bool valid = true; + + ClipperLib::Clipper clipper; + + for(auto& path : shapes) { + auto clipperpath = Slic3rMultiPoint_to_ClipperPath(path.contour); + + if(!clipperpath.empty()) + valid &= clipper.AddPath(clipperpath, ptSubject, closed); + + auto clipperholes = Slic3rMultiPoints_to_ClipperPaths(path.holes); + + for(auto& hole : clipperholes) { + if(!hole.empty()) + valid &= clipper.AddPath(hole, ptSubject, closed); + } + } + + if(!valid) BOOST_LOG_TRIVIAL(warning) << "Unification of invalid shapes!"; + + ClipperLib::PolyTree result; + clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNonZero); + + retv.reserve(static_cast(result.Total())); + + // Now we will recursively traverse the polygon tree and serialize it + // into an ExPolygon with holes. The polygon tree has the clipper-ish + // PolyTree structure which alternates its nodes as contours and holes + + // A "declaration" of function for traversing leafs which are holes + std::function processHole; + + // Process polygon which calls processHoles which than calls processPoly + // again until no leafs are left. + auto processPoly = [&retv, &processHole](ClipperLib::PolyNode *pptr) { + ExPolygon poly; + poly.contour.points = ClipperPath_to_Slic3rPolygon(pptr->Contour); + for(auto h : pptr->Childs) { processHole(h, poly); } + retv.push_back(poly); + }; + + // Body of the processHole function + processHole = [&processPoly](ClipperLib::PolyNode *pptr, ExPolygon& poly) + { + poly.holes.emplace_back(); + poly.holes.back().points = ClipperPath_to_Slic3rPolygon(pptr->Contour); + for(auto c : pptr->Childs) processPoly(c); + }; + + // Wrapper for traversing. + auto traverse = [&processPoly] (ClipperLib::PolyNode *node) + { + for(auto ch : node->Childs) { + processPoly(ch); + } + }; + + // Here is the actual traverse + traverse(&result); + + return retv; +} + +/// Only a debug function to generate top and bottom plates from a 2D shape. +/// It is not used in the algorithm directly. +inline Contour3D roofs(const ExPolygon& poly, coord_t z_distance) { + Polygons triangles = triangulate(poly); + + auto lower = convert(triangles, 0, false); + auto upper = convert(triangles, z_distance, true); + lower.merge(upper); + return lower; +} + +template +Contour3D round_edges(const ExPolygon& base_plate, + double radius_mm, + double degrees, + double ceilheight_mm, + bool dir, + ExP&& last_offset = ExP(), D&& last_height = D()) +{ + auto ob = base_plate; + auto ob_prev = ob; + double wh = ceilheight_mm, wh_prev = wh; + Contour3D curvedwalls; + + int steps = 15; // int(std::ceil(10*std::pow(radius_mm, 1.0/3))); + double stepx = radius_mm / steps; + coord_t s = dir? 1 : -1; + degrees = std::fmod(degrees, 180); + + if(degrees >= 90) { + for(int i = 1; i <= steps; ++i) { + ob = base_plate; + + double r2 = radius_mm * radius_mm; + double xx = i*stepx; + double x2 = xx*xx; + double stepy = std::sqrt(r2 - x2); + + offset(ob, s*mm(xx)); + wh = ceilheight_mm - radius_mm + stepy; + + Contour3D pwalls; + pwalls = walls(ob, ob_prev, wh, wh_prev); + + curvedwalls.merge(pwalls); + ob_prev = ob; + wh_prev = wh; + } + } + + double tox = radius_mm - radius_mm*std::sin(degrees * PI / 180); + int tos = int(tox / stepx); + + for(int i = 1; i <= tos; ++i) { + ob = base_plate; + + double r2 = radius_mm * radius_mm; + double xx = radius_mm - i*stepx; + double x2 = xx*xx; + double stepy = std::sqrt(r2 - x2); + offset(ob, s*mm(xx)); + wh = ceilheight_mm - radius_mm - stepy; + + Contour3D pwalls; + pwalls = walls(ob_prev, ob, wh_prev, wh); + + curvedwalls.merge(pwalls); + ob_prev = ob; + wh_prev = wh; + } + + last_offset = std::move(ob); + last_height = wh; + + return curvedwalls; +} + +/// Generating the concave part of the 3D pool with the bottom plate and the +/// side walls. +Contour3D inner_bed(const ExPolygon& poly, double depth_mm, + double begin_h_mm = 0) { + + Polygons triangles = triangulate(poly); + + coord_t depth = mm(depth_mm); + coord_t begin_h = mm(begin_h_mm); + + auto bottom = convert(triangles, -depth + begin_h, false); + auto lines = poly.lines(); + + // Generate outer walls + auto fp = [](const Point& p, Point::coord_type z) { + return unscale(x(p), y(p), z); + }; + + for(auto& l : lines) { + auto s = coord_t(bottom.points.size()); + + bottom.points.emplace_back(fp(l.a, -depth + begin_h)); + bottom.points.emplace_back(fp(l.b, -depth + begin_h)); + bottom.points.emplace_back(fp(l.a, begin_h)); + bottom.points.emplace_back(fp(l.b, begin_h)); + + bottom.indices.emplace_back(s + 3, s + 1, s); + bottom.indices.emplace_back(s + 2, s + 3, s); + } + + return bottom; +} + +inline Point centroid(Points& pp) { + Point c; + switch(pp.size()) { + case 0: break; + case 1: c = pp.front(); break; + case 2: c = (pp[0] + pp[1]) / 2; break; + default: { + auto MAX = std::numeric_limits::max(); + auto MIN = std::numeric_limits::min(); + Point min = {MAX, MAX}, max = {MIN, MIN}; + + for(auto& p : pp) { + if(p(0) < min(0)) min(0) = p(0); + if(p(1) < min(1)) min(1) = p(1); + if(p(0) > max(0)) max(0) = p(0); + if(p(1) > max(1)) max(1) = p(1); + } + c(0) = min(0) + (max(0) - min(0)) / 2; + c(1) = min(1) + (max(1) - min(1)) / 2; + + // TODO: fails for non convex cluster +// c = std::accumulate(pp.begin(), pp.end(), Point{0, 0}); +// x(c) /= coord_t(pp.size()); y(c) /= coord_t(pp.size()); + break; + } + } + + return c; +} + +inline Point centroid(const ExPolygon& poly) { + return poly.contour.centroid(); +} + +/// A fake concave hull that is constructed by connecting separate shapes +/// with explicit bridges. Bridges are generated from each shape's centroid +/// to the center of the "scene" which is the centroid calculated from the shape +/// centroids (a star is created...) +ExPolygons concave_hull(const ExPolygons& polys, double max_dist_mm = 50) +{ + namespace bgi = boost::geometry::index; + using SpatElement = std::pair; + using SpatIndex = bgi::rtree< SpatElement, bgi::rstar<16, 4> >; + + if(polys.empty()) return ExPolygons(); + + ExPolygons punion = unify(polys); // could be redundant + + if(punion.size() == 1) return punion; + + // We get the centroids of all the islands in the 2D slice + Points centroids; centroids.reserve(punion.size()); + std::transform(punion.begin(), punion.end(), std::back_inserter(centroids), + [](const ExPolygon& poly) { return centroid(poly); }); + + + SpatIndex boxindex; unsigned idx = 0; + std::for_each(punion.begin(), punion.end(), + [&boxindex, &idx](const ExPolygon& expo) { + BoundingBox bb(expo); + boxindex.insert(std::make_pair(bb, idx++)); + }); + + + // Centroid of the centroids of islands. This is where the additional + // connector sticks are routed. + Point cc = centroid(centroids); + + punion.reserve(punion.size() + centroids.size()); + + idx = 0; + std::transform(centroids.begin(), centroids.end(), + std::back_inserter(punion), + [&punion, &boxindex, cc, max_dist_mm, &idx](const Point& c) + { + + double dx = x(c) - x(cc), dy = y(c) - y(cc); + double l = std::sqrt(dx * dx + dy * dy); + double nx = dx / l, ny = dy / l; + double max_dist = mm(max_dist_mm); + + ExPolygon& expo = punion[idx++]; + BoundingBox querybb(expo); + + querybb.offset(max_dist); + std::vector result; + boxindex.query(bgi::intersects(querybb), std::back_inserter(result)); + if(result.size() <= 1) return ExPolygon(); + + ExPolygon r; + auto& ctour = r.contour.points; + + ctour.reserve(3); + ctour.emplace_back(cc); + + Point d(coord_t(mm(1)*nx), coord_t(mm(1)*ny)); + ctour.emplace_back(c + Point( -y(d), x(d) )); + ctour.emplace_back(c + Point( y(d), -x(d) )); + offset(r, mm(1)); + + return r; + }); + + punion = unify(punion); + + return punion; +} + +void base_plate(const TriangleMesh &mesh, ExPolygons &output, float h) +{ + TriangleMesh m = mesh; + TriangleMeshSlicer slicer(&m); + + TriangleMesh upper, lower; + slicer.cut(h, &upper, &lower); + output = lower.horizontal_projection(); + + for(auto& o : output) o = o.simplify(0.1/SCALING_FACTOR).front(); +} + +void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out, + const PoolConfig& cfg) +{ + double mdist = 2*(1.8*cfg.min_wall_thickness_mm + 4*cfg.edge_radius_mm) + + cfg.max_merge_distance_mm; + + auto concavehs = concave_hull(ground_layer, mdist); + for(ExPolygon& concaveh : concavehs) { + if(concaveh.contour.points.empty()) return; + concaveh.holes.clear(); + + const coord_t WALL_THICKNESS = mm(cfg.min_wall_thickness_mm); + + const coord_t WALL_DISTANCE = mm(2*cfg.edge_radius_mm) + + coord_t(0.8*WALL_THICKNESS); + + const coord_t HEIGHT = mm(cfg.min_wall_height_mm); + + auto outer_base = concaveh; + offset(outer_base, WALL_THICKNESS+WALL_DISTANCE); + auto inner_base = outer_base; + offset(inner_base, -WALL_THICKNESS); + inner_base.holes.clear(); outer_base.holes.clear(); + + ExPolygon top_poly; + top_poly.contour = outer_base.contour; + top_poly.holes.emplace_back(inner_base.contour); + auto& tph = top_poly.holes.back().points; + std::reverse(tph.begin(), tph.end()); + + Contour3D pool; + + ExPolygon ob = outer_base; double wh = 0; + + // now we will calculate the angle or portion of the circle from + // pi/2 that will connect perfectly with the bottom plate. + // this is a tangent point calculation problem and the equation can + // be found for example here: + // http://www.ambrsoft.com/TrigoCalc/Circles2/CirclePoint/CirclePointDistance.htm + // the y coordinate would be: + // y = cy + (r^2*py - r*px*sqrt(px^2 + py^2 - r^2) / (px^2 + py^2) + // where px and py are the coordinates of the point outside the circle + // cx and cy are the circle center, r is the radius + // to get the angle we use arcsin function and subtract 90 degrees then + // flip the sign to get the right input to the round_edge function. + double r = cfg.edge_radius_mm; + double cy = 0; + double cx = 0; + double px = cfg.min_wall_thickness_mm; + double py = r - cfg.min_wall_height_mm; + + double pxcx = px - cx; + double pycy = py - cy; + double b_2 = pxcx*pxcx + pycy*pycy; + double r_2 = r*r; + double D = std::sqrt(b_2 - r_2); + double vy = (r_2*pycy - r*pxcx*D) / b_2; + double phi = -(std::asin(vy/r) * 180 / PI - 90); + + auto curvedwalls = round_edges(ob, + r, + phi, // 170 degrees + 0, // z position of the input plane + true, + ob, wh); + + pool.merge(curvedwalls); + + ExPolygon ob_contr = ob; + ob_contr.holes.clear(); + + auto pwalls = walls(ob_contr, inner_base, wh, -cfg.min_wall_height_mm); + pool.merge(pwalls); + + Polygons top_triangles, bottom_triangles; + triangulate(top_poly, top_triangles); + triangulate(inner_base, bottom_triangles); + auto top_plate = convert(top_triangles, 0, false); + auto bottom_plate = convert(bottom_triangles, -HEIGHT, true); + + ob = inner_base; wh = 0; + // rounded edge generation for the inner bed + curvedwalls = round_edges(ob, + cfg.edge_radius_mm, + 90, // 90 degrees + 0, // z position of the input plane + false, + ob, wh); + pool.merge(curvedwalls); + + auto innerbed = inner_bed(ob, cfg.min_wall_height_mm/2 + wh, wh); + + pool.merge(top_plate); + pool.merge(bottom_plate); + pool.merge(innerbed); + + out.merge(mesh(pool)); + } +} + +} +} diff --git a/src/libslic3r/SLA/SLABasePool.hpp b/src/libslic3r/SLA/SLABasePool.hpp new file mode 100644 index 000000000..6abddf5ff --- /dev/null +++ b/src/libslic3r/SLA/SLABasePool.hpp @@ -0,0 +1,37 @@ +#ifndef SLABASEPOOL_HPP +#define SLABASEPOOL_HPP + +#include + +namespace Slic3r { + +class ExPolygon; +class TriangleMesh; + +namespace sla { + +using ExPolygons = std::vector; + +/// Calculate the polygon representing the silhouette from the specified height +void base_plate(const TriangleMesh& mesh, + ExPolygons& output, + float height = 0.1f); + +struct PoolConfig { + double min_wall_thickness_mm = 2; + double min_wall_height_mm = 5; + double max_merge_distance_mm = 50; + double edge_radius_mm = 1; +}; + +/// Calculate the pool for the mesh for SLA printing +void create_base_pool(const ExPolygons& base_plate, + TriangleMesh& output_mesh, + const PoolConfig& = PoolConfig() + ); + +} + +} + +#endif // SLABASEPOOL_HPP diff --git a/src/libslic3r/SLA/SLABoilerPlate.hpp b/src/libslic3r/SLA/SLABoilerPlate.hpp new file mode 100644 index 000000000..c15137ac5 --- /dev/null +++ b/src/libslic3r/SLA/SLABoilerPlate.hpp @@ -0,0 +1,92 @@ +#ifndef SLABOILERPLATE_HPP +#define SLABOILERPLATE_HPP + +#include +#include +#include + +#include "ExPolygon.hpp" +#include "TriangleMesh.hpp" + +namespace Slic3r { +namespace sla { + +using coord_t = Point::coord_type; + +/// get the scaled clipper units for a millimeter value +inline coord_t mm(double v) { return coord_t(v/SCALING_FACTOR); } + +/// Get x and y coordinates (because we are eigenizing...) +inline coord_t x(const Point& p) { return p(0); } +inline coord_t y(const Point& p) { return p(1); } +inline coord_t& x(Point& p) { return p(0); } +inline coord_t& y(Point& p) { return p(1); } + +inline coordf_t x(const Vec3d& p) { return p(0); } +inline coordf_t y(const Vec3d& p) { return p(1); } +inline coordf_t z(const Vec3d& p) { return p(2); } +inline coordf_t& x(Vec3d& p) { return p(0); } +inline coordf_t& y(Vec3d& p) { return p(1); } +inline coordf_t& z(Vec3d& p) { return p(2); } + +inline coord_t& x(Vec3crd& p) { return p(0); } +inline coord_t& y(Vec3crd& p) { return p(1); } +inline coord_t& z(Vec3crd& p) { return p(2); } +inline coord_t x(const Vec3crd& p) { return p(0); } +inline coord_t y(const Vec3crd& p) { return p(1); } +inline coord_t z(const Vec3crd& p) { return p(2); } + +inline void triangulate(const ExPolygon& expoly, Polygons& triangles) { + expoly.triangulate_p2t(&triangles); +} + +inline Polygons triangulate(const ExPolygon& expoly) { + Polygons tri; triangulate(expoly, tri); return tri; +} + +using Indices = std::vector; + +/// Intermediate struct for a 3D mesh +struct Contour3D { + Pointf3s points; + Indices indices; + + void merge(const Contour3D& ctr) { + auto s3 = coord_t(points.size()); + auto s = coord_t(indices.size()); + + points.insert(points.end(), ctr.points.begin(), ctr.points.end()); + indices.insert(indices.end(), ctr.indices.begin(), ctr.indices.end()); + + for(auto n = s; n < indices.size(); n++) { + auto& idx = indices[n]; x(idx) += s3; y(idx) += s3; z(idx) += s3; + } + } +}; + +struct EigenMesh3D { + Eigen::MatrixXd V; + Eigen::MatrixXi F; +}; + +using PointSet = Eigen::MatrixXd; +using ClusterEl = std::vector; +using ClusteredPoints = std::vector; + +/// Convert the triangulation output to an intermediate mesh. +Contour3D convert(const Polygons& triangles, coord_t z, bool dir); + +/// Mesh from an existing contour. +inline TriangleMesh mesh(const Contour3D& ctour) { + return {ctour.points, ctour.indices}; +} + +/// Mesh from an evaporating 3D contour +inline TriangleMesh mesh(Contour3D&& ctour) { + return {std::move(ctour.points), std::move(ctour.indices)}; +} + +} +} + +#endif // SLABOILERPLATE_HPP diff --git a/src/libslic3r/SLA/SLABoostAdapter.hpp b/src/libslic3r/SLA/SLABoostAdapter.hpp new file mode 100644 index 000000000..a31bfe60c --- /dev/null +++ b/src/libslic3r/SLA/SLABoostAdapter.hpp @@ -0,0 +1,132 @@ +#ifndef SLABOOSTADAPTER_HPP +#define SLABOOSTADAPTER_HPP + +#include "SLA/SLABoilerPlate.hpp" +#include + +namespace boost { +namespace geometry { +namespace traits { + +/* ************************************************************************** */ +/* Point concept adaptation ************************************************* */ +/* ************************************************************************** */ + +template<> struct tag { + using type = point_tag; +}; + +template<> struct coordinate_type { + using type = coord_t; +}; + +template<> struct coordinate_system { + using type = cs::cartesian; +}; + +template<> struct dimension: boost::mpl::int_<2> {}; + +template struct access { + static inline coord_t get(Slic3r::Point const& a) { + return a(d); + } + + static inline void set(Slic3r::Point& a, coord_t const& value) { + a(d) = value; + } +}; + +// For Vec2d /////////////////////////////////////////////////////////////////// + +template<> struct tag { + using type = point_tag; +}; + +template<> struct coordinate_type { + using type = double; +}; + +template<> struct coordinate_system { + using type = cs::cartesian; +}; + +template<> struct dimension: boost::mpl::int_<2> {}; + +template struct access { + static inline double get(Slic3r::Vec2d const& a) { + return a(d); + } + + static inline void set(Slic3r::Vec2d& a, double const& value) { + a(d) = value; + } +}; + +// For Vec3d /////////////////////////////////////////////////////////////////// + +template<> struct tag { + using type = point_tag; +}; + +template<> struct coordinate_type { + using type = double; +}; + +template<> struct coordinate_system { + using type = cs::cartesian; +}; + +template<> struct dimension: boost::mpl::int_<3> {}; + +template struct access { + static inline double get(Slic3r::Vec3d const& a) { + return a(d); + } + + static inline void set(Slic3r::Vec3d& a, double const& value) { + a(d) = value; + } +}; + +/* ************************************************************************** */ +/* Box concept adaptation *************************************************** */ +/* ************************************************************************** */ + +template<> struct tag { + using type = box_tag; +}; + +template<> struct point_type { + using type = Slic3r::Point; +}; + +template +struct indexed_access { + static inline coord_t get(Slic3r::BoundingBox const& box) { + return box.min(d); + } + static inline void set(Slic3r::BoundingBox &box, coord_t const& coord) { + box.min(d) = coord; + } +}; + +template +struct indexed_access { + static inline coord_t get(Slic3r::BoundingBox const& box) { + return box.max(d); + } + static inline void set(Slic3r::BoundingBox &box, coord_t const& coord) { + box.max(d) = coord; + } +}; + +} +} + +template<> struct range_value> { + using type = Slic3r::Vec2d; +}; + +} + +#endif // SLABOOSTADAPTER_HPP diff --git a/src/libslic3r/SLA/SLARotfinder.cpp b/src/libslic3r/SLA/SLARotfinder.cpp new file mode 100644 index 000000000..f2dab05ac --- /dev/null +++ b/src/libslic3r/SLA/SLARotfinder.cpp @@ -0,0 +1,156 @@ +#include +#include + +#include +#include "SLABoilerPlate.hpp" +#include "SLARotfinder.hpp" +#include "Model.hpp" + +namespace Slic3r { +namespace sla { + +EigenMesh3D to_eigenmesh(const TriangleMesh& tmesh) { + + const stl_file& stl = tmesh.stl; + + EigenMesh3D outmesh; + auto& V = outmesh.V; + auto& F = outmesh.F; + + V.resize(3*stl.stats.number_of_facets, 3); + F.resize(stl.stats.number_of_facets, 3); + for (unsigned int i=0; ivertex[0](0); V(3*i+0, 1) = + facet->vertex[0](1); V(3*i+0, 2) = facet->vertex[0](2); + V(3*i+1, 0) = facet->vertex[1](0); V(3*i+1, 1) = + facet->vertex[1](1); V(3*i+1, 2) = facet->vertex[1](2); + V(3*i+2, 0) = facet->vertex[2](0); V(3*i+2, 1) = + facet->vertex[2](1); V(3*i+2, 2) = facet->vertex[2](2); + + F(i, 0) = 3*i+0; + F(i, 1) = 3*i+1; + F(i, 2) = 3*i+2; + } + + return outmesh; +} + +inline EigenMesh3D to_eigenmesh(const ModelObject& modelobj) { + return to_eigenmesh(modelobj.raw_mesh()); +} + +std::array find_best_rotation(const ModelObject& modelobj, + float accuracy, + std::function statuscb, + std::function stopcond) +{ + using libnest2d::opt::Method; + using libnest2d::opt::bound; + using libnest2d::opt::Optimizer; + using libnest2d::opt::TOptimizer; + using libnest2d::opt::StopCriteria; + using Quaternion = Eigen::Quaternion; + + static const unsigned MAX_TRIES = 100000; + + // return value + std::array rot; + + // We will use only one instance of this converted mesh to examine different + // rotations + EigenMesh3D emesh = to_eigenmesh(modelobj); + + // For current iteration number + unsigned status = 0; + + // The maximum number of iterations + auto max_tries = unsigned(accuracy * MAX_TRIES); + + // call status callback with zero, because we are at the start + statuscb(status); + + // So this is the object function which is called by the solver many times + // It has to yield a single value representing the current score. We will + // call the status callback in each iteration but the actual value may be + // the same for subsequent iterations (status goes from 0 to 100 but + // iterations can be many more) + auto objfunc = [&emesh, &status, &statuscb, max_tries] + (double rx, double ry, double rz) + { + EigenMesh3D& m = emesh; + + // prepare the rotation transformation + Transform3d rt = Transform3d::Identity(); + + rt.rotate(Eigen::AngleAxisd(rz, Vec3d::UnitZ())); + rt.rotate(Eigen::AngleAxisd(ry, Vec3d::UnitY())); + rt.rotate(Eigen::AngleAxisd(rx, Vec3d::UnitX())); + + double score = 0; + + // For all triangles we calculate the normal and sum up the dot product + // (a scalar indicating how much are two vectors aligned) with each axis + // this will result in a value that is greater if a normal is aligned + // with all axes. If the normal is aligned than the triangle itself is + // orthogonal to the axes and that is good for print quality. + + // TODO: some applications optimize for minimum z-axis cross section + // area. The current function is only an example of how to optimize. + + // Later we can add more criteria like the number of overhangs, etc... + for(int i = 0; i < m.F.rows(); i++) { + auto idx = m.F.row(i); + + Vec3d p1 = m.V.row(idx(0)); + Vec3d p2 = m.V.row(idx(1)); + Vec3d p3 = m.V.row(idx(2)); + + Eigen::Vector3d U = p2 - p1; + Eigen::Vector3d V = p3 - p1; + + // So this is the normal + auto n = U.cross(V).normalized(); + + // rotate the normal with the current rotation given by the solver + n = rt * n; + + // We should score against the alignment with the reference planes + score += std::abs(n.dot(Vec3d::UnitX())); + score += std::abs(n.dot(Vec3d::UnitY())); + score += std::abs(n.dot(Vec3d::UnitZ())); + } + + // report status + statuscb( unsigned(++status * 100.0/max_tries) ); + + return score; + }; + + // Firing up the genetic optimizer. For now it uses the nlopt library. + StopCriteria stc; + stc.max_iterations = max_tries; + stc.relative_score_difference = 1e-3; + stc.stop_condition = stopcond; // stop when stopcond returns true + TOptimizer solver(stc); + + // We are searching rotations around the three axes x, y, z. Thus the + // problem becomes a 3 dimensional optimization task. + // We can specify the bounds for a dimension in the following way: + auto b = bound(-PI/2, PI/2); + + // Now we start the optimization process with initial angles (0, 0, 0) + auto result = solver.optimize_max(objfunc, + libnest2d::opt::initvals(0.0, 0.0, 0.0), + b, b, b); + + // Save the result and fck off + rot[0] = std::get<0>(result.optimum); + rot[1] = std::get<1>(result.optimum); + rot[2] = std::get<2>(result.optimum); + + return rot; +} + +} +} diff --git a/src/libslic3r/SLA/SLARotfinder.hpp b/src/libslic3r/SLA/SLARotfinder.hpp new file mode 100644 index 000000000..bbfab670a --- /dev/null +++ b/src/libslic3r/SLA/SLARotfinder.hpp @@ -0,0 +1,38 @@ +#ifndef SLAROTFINDER_HPP +#define SLAROTFINDER_HPP + +#include +#include + +namespace Slic3r { + +class ModelObject; + +namespace sla { + +/** + * The function should find the best rotation for SLA upside down printing. + * + * @param modelobj The model object representing the 3d mesh. + * @param accuracy The optimization accuracy from 0.0f to 1.0f. Currently, + * the nlopt genetic optimizer is used and the number of iterations is + * accuracy * 100000. This can change in the future. + * @param statuscb A status indicator callback called with the unsigned + * argument spanning from 0 to 100. May not reach 100 if the optimization finds + * an optimum before max iterations are reached. + * @param stopcond A function that if returns true, the search process will be + * terminated and the best solution found will be returned. + * + * @return Returns the rotations around each axis (x, y, z) + */ +std::array find_best_rotation( + const ModelObject& modelobj, + float accuracy = .0f, + std::function statuscb = [] (unsigned) {}, + std::function stopcond = [] () { return false; } + ); + +} +} + +#endif // SLAROTFINDER_HPP diff --git a/src/libslic3r/SLA/SLASpatIndex.hpp b/src/libslic3r/SLA/SLASpatIndex.hpp new file mode 100644 index 000000000..dfe5e03a0 --- /dev/null +++ b/src/libslic3r/SLA/SLASpatIndex.hpp @@ -0,0 +1,112 @@ +#ifndef SPATINDEX_HPP +#define SPATINDEX_HPP + +#include +#include +#include + +#include + +namespace Slic3r { +namespace sla { + +typedef Eigen::Matrix Vec3d; +using SpatElement = std::pair; + +/** + * This class is intended for enhancing range based for loops with indexing. + * So instead of this: + * { int i = 0; for(auto c : container) { process(c, i); ++i; } + * + * you can use this: + * for(auto ic : container) process(ic.value, ic.index); + */ +template class Enumerable { + Container&& m; + using C = typename std::remove_reference::type; + using CC = typename std::remove_const::type; + + template struct get_iterator {}; + template<> struct get_iterator { using type = typename CC::iterator; }; + template<> struct get_iterator { + using type = typename CC::const_iterator; + }; + + template struct _Pair { + Vref value; + size_t index; + _Pair(Vref v, size_t i) : value(v), index(i) {} + operator Vref() { return value; } + }; + + template + class _iterator { + Cit start; + Cit it; + using Pair = _Pair::reference>; + public: + _iterator(Cit b, Cit i): start(b), it(i) {} + _iterator& operator++() { ++it; return *this;} + _iterator operator++(int) { auto tmp = it; ++it; return tmp;} + + bool operator!=(_iterator other) { return it != other.it; } + Pair operator*() { return Pair(*it, it - start); } + using value_type = typename Enumerable::value_type; + }; + +public: + + Enumerable(Container&& c): m(c) {} + + using value_type = typename CC::value_type; + + using iterator = _iterator::type>; + using const_iterator = _iterator; + + iterator begin() { return iterator(m.begin(), m.begin()); } + iterator end() { return iterator(m.begin(), m.end()); } + const_iterator begin() const {return const_iterator(m.cbegin(), m.cbegin());} + const_iterator end() const { return const_iterator(m.cbegin(), m.cend());} + +}; + +template inline Enumerable enumerate(C&& c) { + return Enumerable(std::forward(c)); +} + +class SpatIndex { + class Impl; + + // We use Pimpl because it takes a long time to compile boost headers which + // is the engine of this class. We include it only in the cpp file. + std::unique_ptr m_impl; +public: + + SpatIndex(); + ~SpatIndex(); + + SpatIndex(const SpatIndex&); + SpatIndex(SpatIndex&&); + SpatIndex& operator=(const SpatIndex&); + SpatIndex& operator=(SpatIndex&&); + + void insert(const SpatElement&); + bool remove(const SpatElement&); + + inline void insert(const Vec3d& v, unsigned idx) + { + insert(std::make_pair(v, unsigned(idx))); + } + + std::vector query(std::function); + std::vector nearest(const Vec3d&, unsigned k); + + // For testing + size_t size() const; + bool empty() const { return size() == 0; } +}; + +} +} + +#endif // SPATINDEX_HPP diff --git a/src/libslic3r/SLA/SLASupportTree.cpp b/src/libslic3r/SLA/SLASupportTree.cpp new file mode 100644 index 000000000..1cf3194b3 --- /dev/null +++ b/src/libslic3r/SLA/SLASupportTree.cpp @@ -0,0 +1,1665 @@ +/** + * In this file we will implement the automatic SLA support tree generation. + * + */ + +#include +#include "SLASupportTree.hpp" +#include "SLABoilerPlate.hpp" +#include "SLASpatIndex.hpp" +#include "SLABasePool.hpp" +#include + +#include "Model.hpp" + +/** + * Terminology: + * + * Support point: + * The point on the model surface that needs support. + * + * Pillar: + * A thick column that spans from a support point to the ground and has + * a thick cone shaped base where it touches the ground. + * + * Ground facing support point: + * A support point that can be directly connected with the ground with a pillar + * that does not collide or cut through the model. + * + * Non ground facing support point: + * A support point that cannot be directly connected with the ground (only with + * the model surface). + * + * Head: + * The pinhead that connects to the model surface with the sharp end end + * to a pillar or bridge stick with the dull end. + * + * Headless support point: + * A support point on the model surface for which there is not enough place for + * the head. It is either in a hole or there is some barrier that would collide + * with the head geometry. The headless support point can be ground facing and + * non ground facing as well. + * + * Bridge: + * A stick that connects two pillars or a head with a pillar. + * + * Junction: + * A small ball in the intersection of two or more sticks (pillar, bridge, ...) + * + * CompactBridge: + * A bridge that connects a headless support point with the model surface or a + * nearby pillar. + */ + +namespace Slic3r { +namespace sla { + +using Coordf = double; +using Portion = std::tuple; + +inline Portion make_portion(double a, double b) { + return std::make_tuple(a, b); +} + +template double distance(const Vec& pp1, const Vec& pp2) { + auto p = pp2 - pp1; + return distance(p); +} + +template double distance(const Vec& p) { + return std::sqrt(p.transpose() * p); +} + +/// The horizontally projected 2D boundary of the model as individual line +/// segments. This can be used later to create a spatial index of line segments +/// and query for possible pillar positions for non-ground facing support points +std::vector> model_boundary(const EigenMesh3D& emesh, + double offs = 0.01); + +Contour3D sphere(double rho, Portion portion = make_portion(0.0, 2.0*PI), + double fa=(2*PI/360)) { + + Contour3D ret; + + // prohibit close to zero radius + if(rho <= 1e-6 && rho >= -1e-6) return ret; + + auto& vertices = ret.points; + auto& facets = ret.indices; + + // Algorithm: + // Add points one-by-one to the sphere grid and form facets using relative + // coordinates. Sphere is composed effectively of a mesh of stacked circles. + + // adjust via rounding to get an even multiple for any provided angle. + double angle = (2*PI / floor(2*PI / fa)); + + // Ring to be scaled to generate the steps of the sphere + std::vector ring; + + for (double i = 0; i < 2*PI; i+=angle) ring.emplace_back(i); + + const auto sbegin = size_t(2*std::get<0>(portion)/angle); + const auto send = size_t(2*std::get<1>(portion)/angle); + + const size_t steps = ring.size(); + const double increment = (double)(1.0 / (double)steps); + + // special case: first ring connects to 0,0,0 + // insert and form facets. + if(sbegin == 0) + vertices.emplace_back(Vec3d(0.0, 0.0, -rho + increment*sbegin*2.0*rho)); + + auto id = coord_t(vertices.size()); + for (size_t i = 0; i < ring.size(); i++) { + // Fixed scaling + const double z = -rho + increment*rho*2.0 * (sbegin + 1.0); + // radius of the circle for this step. + const double r = sqrt(abs(rho*rho - z*z)); + Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r); + vertices.emplace_back(Vec3d(b(0), b(1), z)); + + if(sbegin == 0) + facets.emplace_back((i == 0) ? Vec3crd(coord_t(ring.size()), 0, 1) : + Vec3crd(id - 1, 0, id)); + ++ id; + } + + // General case: insert and form facets for each step, + // joining it to the ring below it. + for (size_t s = sbegin + 2; s < send - 1; s++) { + const double z = -rho + increment*(double)s*2.0*rho; + const double r = sqrt(abs(rho*rho - z*z)); + + for (size_t i = 0; i < ring.size(); i++) { + Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r); + vertices.emplace_back(Vec3d(b(0), b(1), z)); + auto id_ringsize = coord_t(id - ring.size()); + if (i == 0) { + // wrap around + facets.emplace_back(Vec3crd(id - 1, id, + id + coord_t(ring.size() - 1))); + facets.emplace_back(Vec3crd(id - 1, id_ringsize, id)); + } else { + facets.emplace_back(Vec3crd(id_ringsize - 1, id_ringsize, id)); + facets.emplace_back(Vec3crd(id - 1, id_ringsize - 1, id)); + } + id++; + } + } + + // special case: last ring connects to 0,0,rho*2.0 + // only form facets. + if(send >= size_t(2*PI / angle)) { + vertices.emplace_back(Vec3d(0.0, 0.0, -rho + increment*send*2.0*rho)); + for (size_t i = 0; i < ring.size(); i++) { + auto id_ringsize = coord_t(id - ring.size()); + if (i == 0) { + // third vertex is on the other side of the ring. + facets.emplace_back(Vec3crd(id - 1, id_ringsize, id)); + } else { + auto ci = coord_t(id_ringsize + i); + facets.emplace_back(Vec3crd(ci - 1, ci, id)); + } + } + } + id++; + + return ret; +} + +Contour3D cylinder(double r, double h, size_t steps) { + Contour3D ret; + + auto& points = ret.points; + auto& indices = ret.indices; + points.reserve(2*steps); + double a = 2*PI/steps; + + Vec3d jp = {0, 0, 0}; + Vec3d endp = {0, 0, h}; + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double ex = endp(X) + r*std::cos(phi); + double ey = endp(Y) + r*std::sin(phi); + points.emplace_back(ex, ey, endp(Z)); + } + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double x = jp(X) + r*std::cos(phi); + double y = jp(Y) + r*std::sin(phi); + points.emplace_back(x, y, jp(Z)); + } + + indices.reserve(2*steps); + auto offs = steps; + for(int i = 0; i < steps - 1; ++i) { + indices.emplace_back(i, i + offs, offs + i + 1); + indices.emplace_back(i, offs + i + 1, i + 1); + } + + auto last = steps - 1; + indices.emplace_back(0, last, offs); + indices.emplace_back(last, offs + last, offs); + + return ret; +} + +struct Head { + Contour3D mesh; + + size_t steps = 45; + Vec3d dir = {0, 0, -1}; + Vec3d tr = {0, 0, 0}; + + double r_back_mm = 1; + double r_pin_mm = 0.5; + double width_mm = 2; + + // For identification purposes. This will be used as the index into the + // container holding the head structures. See SLASupportTree::Impl + long id = -1; + + // If there is a pillar connecting to this head, then the id will be set. + long pillar_id = -1; + + Head(double r_big_mm, + double r_small_mm, + double length_mm, + Vec3d direction = {0, 0, -1}, // direction (normal to the dull end ) + Vec3d offset = {0, 0, 0}, // displacement + const size_t circlesteps = 45): + steps(circlesteps), dir(direction), tr(offset), + r_back_mm(r_big_mm), r_pin_mm(r_small_mm), width_mm(length_mm) + { + + // We create two spheres which will be connected with a robe that fits + // both circles perfectly. + + // Set up the model detail level + const double detail = 2*PI/steps; + + // We don't generate whole circles. Instead, we generate only the + // portions which are visible (not covered by the robe) To know the + // exact portion of the bottom and top circles we need to use some + // rules of tangent circles from which we can derive (using simple + // triangles the following relations: + + // The height of the whole mesh + const double h = r_big_mm + r_small_mm + width_mm; + double phi = PI/2 - std::acos( (r_big_mm - r_small_mm) / h ); + + // To generate a whole circle we would pass a portion of (0, Pi) + // To generate only a half horizontal circle we can pass (0, Pi/2) + // The calculated phi is an offset to the half circles needed to smooth + // the transition from the circle to the robe geometry + + auto&& s1 = sphere(r_big_mm, make_portion(0, PI/2 + phi), detail); + auto&& s2 = sphere(r_small_mm, make_portion(PI/2 + phi, PI), detail); + + for(auto& p : s2.points) z(p) += h; + + mesh.merge(s1); + mesh.merge(s2); + + for(size_t idx1 = s1.points.size() - steps, idx2 = s1.points.size(); + idx1 < s1.points.size() - 1; + idx1++, idx2++) + { + coord_t i1s1 = coord_t(idx1), i1s2 = coord_t(idx2); + coord_t i2s1 = i1s1 + 1, i2s2 = i1s2 + 1; + + mesh.indices.emplace_back(i1s1, i2s1, i2s2); + mesh.indices.emplace_back(i1s1, i2s2, i1s2); + } + + auto i1s1 = coord_t(s1.points.size()) - steps; + auto i2s1 = coord_t(s1.points.size()) - 1; + auto i1s2 = coord_t(s1.points.size()); + auto i2s2 = coord_t(s1.points.size()) + steps - 1; + + mesh.indices.emplace_back(i2s2, i2s1, i1s1); + mesh.indices.emplace_back(i1s2, i2s2, i1s1); + + // To simplify further processing, we translate the mesh so that the + // last vertex of the pointing sphere (the pinpoint) will be at (0,0,0) + for(auto& p : mesh.points) { z(p) -= (h + 0.5 * r_small_mm); } + } + + void transform() + { + using Quaternion = Eigen::Quaternion; + + // We rotate the head to the specified direction The head's pointing + // side is facing upwards so this means that it would hold a support + // point with a normal pointing straight down. This is the reason of + // the -1 z coordinate + auto quatern = Quaternion::FromTwoVectors(Vec3d{0, 0, -1}, dir); + + for(auto& p : mesh.points) p = quatern * p + tr; + } + + double fullwidth() const { + return 1.5 * r_pin_mm + width_mm + 2*r_back_mm; + } + + Vec3d junction_point() const { + return tr + ( 1.5 * r_pin_mm + width_mm + r_back_mm)*dir; + } + + double request_pillar_radius(double radius) const { + const double rmax = r_back_mm /* * 0.65*/ ; + return radius > 0 && radius < rmax ? radius : rmax; + } +}; + +struct Junction { + Contour3D mesh; + double r = 1; + size_t steps = 45; + Vec3d pos; + + long id = -1; + + Junction(const Vec3d& tr, double r_mm, size_t stepnum = 45): + r(r_mm), steps(stepnum), pos(tr) + { + mesh = sphere(r_mm, make_portion(0, PI), 2*PI/steps); + for(auto& p : mesh.points) p += tr; + } +}; + +struct Pillar { + Contour3D mesh; + Contour3D base; + double r = 1; + size_t steps = 0; + Vec3d endpoint; + + long id = -1; + + // If the pillar connects to a head, this is the id of that head + bool starts_from_head = true; // Could start from a junction as well + long start_junction_id = -1; + + Pillar(const Vec3d& jp, const Vec3d& endp, + double radius = 1, size_t st = 45): + r(radius), steps(st), endpoint(endp), starts_from_head(false) + { + auto& points = mesh.points; + auto& indices = mesh.indices; + points.reserve(2*steps); + double a = 2*PI/steps; + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double x = jp(X) + r*std::cos(phi); + double y = jp(Y) + r*std::sin(phi); + points.emplace_back(x, y, jp(Z)); + } + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double ex = endp(X) + r*std::cos(phi); + double ey = endp(Y) + r*std::sin(phi); + points.emplace_back(ex, ey, endp(Z)); + } + + indices.reserve(2*steps); + auto offs = steps; + for(int i = 0; i < steps - 1; ++i) { + indices.emplace_back(i, i + offs, offs + i + 1); + indices.emplace_back(i, offs + i + 1, i + 1); + } + + auto last = steps - 1; + indices.emplace_back(0, last, offs); + indices.emplace_back(last, offs + last, offs); + } + + Pillar(const Junction& junc, const Vec3d& endp): + Pillar(junc.pos, endp, junc.r, junc.steps){} + + Pillar(const Head& head, const Vec3d& endp, double radius = 1): + Pillar(head.junction_point(), endp, head.request_pillar_radius(radius), + head.steps) + { + } + + void add_base(double height = 3, double radius = 2) { + if(height <= 0) return; + + if(radius < r ) radius = r; + + double a = 2*PI/steps; + double z = endpoint(2) + height; + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double x = endpoint(0) + r*std::cos(phi); + double y = endpoint(1) + r*std::sin(phi); + base.points.emplace_back(x, y, z); + } + + for(int i = 0; i < steps; ++i) { + double phi = i*a; + double x = endpoint(0) + radius*std::cos(phi); + double y = endpoint(1) + radius*std::sin(phi); + base.points.emplace_back(x, y, z - height); + } + + auto ep = endpoint; ep(2) += height; + base.points.emplace_back(endpoint); + base.points.emplace_back(ep); + + auto& indices = base.indices; + auto hcenter = base.points.size() - 1; + auto lcenter = base.points.size() - 2; + auto offs = steps; + for(int i = 0; i < steps - 1; ++i) { + indices.emplace_back(i, i + offs, offs + i + 1); + indices.emplace_back(i, offs + i + 1, i + 1); + indices.emplace_back(i, i + 1, hcenter); + indices.emplace_back(lcenter, offs + i + 1, offs + i); + } + + auto last = steps - 1; + indices.emplace_back(0, last, offs); + indices.emplace_back(last, offs + last, offs); + indices.emplace_back(hcenter, last, 0); + indices.emplace_back(offs, offs + last, lcenter); + + } + + bool has_base() const { return !base.points.empty(); } +}; + +// A Bridge between two pillars (with junction endpoints) +struct Bridge { + Contour3D mesh; + double r = 0.8; + + long id = -1; + long start_jid = -1; + long end_jid = -1; + + // We should reduce the radius a tiny bit to help the convex hull algorithm + Bridge(const Vec3d& j1, const Vec3d& j2, + double r_mm = 0.8, size_t steps = 45): + r(r_mm) + { + using Quaternion = Eigen::Quaternion; + Vec3d dir = (j2 - j1).normalized(); + double d = distance(j2, j1); + + mesh = cylinder(r, d, steps); + + auto quater = Quaternion::FromTwoVectors(Vec3d{0,0,1}, dir); + for(auto& p : mesh.points) p = quater * p + j1; + } + + Bridge(const Junction& j1, const Junction& j2, double r_mm = 0.8): + Bridge(j1.pos, j2.pos, r_mm, j1.steps) {} + + Bridge(const Junction& j, const Pillar& cl) {} + +}; + +// A bridge that spans from model surface to model surface with small connecting +// edges on the endpoints. Used for headless support points. +struct CompactBridge { + Contour3D mesh; + long id = -1; + + CompactBridge(const Vec3d& sp, + const Vec3d& ep, + const Vec3d& n, + double r, + size_t steps = 45) + { + Vec3d startp = sp + r * n; + Vec3d dir = (ep - startp).normalized(); + Vec3d endp = ep - r * dir; + + Bridge br(startp, endp, r, steps); + mesh.merge(br.mesh); + + // now add the pins + double fa = 2*PI/steps; + auto upperball = sphere(r, Portion{PI / 2 - fa, PI}, fa); + for(auto& p : upperball.points) p += startp; + + auto lowerball = sphere(r, Portion{0, PI/2 + 2*fa}, fa); + for(auto& p : lowerball.points) p += endp; + + mesh.merge(upperball); + mesh.merge(lowerball); + } +}; + +EigenMesh3D to_eigenmesh(const Contour3D& cntr) { + EigenMesh3D emesh; + + auto& V = emesh.V; + auto& F = emesh.F; + + V.resize(cntr.points.size(), 3); + F.resize(cntr.indices.size(), 3); + + for (int i = 0; i < V.rows(); ++i) { + V.row(i) = cntr.points[i]; + F.row(i) = cntr.indices[i]; + } + + return emesh; +} + +void create_head(TriangleMesh& out, double r1_mm, double r2_mm, double width_mm) +{ + Head head(r1_mm, r2_mm, width_mm, {0, std::sqrt(0.5), -std::sqrt(0.5)}, + {0, 0, 30}); + out.merge(mesh(head.mesh)); + + Pillar cst(head, {0, 0, 0}); + cst.add_base(); + + out.merge(mesh(cst.mesh)); + out.merge(mesh(cst.base)); +} + +// The minimum distance for two support points to remain valid. +static const double /*constexpr*/ D_SP = 0.1; + +enum { // For indexing Eigen vectors as v(X), v(Y), v(Z) instead of numbers + X, Y, Z +}; + +EigenMesh3D to_eigenmesh(const Model& model) { + TriangleMesh combined_mesh; + + for(ModelObject *o : model.objects) { + TriangleMesh tmp = o->raw_mesh(); + for(ModelInstance * inst: o->instances) { + TriangleMesh ttmp(tmp); + inst->transform_mesh(&ttmp); + combined_mesh.merge(ttmp); + } + } + + const stl_file& stl = combined_mesh.stl; + + EigenMesh3D outmesh; + auto& V = outmesh.V; + auto& F = outmesh.F; + + V.resize(3*stl.stats.number_of_facets, 3); + F.resize(stl.stats.number_of_facets, 3); + for (unsigned int i=0; ivertex[0](0); V(3*i+0, 1) = + facet->vertex[0](1); V(3*i+0, 2) = facet->vertex[0](2); + V(3*i+1, 0) = facet->vertex[1](0); V(3*i+1, 1) = + facet->vertex[1](1); V(3*i+1, 2) = facet->vertex[1](2); + V(3*i+2, 0) = facet->vertex[2](0); V(3*i+2, 1) = + facet->vertex[2](1); V(3*i+2, 2) = facet->vertex[2](2); + + F(i, 0) = 3*i+0; + F(i, 1) = 3*i+1; + F(i, 2) = 3*i+2; + } + + return outmesh; +} + +Vec3d model_coord(const ModelInstance& object, const Vec3f& mesh_coord) { + return object.transform_vector(mesh_coord.cast()); +} + +PointSet support_points(const Model& model) { + size_t sum = 0; + for(auto *o : model.objects) + sum += o->instances.size() * o->sla_support_points.size(); + + PointSet ret(sum, 3); + + for(ModelObject *o : model.objects) + for(ModelInstance *inst : o->instances) { + int i = 0; + for(Vec3f& msource : o->sla_support_points) { + ret.row(i++) = model_coord(*inst, msource); + } + } + + return ret; +} + +double ray_mesh_intersect(const Vec3d& s, + const Vec3d& dir, + const EigenMesh3D& m); + +PointSet normals(const PointSet& points, const EigenMesh3D& mesh); + +Vec2d to_vec2(const Vec3d& v3) { + return {v3(X), v3(Y)}; +} + +bool operator==(const SpatElement& e1, const SpatElement& e2) { + return e1.second == e2.second; +} + +// Clustering a set of points by the given criteria +ClusteredPoints cluster( + const PointSet& points, + std::function pred, + unsigned max_points = 0); + +class SLASupportTree::Impl { + std::vector m_heads; + std::vector m_pillars; + std::vector m_junctions; + std::vector m_bridges; + std::vector m_compact_bridges; +public: + + template Head& add_head(Args&&... args) { + m_heads.emplace_back(std::forward(args)...); + m_heads.back().id = long(m_heads.size() - 1); + return m_heads.back(); + } + + template Pillar& add_pillar(long headid, Args&&... args) { + assert(headid >= 0 && headid < m_heads.size()); + Head& head = m_heads[headid]; + m_pillars.emplace_back(head, std::forward(args)...); + Pillar& pillar = m_pillars.back(); + pillar.id = long(m_pillars.size() - 1); + head.pillar_id = pillar.id; + pillar.start_junction_id = head.id; + pillar.starts_from_head = true; + return m_pillars.back(); + } + + const Head& pillar_head(long pillar_id) { + assert(pillar_id > 0 && pillar_id < m_pillars.size()); + Pillar& p = m_pillars[pillar_id]; + assert(p.starts_from_head && p.start_junction_id > 0 && + p.start_junction_id < m_heads.size() ); + return m_heads[p.start_junction_id]; + } + + const Pillar& head_pillar(long headid) { + assert(headid >= 0 && headid < m_heads.size()); + Head& h = m_heads[headid]; + assert(h.pillar_id > 0 && h.pillar_id < m_pillars.size()); + return m_pillars[h.pillar_id]; + } + + template const Junction& add_junction(Args&&... args) { + m_junctions.emplace_back(std::forward(args)...); + m_junctions.back().id = long(m_junctions.size() - 1); + return m_junctions.back(); + } + + template const Bridge& add_bridge(Args&&... args) { + m_bridges.emplace_back(std::forward(args)...); + m_bridges.back().id = long(m_bridges.size() - 1); + return m_bridges.back(); + } + + template + const CompactBridge& add_compact_bridge(Args&&...args) { + m_compact_bridges.emplace_back(std::forward(args)...); + m_compact_bridges.back().id = long(m_compact_bridges.size() - 1); + return m_compact_bridges.back(); + } + + const std::vector& heads() const { return m_heads; } + Head& head(size_t idx) { return m_heads[idx]; } + const std::vector& pillars() const { return m_pillars; } + const std::vector& bridges() const { return m_bridges; } + const std::vector& junctions() const { return m_junctions; } + const std::vector& compact_bridges() const { + return m_compact_bridges; + } +}; + +template +long cluster_centroid(const ClusterEl& clust, + std::function pointfn, + DistFn df) +{ + switch(clust.size()) { + case 0: /* empty cluster */ return -1; + case 1: /* only one element */ return 0; + case 2: /* if two elements, there is no center */ return 0; + default: ; + } + + // The function works by calculating for each point the average distance + // from all the other points in the cluster. We create a selector bitmask of + // the same size as the cluster. The bitmask will have two true bits and + // false bits for the rest of items and we will loop through all the + // permutations of the bitmask (combinations of two points). Get the + // distance for the two points and add the distance to the averages. + // The point with the smallest average than wins. + + std::vector sel(clust.size(), false); // create full zero bitmask + std::fill(sel.end() - 2, sel.end(), true); // insert the two ones + std::vector avgs(clust.size(), 0.0); // store the average distances + + do { + std::array idx; + for(size_t i = 0, j = 0; i < clust.size(); i++) if(sel[i]) idx[j++] = i; + + double d = df(pointfn(clust[idx[0]]), + pointfn(clust[idx[1]])); + + // add the distance to the sums for both associated points + for(auto i : idx) avgs[i] += d; + + // now continue with the next permutation of the bitmask with two 1s + } while(std::next_permutation(sel.begin(), sel.end())); + + // Divide by point size in the cluster to get the average (may be redundant) + for(auto& a : avgs) a /= clust.size(); + + // get the lowest average distance and return the index + auto minit = std::min_element(avgs.begin(), avgs.end()); + return long(minit - avgs.begin()); +} + +/** + * This function will calculate the convex hull of the input point set and + * return the indices of those points belonging to the chull in the right + * (counter clockwise) order. The input is also the set of indices and a + * functor to get the actual point form the index. + * + * I've adapted this algorithm from here: + * https://www.geeksforgeeks.org/convex-hull-set-1-jarviss-algorithm-or-wrapping/ + * and modified it so that it starts with the leftmost lower vertex. Also added + * support for floating point coordinates. + * + * This function is a modded version of the standard convex hull. If the points + * are all collinear with each other, it will return their indices in spatially + * subsequent order (the order they appear on the screen). + */ +ClusterEl pts_convex_hull(const ClusterEl& inpts, + std::function pfn) +{ + using Point = Vec2d; + using std::vector; + + static const double ERR = 1e-3; + + auto orientation = [](const Point& p, const Point& q, const Point& r) + { + double val = (q(Y) - p(Y)) * (r(X) - q(X)) - + (q(X) - p(X)) * (r(Y) - q(Y)); + + if (std::abs(val) < ERR) return 0; // collinear + return (val > ERR)? 1: 2; // clock or counterclockwise + }; + + size_t n = inpts.size(); + + if (n < 3) return inpts; + + // Initialize Result + ClusterEl hull; + vector points; points.reserve(n); + for(auto i : inpts) { + points.emplace_back(pfn(i)); + } + + // Check if the triplet of points is collinear. The standard convex hull + // algorithms are not capable of handling such input properly. + bool collinear = true; + for(auto one = points.begin(), two = std::next(one), three = std::next(two); + three != points.end() && collinear; + ++one, ++two, ++three) + { + // check if the points are collinear + if(orientation(*one, *two, *three) != 0) collinear = false; + } + + // Find the leftmost (bottom) point + int l = 0; + for (int i = 1; i < n; i++) { + if(std::abs(points[i](X) - points[l](X)) < ERR) { + if(points[i](Y) < points[l](Y)) l = i; + } + else if (points[i](X) < points[l](X)) l = i; + } + + if(collinear) { + // fill the output with the spatially ordered set of points. + + // find the direction + Vec2d dir = (points[l] - points[(l+1)%n]).normalized(); + hull = inpts; + auto& lp = points[l]; + std::sort(hull.begin(), hull.end(), + [&lp, points](unsigned i1, unsigned i2) { + // compare the distance from the leftmost point + return distance(lp, points[i1]) < distance(lp, points[i2]); + }); + + return hull; + } + + // TODO: this algorithm is O(m*n) and O(n^2) in the worst case so it needs + // to be replaced with a graham scan or something O(nlogn) + + // Start from leftmost point, keep moving counterclockwise + // until reach the start point again. This loop runs O(h) + // times where h is number of points in result or output. + int p = l; + do + { + // Add current point to result + hull.push_back(inpts[p]); + + // Search for a point 'q' such that orientation(p, x, + // q) is counterclockwise for all points 'x'. The idea + // is to keep track of last visited most counterclock- + // wise point in q. If any point 'i' is more counterclock- + // wise than q, then update q. + int q = (p+1)%n; + for (int i = 0; i < n; i++) + { + // If i is more counterclockwise than current q, then + // update q + if (orientation(points[p], points[i], points[q]) == 2) q = i; + } + + // Now q is the most counterclockwise with respect to p + // Set p as q for next iteration, so that q is added to + // result 'hull' + p = q; + + } while (p != l); // While we don't come to first point + + auto first = hull.front(); + hull.emplace_back(first); + + return hull; +} + +Vec3d dirv(const Vec3d& startp, const Vec3d& endp) { + return (endp - startp).normalized(); +} + +/// Generation of the supports, entry point function. This is called from the +/// SLASupportTree constructor and throws an SLASupportsStoppedException if it +/// gets canceled by the ctl object's stopcondition functor. +bool SLASupportTree::generate(const PointSet &points, + const EigenMesh3D& mesh, + const SupportConfig &cfg, + const Controller &ctl) +{ + PointSet filtered_points; // all valid support points + PointSet head_positions; // support points with pinhead + PointSet head_normals; // head normals + PointSet headless_positions; // headless support points + PointSet headless_normals; // headless support point normals + + using IndexSet = std::vector; + + // Distances from head positions to ground or mesh touch points + std::vector head_heights; + + // Indices of those who touch the ground + IndexSet ground_heads; + + // Indices of those who don't touch the ground + IndexSet noground_heads; + + ClusteredPoints ground_connectors; + + auto gnd_head_pt = [&ground_heads, &head_positions] (size_t idx) { + return Vec3d(head_positions.row(ground_heads[idx])); + }; + + using Result = SLASupportTree::Impl; + + Result& result = *m_impl; + + enum Steps { + BEGIN, + FILTER, + PINHEADS, + CLASSIFY, + ROUTING_GROUND, + ROUTING_NONGROUND, + HEADLESS, + DONE, + HALT, + ABORT, + NUM_STEPS + //... + }; + + auto filterfn = [] ( + const SupportConfig& cfg, + const PointSet& points, + const EigenMesh3D& mesh, + PointSet& filt_pts, + PointSet& head_norm, + PointSet& head_pos, + PointSet& headless_pos, + PointSet& headless_norm) + { + + /* ******************************************************** */ + /* Filtering step */ + /* ******************************************************** */ + + // Get the points that are too close to each other and keep only the + // first one + auto aliases = cluster(points, + [cfg](const SpatElement& p, + const SpatElement& se){ + return distance(p.first, se.first) < D_SP; + }, 2); + + filt_pts.resize(aliases.size(), 3); + int count = 0; + for(auto& a : aliases) { + // Here we keep only the front point of the cluster. TODO: centroid + filt_pts.row(count++) = points.row(a.front()); + } + + // calculate the normals to the triangles belonging to filtered points + auto nmls = sla::normals(filt_pts, mesh); + + head_norm.resize(count, 3); + head_pos.resize(count, 3); + headless_pos.resize(count, 3); + headless_norm.resize(count, 3); + + // Not all of the support points have to be a valid position for + // support creation. The angle may be inappropriate or there may + // not be enough space for the pinhead. Filtering is applied for + // these reasons. + + int pcount = 0, hlcount = 0; + for(int i = 0; i < count; i++) { + auto n = nmls.row(i); + + // for all normals we generate the spherical coordinates and + // saturate the polar angle to 45 degrees from the bottom then + // convert back to standard coordinates to get the new normal. + // Then we just create a quaternion from the two normals + // (Quaternion::FromTwoVectors) and apply the rotation to the + // arrow head. + + double z = n(2); + double r = 1.0; // for normalized vector + double polar = std::acos(z / r); + double azimuth = std::atan2(n(1), n(0)); + + if(polar >= PI / 2) { // skip if the tilt is not sane + + // We saturate the polar angle to 3pi/4 + polar = std::max(polar, 3*PI / 4); + + // Reassemble the now corrected normal + Vec3d nn(std::cos(azimuth) * std::sin(polar), + std::sin(azimuth) * std::sin(polar), + std::cos(polar)); + + // save the head (pinpoint) position + Vec3d hp = filt_pts.row(i); + + // the full width of the head + double w = cfg.head_width_mm + + cfg.head_back_radius_mm + + 2*cfg.head_front_radius_mm; + + // We should shoot a ray in the direction of the pinhead and + // see if there is enough space for it + double t = ray_mesh_intersect(hp + 0.1*nn, nn, mesh); + + if(t > 2*w || std::isinf(t)) { + // 2*w because of lower and upper pinhead + + head_pos.row(pcount) = hp; + + // save the verified and corrected normal + head_norm.row(pcount) = nn; + + ++pcount; + } else { + headless_norm.row(hlcount) = nn; + headless_pos.row(hlcount++) = hp; + } + } + } + + head_pos.conservativeResize(pcount, Eigen::NoChange); + head_norm.conservativeResize(pcount, Eigen::NoChange); + headless_pos.conservativeResize(hlcount, Eigen::NoChange); + headless_norm.conservativeResize(hlcount, Eigen::NoChange); + }; + + // Function to write the pinheads into the result + auto pinheadfn = [] ( + const SupportConfig& cfg, + PointSet& head_pos, + PointSet& nmls, + Result& result + ) + { + + /* ******************************************************** */ + /* Generating Pinheads */ + /* ******************************************************** */ + + for (int i = 0; i < head_pos.rows(); ++i) { + result.add_head( + cfg.head_back_radius_mm, + cfg.head_front_radius_mm, + cfg.head_width_mm, + nmls.row(i), // dir + head_pos.row(i) // displacement + ); + } + }; + + // &filtered_points, &head_positions, &result, &mesh, + // &gndidx, &gndheight, &nogndidx, cfg + auto classifyfn = [] ( + const SupportConfig& cfg, + const EigenMesh3D& mesh, + PointSet& head_pos, + IndexSet& gndidx, + IndexSet& nogndidx, + std::vector& gndheight, + ClusteredPoints& ground_clusters, + Result& result + ) { + + /* ******************************************************** */ + /* Classification */ + /* ******************************************************** */ + + // We should first get the heads that reach the ground directly + gndheight.reserve(head_pos.rows()); + gndidx.reserve(head_pos.rows()); + nogndidx.reserve(head_pos.rows()); + + for(unsigned i = 0; i < head_pos.rows(); i++) { + auto& head = result.heads()[i]; + + Vec3d dir(0, 0, -1); + Vec3d startpoint = head.junction_point(); + + double t = ray_mesh_intersect(startpoint, dir, mesh); + + gndheight.emplace_back(t); + + if(std::isinf(t)) gndidx.emplace_back(i); + else nogndidx.emplace_back(i); + } + + PointSet gnd(gndidx.size(), 3); + + for(size_t i = 0; i < gndidx.size(); i++) + gnd.row(i) = head_pos.row(gndidx[i]); + + // We want to search for clusters of points that are far enough from + // each other in the XY plane to not cross their pillar bases + // These clusters of support points will join in one pillar, possibly in + // their centroid support point. + auto d_base = 2*cfg.base_radius_mm; + ground_clusters = cluster(gnd, + [d_base, &cfg](const SpatElement& p, const SpatElement& s){ + return distance(Vec2d(p.first(X), p.first(Y)), + Vec2d(s.first(X), s.first(Y))) < d_base; + }, 3); // max 3 heads to connect to one centroid + }; + + // Helper function for interconnecting two pillars with zig-zag bridges + auto interconnect = [&cfg]( + const Pillar& pillar, + const Pillar& nextpillar, + const EigenMesh3D& emesh, + Result& result) + { + const Head& phead = result.pillar_head(pillar.id); + const Head& nextphead = result.pillar_head(nextpillar.id); + + double d = 2*pillar.r; + const Vec3d& pp = pillar.endpoint.cwiseProduct(Vec3d{1, 1, 0}); + + Vec3d sj = phead.junction_point(); + sj(Z) = std::min(sj(Z), nextphead.junction_point()(Z)); + Vec3d ej = nextpillar.endpoint; + double pillar_dist = distance(Vec2d{sj(X), sj(Y)}, + Vec2d{ej(X), ej(Y)}); + double zstep = pillar_dist * std::tan(-cfg.tilt); + ej(Z) = sj(Z) + zstep; + + double chkd = ray_mesh_intersect(sj, dirv(sj, ej), emesh); + double bridge_distance = pillar_dist / std::cos(-cfg.tilt); + + // If the pillars are so close that they touch each other, + // there is no need to bridge them together. + if(pillar_dist > 2*cfg.pillar_radius_mm && + bridge_distance < cfg.max_bridge_length_mm) + while(sj(Z) > pillar.endpoint(Z) && + ej(Z) > nextpillar.endpoint(Z)) + { + if(chkd >= bridge_distance) { + result.add_bridge(sj, ej, pillar.r); + + // double bridging: (crosses) + if(bridge_distance > 2*cfg.base_radius_mm) { + // If the columns are close together, no need to + // double bridge them + Vec3d bsj(ej(X), ej(Y), sj(Z)); + Vec3d bej(sj(X), sj(Y), ej(Z)); + + // need to check collision for the cross stick + double backchkd = ray_mesh_intersect(bsj, + dirv(bsj, bej), + emesh); + + if(backchkd >= bridge_distance) { + result.add_bridge(bsj, bej, pillar.r); + } + } + } + sj.swap(ej); + ej(Z) = sj(Z) + zstep; + chkd = ray_mesh_intersect(sj, dirv(sj, ej), emesh); + } + }; + + auto routing_ground_fn = [gnd_head_pt, interconnect]( + const SupportConfig& cfg, + const ClusteredPoints& gnd_clusters, + const IndexSet& gndidx, + const EigenMesh3D& emesh, + Result& result) + { + const double hbr = cfg.head_back_radius_mm; + const double pradius = cfg.pillar_radius_mm; + const double maxbridgelen = cfg.max_bridge_length_mm; + + ClusterEl cl_centroids; + cl_centroids.reserve(gnd_clusters.size()); + + SpatIndex pheadindex; // spatial index for the junctions + for(auto cl : gnd_clusters) { + // place all the centroid head positions into the index. We will + // query for alternative pillar positions. If a sidehead cannot + // connect to the cluster centroid, we have to search for another + // head with a full pillar. Also when there are two elements in the + // cluster, the centroid is arbitrary and the sidehead is allowed to + // connect to a nearby pillar to increase structural stability. + + // get the current cluster centroid + unsigned cid = cluster_centroid(cl, gnd_head_pt, + [](const Vec3d& p1, const Vec3d& p2) + { + return distance(Vec2d(p1(X), p1(Y)), Vec2d(p2(X), p2(Y))); + }); + + cl_centroids.push_back(cid); + + unsigned hid = gndidx[cl[cid]]; // Head index + Head& h = result.head(hid); + h.transform(); + Vec3d p = h.junction_point(); p(Z) = 0; + + pheadindex.insert(p, hid); + } + + // now we will go through the clusters ones again and connect the + // sidepoints with the cluster centroid (which is a ground pillar) + // or a nearby pillar if the centroid is unreachable. + for(auto cle : enumerate(gnd_clusters)) { + auto cl = cle.value; + auto cidx = cl_centroids[cle.index]; + cl_centroids[cle.index] = cl[cidx]; + + long index_to_heads = gndidx[cl[cidx]]; + auto& head = result.head(index_to_heads); + + Vec3d startpoint = head.junction_point(); + auto endpoint = startpoint; endpoint(Z) = 0; + + // Create the central pillar of the cluster with its base on the + // ground + result.add_pillar(index_to_heads, endpoint, pradius) + .add_base(cfg.base_height_mm, cfg.base_radius_mm); + + // Process side point in current cluster + cl.erase(cl.begin() + cidx); // delete the centroid before looping + + // TODO: dont consider the cluster centroid but calculate a central + // position where the pillar can be placed. this way the weight + // is distributed more effectively on the pillar. + + auto search_nearest = + [&cfg, &result, &emesh, maxbridgelen] + (SpatIndex& spindex, const Vec3d& jsh) + { + long nearest_id = -1; + while(nearest_id < 0 && !spindex.empty()) { + // loop until a suitable head is not found + // if there is a pillar closer than the cluster center + // (this may happen as the clustering is not perfect) + // than we will bridge to this closer pillar + + Vec3d qp(jsh(X), jsh(Y), 0); + auto ne = spindex.nearest(qp, 1).front(); + const Head& nearhead = result.heads()[ne.second]; + + Vec3d jh = nearhead.junction_point(); + Vec3d jp = jsh; + double dist2d = distance(qp, ne.first); + + // Bridge endpoint on the main pillar + Vec3d jn(jh(X), jh(Y), jp(Z) + dist2d*std::tan(-cfg.tilt)); + + if(jn(Z) > jh(Z)) { + // if the main head is below the point where the bridge + // would connect, than we must adjust the bridge + // endpoints + double hdiff = jn(Z) - jh(Z); + jp(Z) -= hdiff; + jn(Z) -= hdiff; + } + + double d = distance(jp, jn); + if(jn(Z) <= 0 || d > maxbridgelen) break; + + double chkd = ray_mesh_intersect(jp, dirv(jp, jn), emesh); + if(chkd >= d) nearest_id = ne.second; + + spindex.remove(ne); + } + return nearest_id; + }; + + for(auto c : cl) { + auto& sidehead = result.head(gndidx[c]); + sidehead.transform(); + + Vec3d jsh = sidehead.junction_point(); + Vec3d jp2d = {jsh(X), jsh(Y), 0}; + SpatIndex spindex = pheadindex; + long nearest_id = search_nearest(spindex, jsh); + + // at this point we either have our pillar index or we have + // to connect the sidehead to the ground + if(nearest_id < 0) { + // Could not find a pillar, create one + Vec3d jp = jsh; jp(Z) = 0; + result.add_pillar(sidehead.id, jp, pradius). + add_base(cfg.base_height_mm, cfg.base_radius_mm); + + // connects to ground, eligible for bridging + cl_centroids.emplace_back(sidehead.id); + } else { + // Creating the bridge to the nearest pillar + + const Head& nearhead = result.heads()[nearest_id]; + Vec3d jp = jsh; + Vec3d jh = nearhead.junction_point(); + + double d = distance(Vec2d{jp(X), jp(Y)}, + Vec2d{jh(X), jh(Y)}); + Vec3d jn(jh(X), jh(Y), jp(Z) + d*std::tan(-cfg.tilt)); + + if(jn(Z) > jh(Z)) { + double hdiff = jn(Z) - jh(Z); + jp(Z) -= hdiff; + jn(Z) -= hdiff; + + // pillar without base, this does not connect to ground. + result.add_pillar(sidehead.id, jp, pradius); + } + + if(jp(Z) < jsh(Z)) result.add_junction(jp, hbr); + if(jn(Z) >= jh(Z)) result.add_junction(jn, hbr); + double r_pillar = sidehead.request_pillar_radius(pradius); + result.add_bridge(jp, jn, r_pillar); + } + } + } + + // We will break down the pillar positions in 2D into concentric rings. + // Connecting the pillars belonging to the same ring will prevent + // bridges from crossing each other. After bridging the rings we can + // create bridges between the rings without the possibility of crossing + // bridges. Two pillars will be bridged with X shaped stick pairs. + // If they are really close to each other, than only one stick will be + // used in zig-zag mode. + + // Breaking down the points into rings will be done with a modified + // convex hull algorithm (see pts_convex_hull()), that works for + // collinear points as well. If the points are on the same surface, + // they can be part of an imaginary line segment for which the convex + // hull is not defined. I this case it is enough to sort the points + // spatially and create the bridge stick from the one endpoint to + // another. + + ClusterEl rem = cl_centroids; + ClusterEl ring; + + while(!rem.empty()) { // loop until all the points belong to some ring + std::sort(rem.begin(), rem.end()); + + auto newring = pts_convex_hull(rem, + [gnd_head_pt](unsigned i) { + auto& p = gnd_head_pt(i); + return Vec2d(p(X), p(Y)); // project to 2D in along Z axis + }); + + if(!ring.empty()) { + // inner ring is now in 'newring' and outer ring is in 'ring' + SpatIndex innerring; + for(unsigned i : newring) { + const Pillar& pill = result.head_pillar(gndidx[i]); + innerring.insert(pill.endpoint, pill.id); + } + + // For all pillars in the outer ring find the closest in the + // inner ring and connect them. This will create the spider web + // fashioned connections between pillars + for(unsigned i : ring) { + const Pillar& outerpill = result.head_pillar(gndidx[i]); + auto res = innerring.nearest(outerpill.endpoint, 1); + if(res.empty()) continue; + + auto ne = res.front(); + const Pillar& innerpill = result.pillars()[ne.second]; + interconnect(outerpill, innerpill, emesh, result); + } + } + + // no need for newring anymore in the current iteration + ring.swap(newring); + + /*std::cout << "ring: \n"; + for(auto ri : ring) { + std::cout << ri << " " << " X = " << gnd_head_pt(ri)(X) + << " Y = " << gnd_head_pt(ri)(Y) << std::endl; + } + std::cout << std::endl;*/ + + // now the ring has to be connected with bridge sticks + for(auto it = ring.begin(), next = std::next(it); + next != ring.end(); + ++it, ++next) + { + const Pillar& pillar = result.head_pillar(gndidx[*it]); + const Pillar& nextpillar = result.head_pillar(gndidx[*next]); + interconnect(pillar, nextpillar, emesh, result); + } + + auto sring = ring; ClusterEl tmp; + std::sort(sring.begin(), sring.end()); + std::set_difference(rem.begin(), rem.end(), + sring.begin(), sring.end(), + std::back_inserter(tmp)); + rem.swap(tmp); + } + }; + + auto routing_nongnd_fn = []( + const SupportConfig& cfg, + const std::vector& gndheight, + const IndexSet& nogndidx, + Result& result) + { + // TODO: connect these to the ground pillars if possible + for(auto idx : nogndidx) { + auto& head = result.head(idx); + head.transform(); + + double gh = gndheight[idx]; + Vec3d headend = head.junction_point(); + + Head base_head(cfg.head_back_radius_mm, + cfg.head_front_radius_mm, + cfg.head_width_mm, + {0.0, 0.0, 1.0}, + {headend(X), headend(Y), headend(Z) - gh}); + + base_head.transform(); + + double hl = head.fullwidth() - head.r_back_mm; + + result.add_pillar(idx, + Vec3d{headend(X), headend(Y), headend(Z) - gh + hl}, + cfg.pillar_radius_mm + ).base = base_head.mesh; + } + }; + + auto process_headless = []( + const SupportConfig& cfg, + const PointSet& headless_pts, + const PointSet& headless_norm, + const EigenMesh3D& emesh, + Result& result) + { + // For now we will just generate smaller headless sticks with a sharp + // ending point that connects to the mesh surface. + + const double R = 0.5*cfg.pillar_radius_mm; + const double HWIDTH_MM = R/3; + + // We will sink the pins into the model surface for a distance of 1/3 of + // HWIDTH_MM + for(int i = 0; i < headless_pts.rows(); i++) { + Vec3d sp = headless_pts.row(i); + + Vec3d n = headless_norm.row(i); + sp = sp - n * HWIDTH_MM; + + Vec3d dir = {0, 0, -1}; + Vec3d sj = sp + R * n; + double dist = ray_mesh_intersect(sj, dir, emesh); + + if(std::isinf(dist) || std::isnan(dist)) continue; + + Vec3d ej = sj + (dist + HWIDTH_MM)* dir; + result.add_compact_bridge(sp, ej, n, R); + } + }; + + using std::ref; + using std::cref; + using std::bind; + + // Here we can easily track what goes in and what comes out of each step: + // (see the cref-s as inputs and ref-s as outputs) + std::array, NUM_STEPS> program = { + [] () { + // Begin + // clear up the shared data + }, + + // Filtering unnecessary support points + bind(filterfn, cref(cfg), cref(points), cref(mesh), + ref(filtered_points), ref(head_normals), + ref(head_positions), ref(headless_positions), ref(headless_normals)), + + // Pinhead generation + bind(pinheadfn, cref(cfg), + ref(head_positions), ref(head_normals), ref(result)), + + // Classification of support points + bind(classifyfn, cref(cfg), cref(mesh), + ref(head_positions), ref(ground_heads), ref(noground_heads), + ref(head_heights), ref(ground_connectors), ref(result)), + + // Routing ground connecting clusters + bind(routing_ground_fn, + cref(cfg), cref(ground_connectors), cref(ground_heads), cref(mesh), + ref(result)), + + // routing non ground connecting support points + bind(routing_nongnd_fn, cref(cfg), cref(head_heights), cref(noground_heads), + ref(result)), + + bind(process_headless, + cref(cfg), cref(headless_positions), + cref(headless_normals), cref(mesh), + ref(result)), + [] () { + // Done + }, + [] () { + // Halt + }, + [] () { + // Abort + } + }; + + Steps pc = BEGIN, pc_prev = BEGIN; + + auto progress = [&ctl, &pc, &pc_prev] () { + static const std::array stepstr { + "", + "Filtering", + "Generate pinheads", + "Classification", + "Routing to ground", + "Routing supports to model surface", + "Processing small holes", + "Done", + "Halt", + "Abort" + }; + + static const std::array stepstate { + 0, + 10, + 30, + 50, + 60, + 70, + 80, + 100, + 0, + 0 + }; + + if(ctl.stopcondition()) pc = ABORT; + + switch(pc) { + case BEGIN: pc = FILTER; break; + case FILTER: pc = PINHEADS; break; + case PINHEADS: pc = CLASSIFY; break; + case CLASSIFY: pc = ROUTING_GROUND; break; + case ROUTING_GROUND: pc = ROUTING_NONGROUND; break; + case ROUTING_NONGROUND: pc = HEADLESS; break; + case HEADLESS: pc = DONE; break; + case HALT: pc = pc_prev; break; + case DONE: + case ABORT: break; + } + ctl.statuscb(stepstate[pc], stepstr[pc]); + }; + + // Just here we run the computation... + while(pc < DONE || pc == HALT) { + progress(); + program[pc](); + } + + if(pc == ABORT) throw SLASupportsStoppedException(); + + return pc == ABORT; +} + +void SLASupportTree::merged_mesh(TriangleMesh &outmesh) const +{ + const SLASupportTree::Impl& stree = get(); + + for(auto& head : stree.heads()) { + outmesh.merge(mesh(head.mesh)); + } + + for(auto& stick : stree.pillars()) { + outmesh.merge(mesh(stick.mesh)); + outmesh.merge(mesh(stick.base)); + } + + for(auto& j : stree.junctions()) { + outmesh.merge(mesh(j.mesh)); + } + + for(auto& cb : stree.compact_bridges()) { + outmesh.merge(mesh(cb.mesh)); + } + + for(auto& bs : stree.bridges()) { + outmesh.merge(mesh(bs.mesh)); + } +} + +SlicedSupports SLASupportTree::slice() const +{ + return {}; +} + +SLASupportTree::SLASupportTree(const Model& model, + const SupportConfig& cfg, + const Controller& ctl): m_impl(new Impl()) +{ + generate(support_points(model), to_eigenmesh(model), cfg, ctl); +} + +SLASupportTree::SLASupportTree(const PointSet &points, + const EigenMesh3D& emesh, + const SupportConfig &cfg, + const Controller &ctl): m_impl(new Impl()) +{ + generate(points, emesh, cfg, ctl); +} + +SLASupportTree::SLASupportTree(const SLASupportTree &c): + m_impl( new Impl(*c.m_impl)) {} + +SLASupportTree &SLASupportTree::operator=(const SLASupportTree &c) +{ + m_impl = make_unique(*c.m_impl); + return *this; +} + +SLASupportTree::~SLASupportTree() {} + +void add_sla_supports(Model &model, + const SupportConfig &cfg, + const Controller &ctl) +{ + Benchmark bench; + + bench.start(); + SLASupportTree _stree(model, cfg, ctl); + bench.stop(); + + std::cout << "Support tree creation time: " << bench.getElapsedSec() + << " seconds" << std::endl; + + bench.start(); + SLASupportTree::Impl& stree = _stree.get(); + ModelObject* o = model.add_object(); + o->add_instance(); + + TriangleMesh streemsh; + _stree.merged_mesh(streemsh); + o->add_volume(streemsh); + + bench.stop(); + std::cout << "support tree added to model in: " << bench.getElapsedSec() + << " seconds" << std::endl; + + // TODO this would roughly be the code for the base pool + ExPolygons plate; + auto modelmesh = model.mesh(); + TriangleMesh poolmesh; + sla::PoolConfig poolcfg; + poolcfg.min_wall_height_mm = 1; + poolcfg.edge_radius_mm = 0.1; + poolcfg.min_wall_thickness_mm = 0.8; + + bench.start(); + sla::base_plate(modelmesh, plate); + bench.stop(); + + std::cout << "Base plate calculation time: " << bench.getElapsedSec() + << " seconds." << std::endl; + + bench.start(); + sla::create_base_pool(plate, poolmesh, poolcfg); + bench.stop(); + + std::cout << "Pool generation completed in " << bench.getElapsedSec() + << " second." << std::endl; + + bench.start(); + poolmesh.translate(0, 0, poolcfg.min_wall_height_mm / 2); + o->add_volume(poolmesh); + bench.stop(); + + // TODO: will cause incorrect placement of the model; +// o->translate({0, 0, poolcfg.min_wall_height_mm / 2}); + + std::cout << "Added pool to model in " << bench.getElapsedSec() + << " seconds." << std::endl; + +} + +} +} diff --git a/src/libslic3r/SLA/SLASupportTree.hpp b/src/libslic3r/SLA/SLASupportTree.hpp new file mode 100644 index 000000000..d6e7afca5 --- /dev/null +++ b/src/libslic3r/SLA/SLASupportTree.hpp @@ -0,0 +1,132 @@ +#ifndef SLASUPPORTTREE_HPP +#define SLASUPPORTTREE_HPP + +#include +#include +#include +#include +#include + +namespace Slic3r { + +// Needed types from Point.hpp +typedef int32_t coord_t; +typedef Eigen::Matrix Vec3d; +typedef Eigen::Matrix Vec3f; +typedef Eigen::Matrix Vec3crd; +typedef std::vector Pointf3s; +typedef std::vector Points3; + +class TriangleMesh; +class Model; +class ModelInstance; +class ExPolygon; + +using SliceLayer = std::vector; +using SlicedSupports = std::vector; + +namespace sla { + +struct SupportConfig { + // Radius in mm of the pointing side of the head. + double head_front_radius_mm = 0.2; + + // Radius of the back side of the 3d arrow. + double head_back_radius_mm = 0.5; + + // Width in mm from the back sphere center to the front sphere center. + double head_width_mm = 1.0; + + // Radius in mm of the support pillars. + // Warning: this value will be at most 65% of head_back_radius_mm + double pillar_radius_mm = 0.8; + + // Radius in mm of the pillar base. + double base_radius_mm = 2.0; + + // The height of the pillar base cone in mm. + double base_height_mm = 1.0; + + // The default angle for connecting support sticks and junctions. + double tilt = M_PI/4; + + // The max length of a bridge in mm + double max_bridge_length_mm = 15.0; +}; + +/// A Control structure for the support calculation. Consists of the status +/// indicator callback and the stop condition predicate. +struct Controller { + std::function statuscb = + [](unsigned, const std::string&){}; + + std::function stopcondition = [](){ return false; }; +}; + +/* ************************************************************************** */ +/* TODO: May not be needed: */ +/* ************************************************************************** */ + +void create_head(TriangleMesh&, double r1_mm, double r2_mm, double width_mm); + +/// Add support volumes to the model directly +void add_sla_supports(Model& model, const SupportConfig& cfg = {}, + const Controller& ctl = {}); + +/* ************************************************************************** */ + +using PointSet = Eigen::MatrixXd; +struct EigenMesh3D; + +/// Just a wrapper to the runtime error to be recognizable in try blocks +class SLASupportsStoppedException: public std::runtime_error { +public: + using std::runtime_error::runtime_error; + SLASupportsStoppedException(): std::runtime_error("") {} +}; + +/// The class containing mesh data for the generated supports. +class SLASupportTree { + class Impl; + std::unique_ptr m_impl; + Impl& get() { return *m_impl; } + const Impl& get() const { return *m_impl; } + + friend void add_sla_supports(Model&, + const SupportConfig&, + const Controller&); + + /// Generate the 3D supports for a model intended for SLA print. + bool generate(const PointSet& pts, + const EigenMesh3D& mesh, + const SupportConfig& cfg = {}, + const Controller& ctl = {}); +public: + + // Constructors will throw if the stop condition becomes true. + SLASupportTree(const Model& model, + const SupportConfig& cfg = {}, + const Controller& ctl = {}); + + SLASupportTree(const PointSet& pts, + const EigenMesh3D& em, + const SupportConfig& cfg = {}, + const Controller& ctl = {}); + + SLASupportTree(const SLASupportTree&); + SLASupportTree& operator=(const SLASupportTree&); + + ~SLASupportTree(); + + /// Get the whole mesh united into the output TriangleMesh + void merged_mesh(TriangleMesh& outmesh) const; + + /// Get the sliced 2d layers of the support geometry. + SlicedSupports slice() const; +}; + +} + +} + +#endif // SLASUPPORTTREE_HPP diff --git a/src/libslic3r/SLA/SLASupportTreeIGL.cpp b/src/libslic3r/SLA/SLASupportTreeIGL.cpp new file mode 100644 index 000000000..88782ba06 --- /dev/null +++ b/src/libslic3r/SLA/SLASupportTreeIGL.cpp @@ -0,0 +1,237 @@ +#include "SLA/SLASupportTree.hpp" +#include "SLA/SLABoilerPlate.hpp" +#include "SLA/SLASpatIndex.hpp" + +// HEAVY headers... takes eternity to compile + +// for concave hull merging decisions +#include "SLABoostAdapter.hpp" +#include "boost/geometry/index/rtree.hpp" + +#include +#include +#include "SLASpatIndex.hpp" +#include "ClipperUtils.hpp" + +namespace Slic3r { +namespace sla { + +class SpatIndex::Impl { +public: + using BoostIndex = boost::geometry::index::rtree< SpatElement, + boost::geometry::index::rstar<16, 4> /* ? */ >; + + BoostIndex m_store; +}; + +SpatIndex::SpatIndex(): m_impl(new Impl()) {} +SpatIndex::~SpatIndex() {} + +SpatIndex::SpatIndex(const SpatIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {} +SpatIndex::SpatIndex(SpatIndex&& cpy): m_impl(std::move(cpy.m_impl)) {} + +SpatIndex& SpatIndex::operator=(const SpatIndex &cpy) +{ + m_impl.reset(new Impl(*cpy.m_impl)); + return *this; +} + +SpatIndex& SpatIndex::operator=(SpatIndex &&cpy) +{ + m_impl.swap(cpy.m_impl); + return *this; +} + +void SpatIndex::insert(const SpatElement &el) +{ + m_impl->m_store.insert(el); +} + +bool SpatIndex::remove(const SpatElement& el) +{ + return m_impl->m_store.remove(el); +} + +std::vector +SpatIndex::query(std::function fn) +{ + namespace bgi = boost::geometry::index; + + std::vector ret; + m_impl->m_store.query(bgi::satisfies(fn), std::back_inserter(ret)); + return ret; +} + +std::vector SpatIndex::nearest(const Vec3d &el, unsigned k = 1) +{ + namespace bgi = boost::geometry::index; + std::vector ret; ret.reserve(k); + m_impl->m_store.query(bgi::nearest(el, k), std::back_inserter(ret)); + return ret; +} + +size_t SpatIndex::size() const +{ + return m_impl->m_store.size(); +} + +PointSet normals(const PointSet& points, const EigenMesh3D& mesh) { + Eigen::VectorXd dists; + Eigen::VectorXi I; + PointSet C; + igl::point_mesh_squared_distance( points, mesh.V, mesh.F, dists, I, C); + + PointSet ret(I.rows(), 3); + for(int i = 0; i < I.rows(); i++) { + auto idx = I(i); + auto trindex = mesh.F.row(idx); + + auto& p1 = mesh.V.row(trindex(0)); + auto& p2 = mesh.V.row(trindex(1)); + auto& p3 = mesh.V.row(trindex(2)); + + Eigen::Vector3d U = p2 - p1; + Eigen::Vector3d V = p3 - p1; + ret.row(i) = U.cross(V).normalized(); + } + + return ret; +} + +double ray_mesh_intersect(const Vec3d& s, + const Vec3d& dir, + const EigenMesh3D& m) +{ + igl::Hit hit; + hit.t = std::numeric_limits::infinity(); + igl::ray_mesh_intersect(s, dir, m.V, m.F, hit); + return hit.t; +} + +// Clustering a set of points by the given criteria +ClusteredPoints cluster( + const sla::PointSet& points, + std::function pred, + unsigned max_points = 0) +{ + + namespace bgi = boost::geometry::index; + using Index3D = bgi::rtree< SpatElement, bgi::rstar<16, 4> /* ? */ >; + + // A spatial index for querying the nearest points + Index3D sindex; + + // Build the index + for(unsigned idx = 0; idx < points.rows(); idx++) + sindex.insert( std::make_pair(points.row(idx), idx)); + + using Elems = std::vector; + + // Recursive function for visiting all the points in a given distance to + // each other + std::function group = + [&sindex, &group, pred, max_points](Elems& pts, Elems& cluster) + { + for(auto& p : pts) { + std::vector tmp; + + sindex.query( + bgi::satisfies([p, pred](const SpatElement& se) { + return pred(p, se); + }), + std::back_inserter(tmp) + ); + + auto cmp = [](const SpatElement& e1, const SpatElement& e2){ + return e1.second < e2.second; + }; + + std::sort(tmp.begin(), tmp.end(), cmp); + + Elems newpts; + std::set_difference(tmp.begin(), tmp.end(), + cluster.begin(), cluster.end(), + std::back_inserter(newpts), cmp); + + int c = max_points && newpts.size() + cluster.size() > max_points? + int(max_points - cluster.size()) : int(newpts.size()); + + cluster.insert(cluster.end(), newpts.begin(), newpts.begin() + c); + std::sort(cluster.begin(), cluster.end(), cmp); + + if(!newpts.empty() && (!max_points || cluster.size() < max_points)) + group(newpts, cluster); + } + }; + + std::vector clusters; + for(auto it = sindex.begin(); it != sindex.end();) { + Elems cluster = {}; + Elems pts = {*it}; + group(pts, cluster); + + for(auto& c : cluster) sindex.remove(c); + it = sindex.begin(); + + clusters.emplace_back(cluster); + } + + ClusteredPoints result; + for(auto& cluster : clusters) { + result.emplace_back(); + for(auto c : cluster) result.back().emplace_back(c.second); + } + + return result; +} + +using Segments = std::vector>; + +Segments model_boundary(const EigenMesh3D& emesh, double offs) +{ + Segments ret; + Polygons pp; + pp.reserve(emesh.F.rows()); + + for (int i = 0; i < emesh.F.rows(); i++) { + auto trindex = emesh.F.row(i); + auto& p1 = emesh.V.row(trindex(0)); + auto& p2 = emesh.V.row(trindex(1)); + auto& p3 = emesh.V.row(trindex(2)); + + Polygon p; + p.points.resize(3); + p.points[0] = Point::new_scale(p1(X), p1(Y)); + p.points[1] = Point::new_scale(p2(X), p2(Y)); + p.points[2] = Point::new_scale(p3(X), p3(Y)); + p.make_counter_clockwise(); + pp.emplace_back(p); + } + + ExPolygons merged = union_ex(offset(pp, float(scale_(offs))), true); + + for(auto& expoly : merged) { + auto lines = expoly.lines(); + for(Line& l : lines) { + Vec2d a(l.a(X) * SCALING_FACTOR, l.a(Y) * SCALING_FACTOR); + Vec2d b(l.b(X) * SCALING_FACTOR, l.b(Y) * SCALING_FACTOR); + ret.emplace_back(std::make_pair(a, b)); + } + } + + return ret; +} + +//struct SegmentIndex { + +//}; + +//using SegmentIndexEl = std::pair; + +//SegmentIndexEl + + + + +} +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 000000000..11bdc4b3d --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,3 @@ +# TODO Add individual tests as executables in separate directories + +# add_subirectory() \ No newline at end of file