diff --git a/src/libslic3r/SLA/IndexedMesh.cpp b/src/libslic3r/AABBMesh.cpp similarity index 51% rename from src/libslic3r/SLA/IndexedMesh.cpp rename to src/libslic3r/AABBMesh.cpp index 07c4203ab..835df2ebe 100644 --- a/src/libslic3r/SLA/IndexedMesh.cpp +++ b/src/libslic3r/AABBMesh.cpp @@ -1,5 +1,5 @@ -#include "IndexedMesh.hpp" -#include "Concurrency.hpp" +#include "AABBMesh.hpp" +#include #include #include @@ -12,9 +12,7 @@ namespace Slic3r { -namespace sla { - -class IndexedMesh::AABBImpl { +class AABBMesh::AABBImpl { private: AABBTreeIndirect::Tree3f m_tree; double m_triangle_ray_epsilon; @@ -68,7 +66,7 @@ public: } }; -template void IndexedMesh::init(const M &mesh, bool calculate_epsilon) +template void AABBMesh::init(const M &mesh, bool calculate_epsilon) { BoundingBoxf3 bb = bounding_box(mesh); m_ground_level += bb.min(Z); @@ -77,73 +75,86 @@ template void IndexedMesh::init(const M &mesh, bool calculate_epsilon) m_aabb->init(*m_tm, calculate_epsilon); } -IndexedMesh::IndexedMesh(const indexed_triangle_set& tmesh, bool calculate_epsilon) - : m_aabb(new AABBImpl()), m_tm(&tmesh) +AABBMesh::AABBMesh(const indexed_triangle_set &tmesh, bool calculate_epsilon) + : m_tm(&tmesh) + , m_aabb(new AABBImpl()) + , m_vfidx{tmesh} + , m_fnidx{its_face_neighbors(tmesh)} { init(tmesh, calculate_epsilon); } -IndexedMesh::IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon) - : m_aabb(new AABBImpl()), m_tm(&mesh.its) +AABBMesh::AABBMesh(const TriangleMesh &mesh, bool calculate_epsilon) + : m_tm(&mesh.its) + , m_aabb(new AABBImpl()) + , m_vfidx{mesh.its} + , m_fnidx{its_face_neighbors(mesh.its)} { init(mesh, calculate_epsilon); } -IndexedMesh::~IndexedMesh() {} +AABBMesh::~AABBMesh() {} -IndexedMesh::IndexedMesh(const IndexedMesh &other): - m_tm(other.m_tm), m_ground_level(other.m_ground_level), - m_aabb( new AABBImpl(*other.m_aabb) ) {} +AABBMesh::AABBMesh(const AABBMesh &other) + : m_tm(other.m_tm) + , m_ground_level(other.m_ground_level) + , m_aabb(new AABBImpl(*other.m_aabb)) + , m_vfidx{other.m_vfidx} + , m_fnidx{other.m_fnidx} +{} - -IndexedMesh &IndexedMesh::operator=(const IndexedMesh &other) +AABBMesh &AABBMesh::operator=(const AABBMesh &other) { m_tm = other.m_tm; m_ground_level = other.m_ground_level; - m_aabb.reset(new AABBImpl(*other.m_aabb)); return *this; + m_aabb.reset(new AABBImpl(*other.m_aabb)); + m_vfidx = other.m_vfidx; + m_fnidx = other.m_fnidx; + + return *this; } -IndexedMesh &IndexedMesh::operator=(IndexedMesh &&other) = default; +AABBMesh &AABBMesh::operator=(AABBMesh &&other) = default; -IndexedMesh::IndexedMesh(IndexedMesh &&other) = default; +AABBMesh::AABBMesh(AABBMesh &&other) = default; -const std::vector& IndexedMesh::vertices() const +const std::vector& AABBMesh::vertices() const { return m_tm->vertices; } -const std::vector& IndexedMesh::indices() const +const std::vector& AABBMesh::indices() const { return m_tm->indices; } -const Vec3f& IndexedMesh::vertices(size_t idx) const +const Vec3f& AABBMesh::vertices(size_t idx) const { return m_tm->vertices[idx]; } -const Vec3i& IndexedMesh::indices(size_t idx) const +const Vec3i& AABBMesh::indices(size_t idx) const { return m_tm->indices[idx]; } -Vec3d IndexedMesh::normal_by_face_id(int face_id) const { +Vec3d AABBMesh::normal_by_face_id(int face_id) const { return its_unnormalized_normal(*m_tm, face_id).cast().normalized(); } -IndexedMesh::hit_result -IndexedMesh::query_ray_hit(const Vec3d &s, const Vec3d &dir) const +AABBMesh::hit_result +AABBMesh::query_ray_hit(const Vec3d &s, const Vec3d &dir) const { assert(is_approx(dir.norm(), 1.)); igl::Hit hit{-1, -1, 0.f, 0.f, 0.f}; @@ -171,10 +182,10 @@ IndexedMesh::query_ray_hit(const Vec3d &s, const Vec3d &dir) const return ret; } -std::vector -IndexedMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const +std::vector +AABBMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const { - std::vector outs; + std::vector outs; std::vector hits; m_aabb->intersect_ray(*m_tm, s, dir, hits); @@ -192,7 +203,7 @@ IndexedMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const // Convert the igl::Hit into hit_result outs.reserve(hits.size()); for (const igl::Hit& hit : hits) { - outs.emplace_back(IndexedMesh::hit_result(*this)); + outs.emplace_back(AABBMesh::hit_result(*this)); outs.back().m_t = double(hit.t); outs.back().m_dir = dir; outs.back().m_source = s; @@ -207,8 +218,8 @@ IndexedMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const #ifdef SLIC3R_HOLE_RAYCASTER -IndexedMesh::hit_result IndexedMesh::filter_hits( - const std::vector& object_hits) const +AABBMesh::hit_result IndexedMesh::filter_hits( + const std::vector& object_hits) const { assert(! m_holes.empty()); hit_result out(*this); @@ -304,7 +315,7 @@ IndexedMesh::hit_result IndexedMesh::filter_hits( #endif -double IndexedMesh::squared_distance(const Vec3d &p, int& i, Vec3d& c) const { +double AABBMesh::squared_distance(const Vec3d &p, int& i, Vec3d& c) const { double sqdst = 0; Eigen::Matrix pp = p; Eigen::Matrix cc; @@ -313,143 +324,4 @@ double IndexedMesh::squared_distance(const Vec3d &p, int& i, Vec3d& c) const { return sqdst; } - -static bool point_on_edge(const Vec3d& p, const Vec3d& e1, const Vec3d& e2, - double eps = 0.05) -{ - using Line3D = Eigen::ParametrizedLine; - - auto line = Line3D::Through(e1, e2); - double d = line.distance(p); - return std::abs(d) < eps; -} - -PointSet normals(const PointSet& points, - const IndexedMesh& mesh, - double eps, - std::function thr, // throw on cancel - const std::vector& pt_indices) -{ - if (points.rows() == 0 || mesh.vertices().empty() || mesh.indices().empty()) - return {}; - - std::vector range = pt_indices; - if (range.empty()) { - range.resize(size_t(points.rows()), 0); - std::iota(range.begin(), range.end(), 0); - } - - PointSet ret(range.size(), 3); - - // for (size_t ridx = 0; ridx < range.size(); ++ridx) - ccr::for_each(size_t(0), range.size(), - [&ret, &mesh, &points, thr, eps, &range](size_t ridx) { - thr(); - unsigned el = range[ridx]; - auto eidx = Eigen::Index(el); - int faceid = 0; - Vec3d p; - - mesh.squared_distance(points.row(eidx), faceid, p); - - auto trindex = mesh.indices(faceid); - - const Vec3d &p1 = mesh.vertices(trindex(0)).cast(); - const Vec3d &p2 = mesh.vertices(trindex(1)).cast(); - const Vec3d &p3 = mesh.vertices(trindex(2)).cast(); - - // We should check if the point lies on an edge of the hosting - // triangle. If it does then all the other triangles using the - // same two points have to be searched and the final normal should - // be some kind of aggregation of the participating triangle - // normals. We should also consider the cases where the support - // point lies right on a vertex of its triangle. The procedure is - // the same, get the neighbor triangles and calculate an average - // normal. - - // mark the vertex indices of the edge. ia and ib marks and edge - // ic will mark a single vertex. - int ia = -1, ib = -1, ic = -1; - - if (std::abs((p - p1).norm()) < eps) { - ic = trindex(0); - } else if (std::abs((p - p2).norm()) < eps) { - ic = trindex(1); - } else if (std::abs((p - p3).norm()) < eps) { - ic = trindex(2); - } else if (point_on_edge(p, p1, p2, eps)) { - ia = trindex(0); - ib = trindex(1); - } else if (point_on_edge(p, p2, p3, eps)) { - ia = trindex(1); - ib = trindex(2); - } else if (point_on_edge(p, p1, p3, eps)) { - ia = trindex(0); - ib = trindex(2); - } - - // vector for the neigboring triangles including the detected one. - std::vector neigh; - if (ic >= 0) { // The point is right on a vertex of the triangle - for (size_t n = 0; n < mesh.indices().size(); ++n) { - thr(); - Vec3i ni = mesh.indices(n); - if ((ni(X) == ic || ni(Y) == ic || ni(Z) == ic)) - neigh.emplace_back(n); - } - } else if (ia >= 0 && ib >= 0) { // the point is on and edge - // now get all the neigboring triangles - for (size_t n = 0; n < mesh.indices().size(); ++n) { - thr(); - Vec3i ni = mesh.indices(n); - if ((ni(X) == ia || ni(Y) == ia || ni(Z) == ia) && - (ni(X) == ib || ni(Y) == ib || ni(Z) == ib)) - neigh.emplace_back(n); - } - } - - // Calculate the normals for the neighboring triangles - std::vector neighnorms; - neighnorms.reserve(neigh.size()); - for (size_t &tri_id : neigh) - neighnorms.emplace_back(mesh.normal_by_face_id(tri_id)); - - // Throw out duplicates. They would cause trouble with summing. We - // will use std::unique which works on sorted ranges. We will sort - // by the coefficient-wise sum of the normals. It should force the - // same elements to be consecutive. - std::sort(neighnorms.begin(), neighnorms.end(), - [](const Vec3d &v1, const Vec3d &v2) { - return v1.sum() < v2.sum(); - }); - - auto lend = std::unique(neighnorms.begin(), neighnorms.end(), - [](const Vec3d &n1, const Vec3d &n2) { - // Compare normals for equivalence. - // This is controvers stuff. - auto deq = [](double a, double b) { - return std::abs(a - b) < 1e-3; - }; - return deq(n1(X), n2(X)) && - deq(n1(Y), n2(Y)) && - deq(n1(Z), n2(Z)); - }); - - if (!neighnorms.empty()) { // there were neighbors to count with - // sum up the normals and then normalize the result again. - // This unification seems to be enough. - Vec3d sumnorm(0, 0, 0); - sumnorm = std::accumulate(neighnorms.begin(), lend, sumnorm); - sumnorm.normalize(); - ret.row(long(ridx)) = sumnorm; - } else { // point lies safely within its triangle - Eigen::Vector3d U = p2 - p1; - Eigen::Vector3d V = p3 - p1; - ret.row(long(ridx)) = U.cross(V).normalized(); - } - }); - - return ret; -} - -}} // namespace Slic3r::sla +} // namespace Slic3r diff --git a/src/libslic3r/SLA/IndexedMesh.hpp b/src/libslic3r/AABBMesh.hpp similarity index 68% rename from src/libslic3r/SLA/IndexedMesh.hpp rename to src/libslic3r/AABBMesh.hpp index 9348a97c9..bcc1525de 100644 --- a/src/libslic3r/SLA/IndexedMesh.hpp +++ b/src/libslic3r/AABBMesh.hpp @@ -1,5 +1,5 @@ -#ifndef SLA_INDEXEDMESH_H -#define SLA_INDEXEDMESH_H +#ifndef PRUSASLICER_AABBMESH_H +#define PRUSASLICER_AABBMESH_H #include #include @@ -21,25 +21,22 @@ namespace Slic3r { class TriangleMesh; -namespace sla { - -using PointSet = Eigen::MatrixXd; - -/// An index-triangle structure for libIGL functions. Also serves as an -/// alternative (raw) input format for the SLASupportTree. -// Implemented in libslic3r/SLA/Common.cpp -class IndexedMesh { +// An index-triangle structure coupled with an AABB index to support ray +// casting and other higher level operations. +class AABBMesh { class AABBImpl; const indexed_triangle_set* m_tm; - double m_ground_level = 0, m_gnd_offset = 0; + double m_ground_level = 0/*, m_gnd_offset = 0*/; std::unique_ptr m_aabb; + VertexFaceIndex m_vfidx; // vertex-face index + std::vector m_fnidx; // face-neighbor index #ifdef SLIC3R_HOLE_RAYCASTER // This holds a copy of holes in the mesh. Initialized externally // by load_mesh setter. - std::vector m_holes; + std::vector m_holes; #endif template void init(const M &mesh, bool calculate_epsilon); @@ -48,20 +45,20 @@ public: // calculate_epsilon ... calculate epsilon for triangle-ray intersection from an average triangle edge length. // If set to false, a default epsilon is used, which works for "reasonable" meshes. - explicit IndexedMesh(const indexed_triangle_set &tmesh, bool calculate_epsilon = false); - explicit IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon = false); + explicit AABBMesh(const indexed_triangle_set &tmesh, bool calculate_epsilon = false); + explicit AABBMesh(const TriangleMesh &mesh, bool calculate_epsilon = false); - IndexedMesh(const IndexedMesh& other); - IndexedMesh& operator=(const IndexedMesh&); + AABBMesh(const AABBMesh& other); + AABBMesh& operator=(const AABBMesh&); - IndexedMesh(IndexedMesh &&other); - IndexedMesh& operator=(IndexedMesh &&other); + AABBMesh(AABBMesh &&other); + AABBMesh& operator=(AABBMesh &&other); - ~IndexedMesh(); + ~AABBMesh(); - inline double ground_level() const { return m_ground_level + m_gnd_offset; } - inline void ground_level_offset(double o) { m_gnd_offset = o; } - inline double ground_level_offset() const { return m_gnd_offset; } + inline double ground_level() const { return m_ground_level /*+ m_gnd_offset*/; } +// inline void ground_level_offset(double o) { m_gnd_offset = o; } +// inline double ground_level_offset() const { return m_gnd_offset; } const std::vector& vertices() const; const std::vector& indices() const; @@ -73,15 +70,15 @@ public: // m_t holds a distance from m_source to the intersection. double m_t = infty(); int m_face_id = -1; - const IndexedMesh *m_mesh = nullptr; + const AABBMesh *m_mesh = nullptr; Vec3d m_dir; Vec3d m_source; Vec3d m_normal; - friend class IndexedMesh; + friend class AABBMesh; // A valid object of this class can only be obtained from // IndexedMesh::query_ray_hit method. - explicit inline hit_result(const IndexedMesh& em): m_mesh(&em) {} + explicit inline hit_result(const AABBMesh& em): m_mesh(&em) {} public: // This denotes no hit on the mesh. static inline constexpr double infty() { return std::numeric_limits::infinity(); } @@ -109,7 +106,7 @@ public: #ifdef SLIC3R_HOLE_RAYCASTER // Inform the object about location of holes // creates internal copy of the vector - void load_holes(const std::vector& holes) { + void load_holes(const std::vector& holes) { m_holes = holes; } @@ -118,7 +115,7 @@ public: // This function is currently not used anywhere, it was written when the // holes were subtracted on slices, that is, before we started using CGAL // to actually cut the holes into the mesh. - hit_result filter_hits(const std::vector& obj_hits) const; + hit_result filter_hits(const std::vector& obj_hits) const; #endif // Casting a ray on the mesh, returns the distance where the hit occures. @@ -138,16 +135,12 @@ public: Vec3d normal_by_face_id(int face_id) const; const indexed_triangle_set * get_triangle_mesh() const { return m_tm; } + + const VertexFaceIndex &vertex_face_index() const { return m_vfidx; } + const std::vector &face_neighbor_index() const { return m_fnidx; } }; -// Calculate the normals for the selected points (from 'points' set) on the -// mesh. This will call squared distance for each point. -PointSet normals(const PointSet& points, - const IndexedMesh& convert_mesh, - double eps = 0.05, // min distance from edges - std::function throw_on_cancel = [](){}, - const std::vector& selected_points = {}); -}} // namespace Slic3r::sla +} // namespace Slic3r::sla #endif // INDEXEDMESH_H diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 396eb0764..6710fa3c0 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -18,6 +18,9 @@ set(SLIC3R_SOURCES pchheader.cpp pchheader.hpp AStar.hpp + AABBTreeIndirect.hpp + AABBMesh.hpp + AABBMesh.cpp BoundingBox.cpp BoundingBox.hpp BridgeDetector.cpp @@ -172,6 +175,8 @@ set(SLIC3R_SOURCES ModelArrange.cpp MultiMaterialSegmentation.cpp MultiMaterialSegmentation.hpp + MeshNormals.hpp + MeshNormals.cpp CustomGCode.cpp CustomGCode.hpp Arrange.hpp @@ -279,12 +284,11 @@ set(SLIC3R_SOURCES SLA/SupportTreeBuilder.hpp SLA/SupportTreeMesher.hpp SLA/SupportTreeMesher.cpp - SLA/SupportTreeBuildsteps.hpp - SLA/SupportTreeBuildsteps.cpp + SLA/SupportTreeUtils.hpp SLA/SupportTreeBuilder.cpp - SLA/Concurrency.hpp SLA/SupportTree.hpp SLA/SupportTree.cpp + SLA/SupportTreeStrategies.hpp # SLA/SupportTreeIGL.cpp SLA/Rotfinder.hpp SLA/Rotfinder.cpp @@ -304,11 +308,11 @@ set(SLIC3R_SOURCES SLA/SupportPoint.hpp SLA/SupportPointGenerator.hpp SLA/SupportPointGenerator.cpp - SLA/IndexedMesh.hpp - SLA/IndexedMesh.cpp SLA/Clustering.hpp SLA/Clustering.cpp SLA/ReprojectPointsOnMesh.hpp + SLA/DefaultSupportTree.hpp + SLA/DefaultSupportTree.cpp ) add_library(libslic3r STATIC ${SLIC3R_SOURCES}) diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp index aa09a0d3e..4b685800d 100644 --- a/src/libslic3r/Geometry.hpp +++ b/src/libslic3r/Geometry.hpp @@ -558,6 +558,30 @@ inline bool is_rotation_ninety_degrees(const Vec3d &rotation) return is_rotation_ninety_degrees(rotation.x()) && is_rotation_ninety_degrees(rotation.y()) && is_rotation_ninety_degrees(rotation.z()); } +template +std::pair dir_to_spheric(const Vec<3, T> &n, T norm = 1.) +{ + T z = n.z(); + T r = norm; + T polar = std::acos(z / r); + T azimuth = std::atan2(n(1), n(0)); + return {polar, azimuth}; +} + +template +Vec<3, T> spheric_to_dir(double polar, double azimuth) +{ + return {T(std::cos(azimuth) * std::sin(polar)), + T(std::sin(azimuth) * std::sin(polar)), T(std::cos(polar))}; +} + +template +Vec<3, T> spheric_to_dir(const Pair &v) +{ + double plr = std::get<0>(v), azm = std::get<1>(v); + return spheric_to_dir(plr, azm); +} + } } // namespace Slicer::Geometry #endif diff --git a/src/libslic3r/MeshNormals.cpp b/src/libslic3r/MeshNormals.cpp new file mode 100644 index 000000000..a6c84dafa --- /dev/null +++ b/src/libslic3r/MeshNormals.cpp @@ -0,0 +1,158 @@ +#include "MeshNormals.hpp" + +namespace Slic3r { + +static bool point_on_edge(const Vec3d& p, const Vec3d& e1, const Vec3d& e2, + double epsSq = 0.05) +{ + using Line3D = Eigen::ParametrizedLine; + + auto line = Line3D::Through(e1, e2); + return line.squaredDistance(p) < epsSq; +} + +Vec3d get_normal(const AABBMesh &mesh, + const Vec3d &picking_point, + double eps) +{ + Vec3d ret = Vec3d::Zero(); + + int faceid = 0; + Vec3d p; + + mesh.squared_distance(picking_point, faceid, p); + + auto trindex = mesh.indices(faceid); + + const Vec3d &p1 = mesh.vertices(trindex(0)).cast(); + const Vec3d &p2 = mesh.vertices(trindex(1)).cast(); + const Vec3d &p3 = mesh.vertices(trindex(2)).cast(); + + // We should check if the point lies on an edge of the hosting + // triangle. If it does then all the other triangles using the + // same two points have to be searched and the final normal should + // be some kind of aggregation of the participating triangle + // normals. We should also consider the cases where the support + // point lies right on a vertex of its triangle. The procedure is + // the same, get the neighbor triangles and calculate an average + // normal. + + // Mark the vertex indices of the edge. ia and ib marks an edge. + // ic will mark a single vertex. + int vertex_idx = -1; + int edge_idx = -1; + double epsSq = eps * eps; + if ((p - p1).squaredNorm() < epsSq) { + vertex_idx = trindex(0); + } else if ((p - p2).squaredNorm() < epsSq) { + vertex_idx = trindex(1); + } else if ((p - p3).squaredNorm() < epsSq) { + vertex_idx = trindex(2); + } else if (point_on_edge(p, p1, p2, epsSq)) { + edge_idx = 0; + } else if (point_on_edge(p, p2, p3, epsSq)) { + edge_idx = 1; + } else if (point_on_edge(p, p1, p3, epsSq)) { + edge_idx = 2; + } + + // vector for the neigboring triangles including the detected one. + constexpr size_t MAX_EXPECTED_NEIGHBORS = 10; + boost::container::small_vector neigh; + + auto &vfidx = mesh.vertex_face_index(); + auto cmpfn = [](const Vec3d &v1, const Vec3d &v2) { return v1.sum() < v2.sum(); }; + auto eqfn = [](const Vec3d &n1, const Vec3d &n2) { + // Compare normals for equivalence. + // This is controvers stuff. + auto deq = [](double a, double b) { + return std::abs(a - b) < 1e-3; + }; + return deq(n1(X), n2(X)) && + deq(n1(Y), n2(Y)) && + deq(n1(Z), n2(Z)); + }; + + if (vertex_idx >= 0) { // The point is right on a vertex of the triangle + neigh.reserve(vfidx.count(vertex_idx)); + + auto from = vfidx.begin(vertex_idx); + auto to = vfidx.end(vertex_idx); + for (auto it = from; it != to; ++it) { + Vec3d nrm = mesh.normal_by_face_id(*it); + auto oit = std::lower_bound(neigh.begin(), neigh.end(), nrm, cmpfn); + if (oit == neigh.end() || !eqfn(*oit, nrm)) + neigh.insert(oit, mesh.normal_by_face_id(*it)); + } + } else if (edge_idx >= 0) { // the point is on and edge + size_t neighbor_face = mesh.face_neighbor_index()[faceid](edge_idx); + + neigh.emplace_back(mesh.normal_by_face_id(faceid)); + neigh.emplace_back(mesh.normal_by_face_id(neighbor_face)); + } + + if (!neigh.empty()) { // there were neighbors to count with + // sum up the normals and then normalize the result again. + // This unification seems to be enough. + Vec3d sumnorm(0, 0, 0); + sumnorm = std::accumulate(neigh.begin(), neigh.end(), sumnorm); + sumnorm.normalize(); + ret = sumnorm; + } else { // point lies safely within its triangle + Eigen::Vector3d U = p2 - p1; + Eigen::Vector3d V = p3 - p1; + ret = U.cross(V).normalized(); + } + + return ret; +} + +template +Eigen::MatrixXd normals(Ex ex_policy, + const PointSet &points, + const AABBMesh &mesh, + double eps, + std::function thr, // throw on cancel + const std::vector &pt_indices) +{ + if (points.rows() == 0 || mesh.vertices().empty() || + mesh.indices().empty()) + return {}; + + std::vector range = pt_indices; + if (range.empty()) { + range.resize(size_t(points.rows()), 0); + std::iota(range.begin(), range.end(), 0); + } + + PointSet ret(range.size(), 3); + + + execution::for_each(ex_policy, size_t(0), range.size(), + [&ret, &mesh, &points, thr, eps, &range](size_t ridx) { + thr(); + unsigned el = range[ridx]; + auto eidx = Eigen::Index(el); + auto picking_point = points.row(eidx); + + ret.row(ridx) = get_normal(mesh, picking_point, eps); + }); + + return ret; +} + +template Eigen::MatrixXd normals(ExecutionSeq policy, + const PointSet &points, + const AABBMesh &convert_mesh, + double eps, + std::function throw_on_cancel, + const std::vector &selected_points); + +template Eigen::MatrixXd normals(ExecutionTBB policy, + const PointSet &points, + const AABBMesh &convert_mesh, + double eps, + std::function throw_on_cancel, + const std::vector &selected_points); + +} // namespace Slic3r diff --git a/src/libslic3r/MeshNormals.hpp b/src/libslic3r/MeshNormals.hpp new file mode 100644 index 000000000..104f45c24 --- /dev/null +++ b/src/libslic3r/MeshNormals.hpp @@ -0,0 +1,52 @@ +#ifndef MESHNORMALS_HPP +#define MESHNORMALS_HPP + +#include "AABBMesh.hpp" + +#include "libslic3r/Execution/ExecutionSeq.hpp" +#include "libslic3r/Execution/ExecutionTBB.hpp" + +namespace Slic3r { + +// Get a good approximation of the normal for any picking point on the mesh. +// For points projecting to a face, this is the face normal, but when the +// picking point is on an edge or a vertex of the mesh, the normal is the +// normalized sum of each unique face normal (works nicely). The eps parameter +// gives a tolerance for how close a sample point has to be to an edge or +// vertex to start considering neighboring faces for the resulting normal. +Vec3d get_normal(const AABBMesh &mesh, + const Vec3d &picking_point, + double eps = 0.05); + +using PointSet = Eigen::MatrixXd; + +// Calculate the normals for the selected points (from 'points' set) on the +// mesh. This will call squared distance for each point. +template +Eigen::MatrixXd normals( + Ex ex_policy, + const PointSet &points, + const AABBMesh &convert_mesh, + double eps = 0.05, // min distance from edges + std::function throw_on_cancel = []() {}, + const std::vector &selected_points = {}); + +extern template Eigen::MatrixXd normals( + ExecutionSeq policy, + const PointSet &points, + const AABBMesh &convert_mesh, + double eps, + std::function throw_on_cancel, + const std::vector &selected_points); + +extern template Eigen::MatrixXd normals( + ExecutionTBB policy, + const PointSet &points, + const AABBMesh &convert_mesh, + double eps, + std::function throw_on_cancel, + const std::vector &selected_points); + +} // namespace Slic3r + +#endif // MESHNORMALS_HPP diff --git a/src/libslic3r/PrintConfig.cpp b/src/libslic3r/PrintConfig.cpp index b21ed5631..b39ce76d6 100644 --- a/src/libslic3r/PrintConfig.cpp +++ b/src/libslic3r/PrintConfig.cpp @@ -2,6 +2,8 @@ #include "Config.hpp" #include "I18N.hpp" +#include "SLA/SupportTree.hpp" + #include #include #include @@ -161,9 +163,9 @@ static const t_config_enum_values s_keys_map_SLADisplayOrientation = { CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(SLADisplayOrientation) static const t_config_enum_values s_keys_map_SLAPillarConnectionMode = { - {"zigzag", slapcmZigZag}, - {"cross", slapcmCross}, - {"dynamic", slapcmDynamic} + {"zigzag", int(SLAPillarConnectionMode::zigzag)}, + {"cross", int(SLAPillarConnectionMode::cross)}, + {"dynamic", int(SLAPillarConnectionMode::dynamic)} }; CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(SLAPillarConnectionMode) @@ -3498,14 +3500,14 @@ void PrintConfigDef::init_sla_params() " will automatically switch between the first two depending" " on the distance of the two pillars."); def->enum_keys_map = &ConfigOptionEnum::get_enum_values(); - def->enum_values.push_back("zigzag"); - def->enum_values.push_back("cross"); - def->enum_values.push_back("dynamic"); - def->enum_labels.push_back(L("Zig-Zag")); - def->enum_labels.push_back(L("Cross")); - def->enum_labels.push_back(L("Dynamic")); + def->enum_keys_map = &ConfigOptionEnum::get_enum_values(); + def->enum_values = ConfigOptionEnum::get_enum_names(); + def->enum_labels = ConfigOptionEnum::get_enum_names(); + def->enum_labels[0] = L("Zig-Zag"); + def->enum_labels[1] = L("Cross"); + def->enum_labels[2] = L("Dynamic"); def->mode = comAdvanced; - def->set_default_value(new ConfigOptionEnum(slapcmDynamic)); + def->set_default_value(new ConfigOptionEnum(SLAPillarConnectionMode::dynamic)); def = this->add("support_buildplate_only", coBool); def->label = L("Support on build plate only"); diff --git a/src/libslic3r/PrintConfig.hpp b/src/libslic3r/PrintConfig.hpp index 0c1060b7d..424746ec8 100644 --- a/src/libslic3r/PrintConfig.hpp +++ b/src/libslic3r/PrintConfig.hpp @@ -18,6 +18,7 @@ #include "libslic3r.h" #include "Config.hpp" +#include "SLA/SupportTreeStrategies.hpp" #include #include @@ -110,11 +111,8 @@ enum SLADisplayOrientation { sladoPortrait }; -enum SLAPillarConnectionMode { - slapcmZigZag, - slapcmCross, - slapcmDynamic -}; +using SLASupportTreeType = sla::SupportTreeType; +using SLAPillarConnectionMode = sla::PillarConnectionMode; enum BrimType { btNoBrim, diff --git a/src/libslic3r/SLA/Concurrency.hpp b/src/libslic3r/SLA/Concurrency.hpp deleted file mode 100644 index 7299101b3..000000000 --- a/src/libslic3r/SLA/Concurrency.hpp +++ /dev/null @@ -1,70 +0,0 @@ -#ifndef SLA_CONCURRENCY_H -#define SLA_CONCURRENCY_H - -// FIXME: Deprecated - -#include -#include - -namespace Slic3r { -namespace sla { - -// Set this to true to enable full parallelism in this module. -// Only the well tested parts will be concurrent if this is set to false. -const constexpr bool USE_FULL_CONCURRENCY = true; - -template struct _ccr {}; - -template<> struct _ccr -{ - using SpinningMutex = execution::SpinningMutex; - using BlockingMutex = execution::BlockingMutex; - - template - static void for_each(It from, It to, Fn &&fn, size_t granularity = 1) - { - execution::for_each(ex_tbb, from, to, std::forward(fn), granularity); - } - - template - static auto reduce(Args&&...args) - { - return execution::reduce(ex_tbb, std::forward(args)...); - } - - static size_t max_concurreny() - { - return execution::max_concurrency(ex_tbb); - } -}; - -template<> struct _ccr -{ - using SpinningMutex = execution::SpinningMutex; - using BlockingMutex = execution::BlockingMutex; - - template - static void for_each(It from, It to, Fn &&fn, size_t granularity = 1) - { - execution::for_each(ex_seq, from, to, std::forward(fn), granularity); - } - - template - static auto reduce(Args&&...args) - { - return execution::reduce(ex_seq, std::forward(args)...); - } - - static size_t max_concurreny() - { - return execution::max_concurrency(ex_seq); - } -}; - -using ccr = _ccr; -using ccr_seq = _ccr; -using ccr_par = _ccr; - -}} // namespace Slic3r::sla - -#endif // SLACONCURRENCY_H diff --git a/src/libslic3r/SLA/DefaultSupportTree.cpp b/src/libslic3r/SLA/DefaultSupportTree.cpp new file mode 100644 index 000000000..8634fc3bb --- /dev/null +++ b/src/libslic3r/SLA/DefaultSupportTree.cpp @@ -0,0 +1,983 @@ +#include "DefaultSupportTree.hpp" + +#include +#include +#include +#include + +namespace Slic3r { namespace sla { + +using Slic3r::opt::initvals; +using Slic3r::opt::bounds; +using Slic3r::opt::StopCriteria; +using Slic3r::opt::Optimizer; +using Slic3r::opt::AlgNLoptSubplex; +using Slic3r::opt::AlgNLoptGenetic; + +DefaultSupportTree::DefaultSupportTree(SupportTreeBuilder & builder, + const SupportableMesh &sm) + : m_sm(sm) + , m_support_nmls(sm.pts.size(), 3) + , m_builder(builder) + , m_points(sm.pts.size(), 3) + , m_thr(builder.ctl().cancelfn) +{ + // Prepare the support points in Eigen/IGL format as well, we will use + // it mostly in this form. + + long i = 0; + for (const SupportPoint &sp : m_sm.pts) { + m_points.row(i).x() = double(sp.pos.x()); + m_points.row(i).y() = double(sp.pos.y()); + m_points.row(i).z() = double(sp.pos.z()); + ++i; + } +} + +bool DefaultSupportTree::execute(SupportTreeBuilder &builder, + const SupportableMesh &sm) +{ + if(sm.pts.empty()) return false; + + DefaultSupportTree alg(builder, sm); + + // Let's define the individual steps of the processing. We can experiment + // later with the ordering and the dependencies between them. + enum Steps { + BEGIN, + PINHEADS, + CLASSIFY, + ROUTING_GROUND, + ROUTING_NONGROUND, + CASCADE_PILLARS, + MERGE_RESULT, + DONE, + ABORT, + NUM_STEPS + //... + }; + + // Collect the algorithm steps into a nice sequence + std::array, NUM_STEPS> program = { + [] () { + // Begin... + // Potentially clear up the shared data (not needed for now) + }, + + std::bind(&DefaultSupportTree::add_pinheads, &alg), + + std::bind(&DefaultSupportTree::classify, &alg), + + std::bind(&DefaultSupportTree::routing_to_ground, &alg), + + std::bind(&DefaultSupportTree::routing_to_model, &alg), + + std::bind(&DefaultSupportTree::interconnect_pillars, &alg), + + std::bind(&DefaultSupportTree::merge_result, &alg), + + [] () { + // Done + }, + + [] () { + // Abort + } + }; + + Steps pc = BEGIN; + + if(sm.cfg.ground_facing_only) { + program[ROUTING_NONGROUND] = []() { + BOOST_LOG_TRIVIAL(info) + << "Skipping model-facing supports as requested."; + }; + } + + // Let's define a simple automaton that will run our program. + auto progress = [&builder, &pc] () { + static const std::array stepstr { + "Starting", + "Generate pinheads", + "Classification", + "Routing to ground", + "Routing supports to model surface", + "Interconnecting pillars", + "Merging support mesh", + "Done", + "Abort" + }; + + static const std::array stepstate { + 0, + 30, + 50, + 60, + 70, + 80, + 99, + 100, + 0 + }; + + if(builder.ctl().stopcondition()) pc = ABORT; + + switch(pc) { + case BEGIN: 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 = CASCADE_PILLARS; break; + case CASCADE_PILLARS: pc = MERGE_RESULT; break; + case MERGE_RESULT: pc = DONE; break; + case DONE: + case ABORT: break; + default: ; + } + + builder.ctl().statuscb(stepstate[pc], stepstr[pc]); + }; + + // Just here we run the computation... + while(pc < DONE) { + progress(); + program[pc](); + } + + return pc == ABORT; +} + +AABBMesh::hit_result DefaultSupportTree::pinhead_mesh_intersect( + const Vec3d &s, + const Vec3d &dir, + double r_pin, + double r_back, + double width, + double sd) +{ + return sla::pinhead_mesh_hit(suptree_ex_policy, m_sm.emesh, s, dir, r_pin, r_back, width, sd); +} + +AABBMesh::hit_result DefaultSupportTree::bridge_mesh_intersect( + const Vec3d &src, const Vec3d &dir, double r, double sd) +{ + return sla::beam_mesh_hit(suptree_ex_policy, m_sm.emesh, {src, dir, r}, sd); +} + +bool DefaultSupportTree::interconnect(const Pillar &pillar, + const Pillar &nextpillar) +{ + // We need to get the starting point of the zig-zag pattern. We have to + // be aware that the two head junctions are at different heights. We + // may start from the lowest junction and call it a day but this + // strategy would leave unconnected a lot of pillar duos where the + // shorter pillar is too short to start a new bridge but the taller + // pillar could still be bridged with the shorter one. + bool was_connected = false; + + Vec3d supper = pillar.startpoint(); + Vec3d slower = nextpillar.startpoint(); + Vec3d eupper = pillar.endpoint(); + Vec3d elower = nextpillar.endpoint(); + + double zmin = ground_level(this->m_sm) + m_sm.cfg.base_height_mm; + eupper.z() = std::max(eupper.z(), zmin); + elower.z() = std::max(elower.z(), zmin); + + // The usable length of both pillars should be positive + if(slower.z() - elower.z() < 0) return false; + if(supper.z() - eupper.z() < 0) return false; + + double pillar_dist = distance(Vec2d{slower.x(), slower.y()}, + Vec2d{supper.x(), supper.y()}); + double bridge_distance = pillar_dist / std::cos(-m_sm.cfg.bridge_slope); + double zstep = pillar_dist * std::tan(-m_sm.cfg.bridge_slope); + + if(pillar_dist < 2 * m_sm.cfg.head_back_radius_mm || + pillar_dist > m_sm.cfg.max_pillar_link_distance_mm) return false; + + if(supper.z() < slower.z()) supper.swap(slower); + if(eupper.z() < elower.z()) eupper.swap(elower); + + double startz = 0, endz = 0; + + startz = slower.z() - zstep < supper.z() ? slower.z() - zstep : slower.z(); + endz = eupper.z() + zstep > elower.z() ? eupper.z() + zstep : eupper.z(); + + if(slower.z() - eupper.z() < std::abs(zstep)) { + // no space for even one cross + + // Get max available space + startz = std::min(supper.z(), slower.z() - zstep); + endz = std::max(eupper.z() + zstep, elower.z()); + + // Align to center + double available_dist = (startz - endz); + double rounds = std::floor(available_dist / std::abs(zstep)); + startz -= 0.5 * (available_dist - rounds * std::abs(zstep)); + } + + auto pcm = m_sm.cfg.pillar_connection_mode; + bool docrosses = + pcm == PillarConnectionMode::cross || + (pcm == PillarConnectionMode::dynamic && + pillar_dist > 2*m_sm.cfg.base_radius_mm); + + // 'sj' means starting junction, 'ej' is the end junction of a bridge. + // They will be swapped in every iteration thus the zig-zag pattern. + // According to a config parameter, a second bridge may be added which + // results in a cross connection between the pillars. + Vec3d sj = supper, ej = slower; sj.z() = startz; ej.z() = sj.z() + zstep; + + // TODO: This is a workaround to not have a faulty last bridge + while(ej.z() >= eupper.z() /*endz*/) { + if(bridge_mesh_distance(sj, dirv(sj, ej), pillar.r_start) >= bridge_distance) + { + m_builder.add_crossbridge(sj, ej, pillar.r_start); + was_connected = true; + } + + // double bridging: (crosses) + if(docrosses) { + Vec3d sjback(ej.x(), ej.y(), sj.z()); + Vec3d ejback(sj.x(), sj.y(), ej.z()); + if (sjback.z() <= slower.z() && ejback.z() >= eupper.z() && + bridge_mesh_distance(sjback, dirv(sjback, ejback), + pillar.r_start) >= bridge_distance) { + // need to check collision for the cross stick + m_builder.add_crossbridge(sjback, ejback, pillar.r_start); + was_connected = true; + } + } + + sj.swap(ej); + ej.z() = sj.z() + zstep; + } + + return was_connected; +} + +bool DefaultSupportTree::connect_to_nearpillar(const Head &head, + long nearpillar_id) +{ + auto nearpillar = [this, nearpillar_id]() -> const Pillar& { + return m_builder.pillar(nearpillar_id); + }; + + if (m_builder.bridgecount(nearpillar()) > m_sm.cfg.max_bridges_on_pillar) + return false; + + Vec3d headjp = head.junction_point(); + Vec3d nearjp_u = nearpillar().startpoint(); + Vec3d nearjp_l = nearpillar().endpoint(); + + double r = head.r_back_mm; + double d2d = distance(to_2d(headjp), to_2d(nearjp_u)); + double d3d = distance(headjp, nearjp_u); + + double hdiff = nearjp_u.z() - headjp.z(); + double slope = std::atan2(hdiff, d2d); + + Vec3d bridgestart = headjp; + Vec3d bridgeend = nearjp_u; + double max_len = r * m_sm.cfg.max_bridge_length_mm / m_sm.cfg.head_back_radius_mm; + double max_slope = m_sm.cfg.bridge_slope; + double zdiff = 0.0; + + // check the default situation if feasible for a bridge + if(d3d > max_len || slope > -max_slope) { + // not feasible to connect the two head junctions. We have to search + // for a suitable touch point. + + double Zdown = headjp.z() + d2d * std::tan(-max_slope); + Vec3d touchjp = bridgeend; touchjp.z() = Zdown; + double D = distance(headjp, touchjp); + zdiff = Zdown - nearjp_u.z(); + + if(zdiff > 0) { + Zdown -= zdiff; + bridgestart.z() -= zdiff; + touchjp.z() = Zdown; + + double t = bridge_mesh_distance(headjp, DOWN, r); + + // We can't insert a pillar under the source head to connect + // with the nearby pillar's starting junction + if(t < zdiff) return false; + } + + if(Zdown <= nearjp_u.z() && Zdown >= nearjp_l.z() && D < max_len) + bridgeend.z() = Zdown; + else + return false; + } + + // There will be a minimum distance from the ground where the + // bridge is allowed to connect. This is an empiric value. + double minz = ground_level(m_sm) + 4 * head.r_back_mm; + if(bridgeend.z() < minz) return false; + + double t = bridge_mesh_distance(bridgestart, dirv(bridgestart, bridgeend), r); + + // Cannot insert the bridge. (further search might not worth the hassle) + if(t < distance(bridgestart, bridgeend)) return false; + + std::lock_guard lk(m_bridge_mutex); + + if (m_builder.bridgecount(nearpillar()) < m_sm.cfg.max_bridges_on_pillar) { + // A partial pillar is needed under the starting head. + if(zdiff > 0) { + m_builder.add_pillar(head.id, headjp.z() - bridgestart.z()); + m_builder.add_junction(bridgestart, r); + m_builder.add_bridge(bridgestart, bridgeend, r); + } else { + m_builder.add_bridge(head.id, bridgeend); + } + + m_builder.increment_bridges(nearpillar()); + } else return false; + + return true; +} + +bool DefaultSupportTree::create_ground_pillar(const Vec3d &hjp, + const Vec3d &sourcedir, + double radius, + long head_id) +{ + auto [ret, pillar_id] = sla::create_ground_pillar(suptree_ex_policy, + m_builder, m_sm, hjp, + sourcedir, radius, radius, + head_id); + + if (pillar_id >= 0) // Save the pillar endpoint in the spatial index + m_pillar_index.guarded_insert(m_builder.pillar(pillar_id).endpt, + unsigned(pillar_id)); + + return ret; +} + +void DefaultSupportTree::add_pinheads() +{ + // The minimum distance for two support points to remain valid. + const double /*constexpr*/ D_SP = 0.1; + + // Get the points that are too close to each other and keep only the + // first one + auto aliases = cluster(m_points, D_SP, 2); + + PtIndices filtered_indices; + filtered_indices.reserve(aliases.size()); + m_iheads.reserve(aliases.size()); + m_iheadless.reserve(aliases.size()); + for(auto& a : aliases) { + // Here we keep only the front point of the cluster. + filtered_indices.emplace_back(a.front()); + } + + // calculate the normals to the triangles for filtered points + auto nmls = normals(suptree_ex_policy, m_points, m_sm.emesh, + m_sm.cfg.head_front_radius_mm, m_thr, + filtered_indices); + + // 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. + + auto heads = reserve_vector(m_sm.pts.size()); + for (const SupportPoint &sp : m_sm.pts) { + m_thr(); + heads.emplace_back( + NaNd, + sp.head_front_radius, + 0., + m_sm.cfg.head_penetration_mm, + Vec3d::Zero(), // dir + sp.pos.cast() // displacement + ); + } + + std::function filterfn; + filterfn = [this, &nmls, &heads, &filterfn](unsigned fidx, size_t i, double back_r) { + m_thr(); + + Vec3d n = nmls.row(Eigen::Index(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. + + auto [polar, azimuth] = dir_to_spheric(n); + + // skip if the tilt is not sane + if (polar < PI - m_sm.cfg.normal_cutoff_angle) return; + + // We saturate the polar angle to 3pi/4 + polar = std::max(polar, PI - m_sm.cfg.bridge_slope); + + // save the head (pinpoint) position + Vec3d hp = m_points.row(fidx); + + double lmin = m_sm.cfg.head_width_mm, lmax = lmin; + + if (back_r < m_sm.cfg.head_back_radius_mm) { + lmin = 0., lmax = m_sm.cfg.head_penetration_mm; + } + + // The distance needed for a pinhead to not collide with model. + double w = lmin + 2 * back_r + 2 * m_sm.cfg.head_front_radius_mm - + m_sm.cfg.head_penetration_mm; + + double pin_r = double(m_sm.pts[fidx].head_front_radius); + + // Reassemble the now corrected normal + auto nn = spheric_to_dir(polar, azimuth).normalized(); + + // check available distance + AABBMesh::hit_result t = pinhead_mesh_intersect(hp, nn, pin_r, + back_r, w); + + if (t.distance() < w) { + // Let's try to optimize this angle, there might be a + // viable normal that doesn't collide with the model + // geometry and its very close to the default. + + Optimizer solver(get_criteria(m_sm.cfg)); + solver.seed(0); // we want deterministic behavior + + auto oresult = solver.to_max().optimize( + [this, pin_r, back_r, hp](const opt::Input<3> &input) + { + auto &[plr, azm, l] = input; + + auto dir = spheric_to_dir(plr, azm).normalized(); + + return pinhead_mesh_intersect( + hp, dir, pin_r, back_r, l).distance(); + }, + initvals({polar, azimuth, (lmin + lmax) / 2.}), // start with what we have + bounds({ + {PI - m_sm.cfg.bridge_slope, PI}, // Must not exceed the slope limit + {-PI, PI}, // azimuth can be a full search + {lmin, lmax} + })); + + if(oresult.score > w) { + polar = std::get<0>(oresult.optimum); + azimuth = std::get<1>(oresult.optimum); + nn = spheric_to_dir(polar, azimuth).normalized(); + lmin = std::get<2>(oresult.optimum); + t = AABBMesh::hit_result(oresult.score); + } + } + + if (t.distance() > w && hp.z() + w * nn.z() >= ground_level(m_sm)) { + Head &h = heads[fidx]; + h.id = fidx; + h.dir = nn; + h.width_mm = lmin; + h.r_back_mm = back_r; + } else if (back_r > m_sm.cfg.head_fallback_radius_mm) { + filterfn(fidx, i, m_sm.cfg.head_fallback_radius_mm); + } + }; + + execution::for_each( + suptree_ex_policy, size_t(0), filtered_indices.size(), + [this, &filterfn, &filtered_indices](size_t i) { + filterfn(filtered_indices[i], i, m_sm.cfg.head_back_radius_mm); + }, + execution::max_concurrency(suptree_ex_policy)); + + for (size_t i = 0; i < heads.size(); ++i) + if (heads[i].is_valid()) { + m_builder.add_head(i, heads[i]); + m_iheads.emplace_back(i); + } + + m_thr(); +} + +void DefaultSupportTree::classify() +{ + // We should first get the heads that reach the ground directly + PtIndices ground_head_indices; + ground_head_indices.reserve(m_iheads.size()); + m_iheads_onmodel.reserve(m_iheads.size()); + + // First we decide which heads reach the ground and can be full + // pillars and which shall be connected to the model surface (or + // search a suitable path around the surface that leads to the + // ground -- TODO) + for(unsigned i : m_iheads) { + m_thr(); + + Head &head = m_builder.head(i); + double r = head.r_back_mm; + Vec3d headjp = head.junction_point(); + + // collision check + auto hit = bridge_mesh_intersect(headjp, DOWN, r); + + if(std::isinf(hit.distance())) ground_head_indices.emplace_back(i); + else if(m_sm.cfg.ground_facing_only) head.invalidate(); + else m_iheads_onmodel.emplace_back(i); + + m_head_to_ground_scans[i] = hit; + } + + // 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 pointfn = [this](unsigned i) { + return m_builder.head(i).junction_point(); + }; + + auto predicate = [this](const PointIndexEl &e1, + const PointIndexEl &e2) { + double d2d = distance(to_2d(e1.first), to_2d(e2.first)); + double d3d = distance(e1.first, e2.first); + return d2d < 2 * m_sm.cfg.base_radius_mm + && d3d < m_sm.cfg.max_bridge_length_mm; + }; + + m_pillar_clusters = cluster(ground_head_indices, pointfn, predicate, + m_sm.cfg.max_bridges_on_pillar); +} + +void DefaultSupportTree::routing_to_ground() +{ + ClusterEl cl_centroids; + cl_centroids.reserve(m_pillar_clusters.size()); + + for (auto &cl : m_pillar_clusters) { + m_thr(); + + // 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. + + if (cl.empty()) continue; + + // get the current cluster centroid + auto & thr = m_thr; + const auto &points = m_points; + + long lcid = cluster_centroid( + cl, [&points](size_t idx) { return points.row(long(idx)); }, + [thr](const Vec3d &p1, const Vec3d &p2) { + thr(); + return distance(Vec2d(p1.x(), p1.y()), Vec2d(p2.x(), p2.y())); + }); + + assert(lcid >= 0); + unsigned hid = cl[size_t(lcid)]; // Head ID + + cl_centroids.emplace_back(hid); + + Head &h = m_builder.head(hid); + + if (!create_ground_pillar(h.junction_point(), h.dir, h.r_back_mm, h.id)) { + BOOST_LOG_TRIVIAL(warning) + << "Pillar cannot be created for support point id: " << hid; + m_iheads_onmodel.emplace_back(h.id); + continue; + } + } + + // 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. + size_t ci = 0; + for (auto cl : m_pillar_clusters) { + m_thr(); + + auto cidx = cl_centroids[ci++]; + + auto q = m_pillar_index.query(m_builder.head(cidx).junction_point(), 1); + if (!q.empty()) { + long centerpillarID = q.front().second; + for (auto c : cl) { + m_thr(); + if (c == cidx) continue; + + auto &sidehead = m_builder.head(c); + + if (!connect_to_nearpillar(sidehead, centerpillarID) && + !search_pillar_and_connect(sidehead)) { + Vec3d pstart = sidehead.junction_point(); + // Vec3d pend = Vec3d{pstart.x(), pstart.y(), gndlvl}; + // Could not find a pillar, create one + create_ground_pillar(pstart, sidehead.dir, sidehead.r_back_mm, sidehead.id); + } + } + } + } +} + +bool DefaultSupportTree::connect_to_ground(Head &head) +{ + auto [ret, pillar_id] = sla::search_ground_route(suptree_ex_policy, + m_builder, m_sm, + {head.junction_point(), + head.r_back_mm}, + head.r_back_mm, + head.dir); + + if (pillar_id >= 0) { + // Save the pillar endpoint in the spatial index + m_pillar_index.guarded_insert(m_builder.pillar(pillar_id).endpt, + unsigned(pillar_id)); + + head.pillar_id = pillar_id; + } + + return ret; +} + +bool DefaultSupportTree::connect_to_model_body(Head &head) +{ + if (head.id <= SupportTreeNode::ID_UNSET) return false; + + auto it = m_head_to_ground_scans.find(unsigned(head.id)); + if (it == m_head_to_ground_scans.end()) return false; + + auto &hit = it->second; + + if (!hit.is_hit()) { + // TODO scan for potential anchor points on model surface + return false; + } + + Vec3d hjp = head.junction_point(); + double zangle = std::asin(hit.direction().z()); + zangle = std::max(zangle, PI/4); + double h = std::sin(zangle) * head.fullwidth(); + + // The width of the tail head that we would like to have... + h = std::min(hit.distance() - head.r_back_mm, h); + + // If this is a mini pillar dont bother with the tail width, can be 0. + if (head.r_back_mm < m_sm.cfg.head_back_radius_mm) h = std::max(h, 0.); + else if (h <= 0.) return false; + + Vec3d endp{hjp.x(), hjp.y(), hjp.z() - hit.distance() + h}; + auto center_hit = m_sm.emesh.query_ray_hit(hjp, DOWN); + + double hitdiff = center_hit.distance() - hit.distance(); + Vec3d hitp = std::abs(hitdiff) < 2*head.r_back_mm? + center_hit.position() : hit.position(); + + long pillar_id = m_builder.add_pillar(head.id, hjp.z() - endp.z()); + Pillar &pill = m_builder.pillar(pillar_id); + + Vec3d taildir = endp - hitp; + double dist = (hitp - endp).norm() + m_sm.cfg.head_penetration_mm; + double w = dist - 2 * head.r_pin_mm - head.r_back_mm; + + if (w < 0.) { + BOOST_LOG_TRIVIAL(warning) << "Pinhead width is negative!"; + w = 0.; + } + + m_builder.add_anchor(head.r_back_mm, head.r_pin_mm, w, + m_sm.cfg.head_penetration_mm, taildir, hitp); + + m_pillar_index.guarded_insert(pill.endpoint(), pill.id); + + return true; +} + +bool DefaultSupportTree::search_pillar_and_connect(const Head &source) +{ + // Hope that a local copy takes less time than the whole search loop. + // We also need to remove elements progressively from the copied index. + PointIndex spindex = m_pillar_index.guarded_clone(); + + long nearest_id = SupportTreeNode::ID_UNSET; + + Vec3d querypt = source.junction_point(); + + while(nearest_id < 0 && !spindex.empty()) { m_thr(); + // 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(querypt.x(), querypt.y(), ground_level(m_sm)); + auto qres = spindex.nearest(qp, 1); + if(qres.empty()) break; + + auto ne = qres.front(); + nearest_id = ne.second; + + if(nearest_id >= 0) { + if (size_t(nearest_id) < m_builder.pillarcount()) { + if(!connect_to_nearpillar(source, nearest_id) || + m_builder.pillar(nearest_id).r_start < source.r_back_mm) { + nearest_id = SupportTreeNode::ID_UNSET; // continue searching + spindex.remove(ne); // without the current pillar + } + } + } + } + + return nearest_id >= 0; +} + +void DefaultSupportTree::routing_to_model() +{ + // We need to check if there is an easy way out to the bed surface. + // If it can be routed there with a bridge shorter than + // min_bridge_distance. + + execution::for_each( + suptree_ex_policy, m_iheads_onmodel.begin(), m_iheads_onmodel.end(), + [this](const unsigned idx) { + m_thr(); + + auto &head = m_builder.head(idx); + + // Search nearby pillar + if (search_pillar_and_connect(head)) { return; } + + // Cannot connect to nearby pillar. We will try to search for + // a route to the ground. + if (connect_to_ground(head)) { return; } + + // No route to the ground, so connect to the model body as a last resort + if (connect_to_model_body(head)) { return; } + + // We have failed to route this head. + BOOST_LOG_TRIVIAL(warning) + << "Failed to route model facing support point. ID: " << idx; + + head.invalidate(); + }, + execution::max_concurrency(suptree_ex_policy)); +} + +void DefaultSupportTree::interconnect_pillars() +{ + // Now comes the algorithm that connects pillars with each other. + // Ideally every pillar should be connected with at least one of its + // neighbors if that neighbor is within max_pillar_link_distance + + // Pillars with height exceeding H1 will require at least one neighbor + // to connect with. Height exceeding H2 require two neighbors. + double H1 = m_sm.cfg.max_solo_pillar_height_mm; + double H2 = m_sm.cfg.max_dual_pillar_height_mm; + double d = m_sm.cfg.max_pillar_link_distance_mm; + + //A connection between two pillars only counts if the height ratio is + // bigger than 50% + double min_height_ratio = 0.5; + + std::set pairs; + + // A function to connect one pillar with its neighbors. THe number of + // neighbors is given in the configuration. This function if called + // for every pillar in the pillar index. A pair of pillar will not + // be connected multiple times this is ensured by the 'pairs' set which + // remembers the processed pillar pairs + auto cascadefn = + [this, d, &pairs, min_height_ratio, H1] (const PointIndexEl& el) + { + Vec3d qp = el.first; // endpoint of the pillar + + const Pillar& pillar = m_builder.pillar(el.second); // actual pillar + + // Get the max number of neighbors a pillar should connect to + unsigned neighbors = m_sm.cfg.pillar_cascade_neighbors; + + // connections are already enough for the pillar + if(pillar.links >= neighbors) return; + + double max_d = d * pillar.r_start / m_sm.cfg.head_back_radius_mm; + // Query all remaining points within reach + auto qres = m_pillar_index.query([qp, max_d](const PointIndexEl& e){ + return distance(e.first, qp) < max_d; + }); + + // sort the result by distance (have to check if this is needed) + std::sort(qres.begin(), qres.end(), + [qp](const PointIndexEl& e1, const PointIndexEl& e2){ + return distance(e1.first, qp) < distance(e2.first, qp); + }); + + for(auto& re : qres) { // process the queried neighbors + + if(re.second == el.second) continue; // Skip self + + auto a = el.second, b = re.second; + + // Get unique hash for the given pair (order doesn't matter) + auto hashval = pairhash(a, b); + + // Search for the pair amongst the remembered pairs + if(pairs.find(hashval) != pairs.end()) continue; + + const Pillar& neighborpillar = m_builder.pillar(re.second); + + // this neighbor is occupied, skip + if (neighborpillar.links >= neighbors) continue; + if (neighborpillar.r_start < pillar.r_start) continue; + + if(interconnect(pillar, neighborpillar)) { + pairs.insert(hashval); + + // If the interconnection length between the two pillars is + // less than 50% of the longer pillar's height, don't count + if(pillar.height < H1 || + neighborpillar.height / pillar.height > min_height_ratio) + m_builder.increment_links(pillar); + + if(neighborpillar.height < H1 || + pillar.height / neighborpillar.height > min_height_ratio) + m_builder.increment_links(neighborpillar); + + } + + // connections are enough for one pillar + if(pillar.links >= neighbors) break; + } + }; + + // Run the cascade for the pillars in the index + m_pillar_index.foreach(cascadefn); + + // We would be done here if we could allow some pillars to not be + // connected with any neighbors. But this might leave the support tree + // unprintable. + // + // The current solution is to insert additional pillars next to these + // lonely pillars. One or even two additional pillar might get inserted + // depending on the length of the lonely pillar. + + size_t pillarcount = m_builder.pillarcount(); + + // Again, go through all pillars, this time in the whole support tree + // not just the index. + for(size_t pid = 0; pid < pillarcount; pid++) { + auto pillar = [this, pid]() { return m_builder.pillar(pid); }; + + // Decide how many additional pillars will be needed: + + unsigned needpillars = 0; + if (pillar().bridges > m_sm.cfg.max_bridges_on_pillar) + needpillars = 3; + else if (pillar().links < 2 && pillar().height > H2) { + // Not enough neighbors to support this pillar + needpillars = 2; + } else if (pillar().links < 1 && pillar().height > H1) { + // No neighbors could be found and the pillar is too long. + needpillars = 1; + } + + needpillars = std::max(pillar().links, needpillars) - pillar().links; + if (needpillars == 0) continue; + + // Search for new pillar locations: + + bool found = false; + double alpha = 0; // goes to 2Pi + double r = 2 * m_sm.cfg.base_radius_mm; + Vec3d pillarsp = pillar().startpoint(); + + // temp value for starting point detection + Vec3d sp(pillarsp.x(), pillarsp.y(), pillarsp.z() - r); + + // A vector of bool for placement feasbility + std::vector canplace(needpillars, false); + std::vector spts(needpillars); // vector of starting points + + double gnd = ground_level(m_sm); + double min_dist = m_sm.cfg.pillar_base_safety_distance_mm + + m_sm.cfg.base_radius_mm + EPSILON; + + while(!found && alpha < 2*PI) { + for (unsigned n = 0; + n < needpillars && (!n || canplace[n - 1]); + n++) + { + double a = alpha + n * PI / 3; + Vec3d s = sp; + s.x() += std::cos(a) * r; + s.y() += std::sin(a) * r; + spts[n] = s; + + // Check the path vertically down + Vec3d check_from = s + Vec3d{0., 0., pillar().r_start}; + auto hr = bridge_mesh_intersect(check_from, DOWN, pillar().r_start); + Vec3d gndsp{s.x(), s.y(), gnd}; + + // If the path is clear, check for pillar base collisions + canplace[n] = std::isinf(hr.distance()) && + std::sqrt(m_sm.emesh.squared_distance(gndsp)) > + min_dist; + } + + found = std::all_of(canplace.begin(), canplace.end(), + [](bool v) { return v; }); + + // 20 angles will be tried... + alpha += 0.1 * PI; + } + + std::vector newpills; + newpills.reserve(needpillars); + + if (found) + for (unsigned n = 0; n < needpillars; n++) { + Vec3d s = spts[n]; + Pillar p(Vec3d{s.x(), s.y(), gnd}, s.z() - gnd, pillar().r_start); + + if (interconnect(pillar(), p)) { + Pillar &pp = m_builder.pillar(m_builder.add_pillar(p)); + + add_pillar_base(pp.id); + + m_pillar_index.insert(pp.endpoint(), unsigned(pp.id)); + + m_builder.add_junction(s, pillar().r_start); + double t = bridge_mesh_distance(pillarsp, dirv(pillarsp, s), + pillar().r_start); + if (distance(pillarsp, s) < t) + m_builder.add_bridge(pillarsp, s, pillar().r_start); + + if (pillar().endpoint().z() > ground_level(m_sm) + pillar().r_start) + m_builder.add_junction(pillar().endpoint(), pillar().r_start); + + newpills.emplace_back(pp.id); + m_builder.increment_links(pillar()); + m_builder.increment_links(pp); + } + } + + if(!newpills.empty()) { + for(auto it = newpills.begin(), nx = std::next(it); + nx != newpills.end(); ++it, ++nx) { + const Pillar& itpll = m_builder.pillar(*it); + const Pillar& nxpll = m_builder.pillar(*nx); + if(interconnect(itpll, nxpll)) { + m_builder.increment_links(itpll); + m_builder.increment_links(nxpll); + } + } + + m_pillar_index.foreach(cascadefn); + } + } +} + +}} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.hpp b/src/libslic3r/SLA/DefaultSupportTree.hpp similarity index 55% rename from src/libslic3r/SLA/SupportTreeBuildsteps.hpp rename to src/libslic3r/SLA/DefaultSupportTree.hpp index 478e2df1f..ff0269978 100644 --- a/src/libslic3r/SLA/SupportTreeBuildsteps.hpp +++ b/src/libslic3r/SLA/DefaultSupportTree.hpp @@ -1,124 +1,18 @@ -#ifndef SLASUPPORTTREEALGORITHM_H -#define SLASUPPORTTREEALGORITHM_H +#ifndef LEGACYSUPPORTTREE_HPP +#define LEGACYSUPPORTTREE_HPP -#include -#include +#include "SupportTreeUtils.hpp" -#include -#include #include +#include -namespace Slic3r { -namespace sla { +namespace Slic3r { namespace sla { -// The minimum distance for two support points to remain valid. -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 -}; - -inline Vec2d to_vec2(const Vec3d &v3) { return {v3(X), v3(Y)}; } - -inline std::pair dir_to_spheric(const Vec3d &n, double norm = 1.) -{ - double z = n.z(); - double r = norm; - double polar = std::acos(z / r); - double azimuth = std::atan2(n(1), n(0)); - return {polar, azimuth}; -} - -inline Vec3d spheric_to_dir(double polar, double azimuth) -{ - return {std::cos(azimuth) * std::sin(polar), - std::sin(azimuth) * std::sin(polar), std::cos(polar)}; -} - -inline Vec3d spheric_to_dir(const std::tuple &v) -{ - auto [plr, azm] = v; - return spheric_to_dir(plr, azm); -} - -inline Vec3d spheric_to_dir(const std::pair &v) -{ - return spheric_to_dir(v.first, v.second); -} - -inline Vec3d spheric_to_dir(const std::array &v) -{ - return spheric_to_dir(v[0], v[1]); -} - -// Give points on a 3D ring with given center, radius and orientation -// method based on: -// https://math.stackexchange.com/questions/73237/parametric-equation-of-a-circle-in-3d-space -template -class PointRing { - std::array m_phis; - - // Two vectors that will be perpendicular to each other and to the - // axis. Values for a(X) and a(Y) are now arbitrary, a(Z) is just a - // placeholder. - // a and b vectors are perpendicular to the ring direction and to each other. - // Together they define the plane where we have to iterate with the - // given angles in the 'm_phis' vector - Vec3d a = {0, 1, 0}, b; - double m_radius = 0.; - - static inline bool constexpr is_one(double val) - { - return std::abs(std::abs(val) - 1) < 1e-20; - } - -public: - - PointRing(const Vec3d &n) - { - m_phis = linspace_array(0., 2 * PI); - - // We have to address the case when the direction vector v (same as - // dir) is coincident with one of the world axes. In this case two of - // its components will be completely zero and one is 1.0. Our method - // becomes dangerous here due to division with zero. Instead, vector - // 'a' can be an element-wise rotated version of 'v' - if(is_one(n(X)) || is_one(n(Y)) || is_one(n(Z))) { - a = {n(Z), n(X), n(Y)}; - b = {n(Y), n(Z), n(X)}; - } - else { - a(Z) = -(n(Y)*a(Y)) / n(Z); a.normalize(); - b = a.cross(n); - } - } - - Vec3d get(size_t idx, const Vec3d src, double r) const - { - double phi = m_phis[idx]; - double sinphi = std::sin(phi); - double cosphi = std::cos(phi); - - double rpscos = r * cosphi; - double rpssin = r * sinphi; - - // Point on the sphere - return {src(X) + rpscos * a(X) + rpssin * b(X), - src(Y) + rpscos * a(Y) + rpssin * b(Y), - src(Z) + rpscos * a(Z) + rpssin * b(Z)}; - } -}; - -//IndexedMesh::hit_result query_hit(const SupportableMesh &msh, const Bridge &br, double safety_d = std::nan("")); -//IndexedMesh::hit_result query_hit(const SupportableMesh &msh, const Head &br, double safety_d = std::nan("")); - -inline Vec3d dirv(const Vec3d& startp, const Vec3d& endp) { - return (endp - startp).normalized(); -} +inline constexpr const auto &suptree_ex_policy = ex_tbb; class PillarIndex { PointIndex m_index; - using Mutex = ccr::BlockingMutex; + using Mutex = execution::BlockingMutex; mutable Mutex m_mutex; public: @@ -161,45 +55,19 @@ public: } }; -// Helper function for pillar interconnection where pairs of already connected -// pillars should be checked for not to be processed again. This can be done -// in constant time with a set of hash values uniquely representing a pair of -// integers. The order of numbers within the pair should not matter, it has -// the same unique hash. The hash value has to have twice as many bits as the -// arguments need. If the same integral type is used for args and return val, -// make sure the arguments use only the half of the type's bit depth. -template> -IntegerOnly pairhash(I a, I b) -{ - using std::ceil; using std::log2; using std::max; using std::min; - static const auto constexpr Ibits = int(sizeof(I) * CHAR_BIT); - static const auto constexpr DoubleIbits = int(sizeof(DoubleI) * CHAR_BIT); - static const auto constexpr shift = DoubleIbits / 2 < Ibits ? Ibits / 2 : Ibits; - - I g = min(a, b), l = max(a, b); - - // Assume the hash will fit into the output variable - assert((g ? (ceil(log2(g))) : 0) <= shift); - assert((l ? (ceil(log2(l))) : 0) <= shift); - - return (DoubleI(g) << shift) + l; -} - -class SupportTreeBuildsteps { - const SupportTreeConfig& m_cfg; - const IndexedMesh& m_mesh; - const std::vector& m_support_pts; +class DefaultSupportTree { + const SupportableMesh &m_sm; using PtIndices = std::vector; PtIndices m_iheads; // support points with pinhead PtIndices m_iheads_onmodel; PtIndices m_iheadless; // headless support points - - std::map m_head_to_ground_scans; + + std::map m_head_to_ground_scans; // normals for support points from model faces. - PointSet m_support_nmls; + Eigen::MatrixXd m_support_nmls; // Clusters of points which can reach the ground directly and can be // bridged to one central pillar @@ -210,7 +78,7 @@ class SupportTreeBuildsteps { SupportTreeBuilder& m_builder; // support points in Eigen/IGL format - PointSet m_points; + Eigen::MatrixXd m_points; // throw if canceled: It will be called many times so a shorthand will // come in handy. @@ -220,12 +88,12 @@ class SupportTreeBuildsteps { PillarIndex m_pillar_index; // When bridging heads to pillars... TODO: find a cleaner solution - ccr::BlockingMutex m_bridge_mutex; + execution::BlockingMutex m_bridge_mutex; - inline IndexedMesh::hit_result ray_mesh_intersect(const Vec3d& s, + inline AABBMesh::hit_result ray_mesh_intersect(const Vec3d& s, const Vec3d& dir) { - return m_mesh.query_ray_hit(s, dir); + return m_sm.emesh.query_ray_hit(s, dir); } // This function will test if a future pinhead would not collide with the @@ -239,7 +107,7 @@ class SupportTreeBuildsteps { // point was inside the model, an "invalid" hit_result will be returned // with a zero distance value instead of a NAN. This way the result can // be used safely for comparison with other distances. - IndexedMesh::hit_result pinhead_mesh_intersect( + AABBMesh::hit_result pinhead_mesh_intersect( const Vec3d& s, const Vec3d& dir, double r_pin, @@ -247,16 +115,15 @@ class SupportTreeBuildsteps { double width, double safety_d); - IndexedMesh::hit_result pinhead_mesh_intersect( - const Vec3d& s, - const Vec3d& dir, - double r_pin, - double r_back, - double width) + AABBMesh::hit_result pinhead_mesh_intersect(const Vec3d &s, + const Vec3d &dir, + double r_pin, + double r_back, + double width) { return pinhead_mesh_intersect(s, dir, r_pin, r_back, width, - r_back * m_cfg.safety_distance_mm / - m_cfg.head_back_radius_mm); + r_back * m_sm.cfg.safety_distance_mm / + m_sm.cfg.head_back_radius_mm); } // Checking bridge (pillar and stick as well) intersection with the model. @@ -267,22 +134,22 @@ class SupportTreeBuildsteps { // point was inside the model, an "invalid" hit_result will be returned // with a zero distance value instead of a NAN. This way the result can // be used safely for comparison with other distances. - IndexedMesh::hit_result bridge_mesh_intersect( + AABBMesh::hit_result bridge_mesh_intersect( const Vec3d& s, const Vec3d& dir, double r, double safety_d); - IndexedMesh::hit_result bridge_mesh_intersect( + AABBMesh::hit_result bridge_mesh_intersect( const Vec3d& s, const Vec3d& dir, double r) { return bridge_mesh_intersect(s, dir, r, - r * m_cfg.safety_distance_mm / - m_cfg.head_back_radius_mm); + r * m_sm.cfg.safety_distance_mm / + m_sm.cfg.head_back_radius_mm); } - + template inline double bridge_mesh_distance(Args&&...args) { return bridge_mesh_intersect(std::forward(args)...).distance(); @@ -293,18 +160,17 @@ class SupportTreeBuildsteps { // For connecting a head to a nearby pillar. bool connect_to_nearpillar(const Head& head, long nearpillar_id); - + // Find route for a head to the ground. Inserts additional bridge from the // head to the pillar if cannot create pillar directly. // The optional dir parameter is the direction of the bridge which is the // direction of the pinhead if omitted. - bool connect_to_ground(Head& head, const Vec3d &dir); inline bool connect_to_ground(Head& head); - + bool connect_to_model_body(Head &head); bool search_pillar_and_connect(const Head& source); - + // This is a proxy function for pillar creation which will mind the gap // between the pad and the model bottom in zero elevation mode. // jp is the starting junction point which needs to be routed down. @@ -317,16 +183,19 @@ class SupportTreeBuildsteps { void add_pillar_base(long pid) { - m_builder.add_pillar_base(pid, m_cfg.base_height_mm, m_cfg.base_radius_mm); + m_builder.add_pillar_base(pid, m_sm.cfg.base_height_mm, m_sm.cfg.base_radius_mm); } std::optional search_widening_path(const Vec3d &jp, const Vec3d &dir, double radius, - double new_radius); + double new_radius) + { + return sla::search_widening_path(suptree_ex_policy, m_sm, jp, dir, radius, new_radius); + } public: - SupportTreeBuildsteps(SupportTreeBuilder & builder, const SupportableMesh &sm); + DefaultSupportTree(SupportTreeBuilder & builder, const SupportableMesh &sm); // Now let's define the individual steps of the support generation algorithm @@ -335,10 +204,6 @@ public: // pinhead is applicable and adjust its angle at each support point. We // will also merge the support points that are just too close and can // be considered as one. - void filter(); - - // Pinhead creation: based on the filtering results, the Head objects - // will be constructed (together with their triangle meshes). void add_pinheads(); // Further classification of the support points with pinheads. If the @@ -372,7 +237,11 @@ public: static bool execute(SupportTreeBuilder & builder, const SupportableMesh &sm); }; -} +inline void create_default_tree(SupportTreeBuilder &builder, const SupportableMesh &sm) +{ + DefaultSupportTree::execute(builder, sm); } -#endif // SLASUPPORTTREEALGORITHM_H +}} // namespace Slic3r::sla + +#endif // LEGACYSUPPORTTREE_HPP diff --git a/src/libslic3r/SLA/Hollowing.cpp b/src/libslic3r/SLA/Hollowing.cpp index 8cda341b5..74eb07695 100644 --- a/src/libslic3r/SLA/Hollowing.cpp +++ b/src/libslic3r/SLA/Hollowing.cpp @@ -5,10 +5,11 @@ #include #include #include -#include +#include #include #include #include +#include #include @@ -219,7 +220,7 @@ bool DrainHole::get_intersections(const Vec3f& s, const Vec3f& dir, const Eigen::ParametrizedLine ray(s, dir.normalized()); for (size_t i=0; i<2; ++i) - out[i] = std::make_pair(sla::IndexedMesh::hit_result::infty(), Vec3d::Zero()); + out[i] = std::make_pair(AABBMesh::hit_result::infty(), Vec3d::Zero()); const float sqr_radius = pow(radius, 2.f); @@ -448,7 +449,7 @@ void remove_inside_triangles(TriangleMesh &mesh, const Interior &interior, }; // TODO: Parallel mode not working yet - using exec_policy = ccr_seq; + constexpr auto &exec_policy = ex_seq; // Info about the needed modifications on the input mesh. struct MeshMods { @@ -456,7 +457,7 @@ void remove_inside_triangles(TriangleMesh &mesh, const Interior &interior, // Just a thread safe wrapper for a vector of triangles. struct { std::vector> data; - exec_policy::SpinningMutex mutex; + execution::SpinningMutex mutex; void emplace_back(const std::array &pts) { @@ -533,27 +534,28 @@ void remove_inside_triangles(TriangleMesh &mesh, const Interior &interior, interior.reset_accessor(); - exec_policy::for_each(size_t(0), faces.size(), [&] (size_t face_idx) { - const Vec3i &face = faces[face_idx]; + execution::for_each( + exec_policy, size_t(0), faces.size(), + [&](size_t face_idx) { + const Vec3i &face = faces[face_idx]; - // If the triangle is excluded, we need to keep it. - if (is_excluded(face_idx)) - return; + // If the triangle is excluded, we need to keep it. + if (is_excluded(face_idx)) return; - std::array pts = - { vertices[face(0)], vertices[face(1)], vertices[face(2)] }; + std::array pts = {vertices[face(0)], vertices[face(1)], + vertices[face(2)]}; - BoundingBoxf3 facebb { pts.begin(), pts.end() }; + BoundingBoxf3 facebb{pts.begin(), pts.end()}; - // Face is certainly outside the cavity - if (! facebb.intersects(bb)) return; + // Face is certainly outside the cavity + if (!facebb.intersects(bb)) return; - DivFace df{face, pts, long(face_idx)}; + DivFace df{face, pts, long(face_idx)}; - if (divfn(df)) - divide_triangle(df, divfn); - - }, exec_policy::max_concurreny()); + if (divfn(df)) divide_triangle(df, divfn); + }, + execution::max_concurrency(exec_policy) + ); auto new_faces = reserve_vector(faces.size() + mesh_mods.new_triangles.size()); diff --git a/src/libslic3r/SLA/RasterBase.hpp b/src/libslic3r/SLA/RasterBase.hpp index 657fc865c..ae8c0aefa 100644 --- a/src/libslic3r/SLA/RasterBase.hpp +++ b/src/libslic3r/SLA/RasterBase.hpp @@ -9,7 +9,6 @@ #include #include -//#include namespace Slic3r { diff --git a/src/libslic3r/SLA/ReprojectPointsOnMesh.hpp b/src/libslic3r/SLA/ReprojectPointsOnMesh.hpp index 3ad7d62b1..fc033b832 100644 --- a/src/libslic3r/SLA/ReprojectPointsOnMesh.hpp +++ b/src/libslic3r/SLA/ReprojectPointsOnMesh.hpp @@ -2,9 +2,9 @@ #define REPROJECTPOINTSONMESH_HPP #include "libslic3r/Point.hpp" +#include "libslic3r/AABBMesh.hpp" #include "SupportPoint.hpp" #include "Hollowing.hpp" -#include "IndexedMesh.hpp" #include "libslic3r/Model.hpp" #include @@ -15,7 +15,7 @@ template Vec3d pos(const Pt &p) { return p.pos.template cast() template void pos(Pt &p, const Vec3d &pp) { p.pos = pp.cast(); } template -void reproject_support_points(const IndexedMesh &mesh, std::vector &pts) +void reproject_support_points(const AABBMesh &mesh, std::vector &pts) { tbb::parallel_for(size_t(0), pts.size(), [&mesh, &pts](size_t idx) { int junk; @@ -33,7 +33,7 @@ inline void reproject_points_and_holes(ModelObject *object) if (!object || (!has_holes && !has_sppoints)) return; TriangleMesh rmsh = object->raw_mesh(); - IndexedMesh emesh{rmsh}; + AABBMesh emesh{rmsh}; if (has_sppoints) reproject_support_points(emesh, object->sla_support_points); diff --git a/src/libslic3r/SLA/SupportPointGenerator.cpp b/src/libslic3r/SLA/SupportPointGenerator.cpp index 4c1af03eb..6d159b37b 100644 --- a/src/libslic3r/SLA/SupportPointGenerator.cpp +++ b/src/libslic3r/SLA/SupportPointGenerator.cpp @@ -1,7 +1,7 @@ #include #include "SupportPointGenerator.hpp" -#include "Concurrency.hpp" +#include "Execution/ExecutionTBB.hpp" #include "Geometry/ConvexHull.hpp" #include "Model.hpp" #include "ExPolygon.hpp" @@ -50,7 +50,7 @@ float SupportPointGenerator::distance_limit(float angle) const }*/ SupportPointGenerator::SupportPointGenerator( - const sla::IndexedMesh &emesh, + const AABBMesh &emesh, const std::vector &slices, const std::vector & heights, const Config & config, @@ -64,7 +64,7 @@ SupportPointGenerator::SupportPointGenerator( } SupportPointGenerator::SupportPointGenerator( - const IndexedMesh &emesh, + const AABBMesh &emesh, const SupportPointGenerator::Config &config, std::function throw_on_cancel, std::function statusfn) @@ -89,7 +89,7 @@ void SupportPointGenerator::project_onto_mesh(std::vector& po // Use a reasonable granularity to account for the worker thread synchronization cost. static constexpr size_t gransize = 64; - ccr_par::for_each(size_t(0), points.size(), [this, &points](size_t idx) + execution::for_each(ex_tbb, size_t(0), points.size(), [this, &points](size_t idx) { if ((idx % 16) == 0) // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves. @@ -97,8 +97,8 @@ void SupportPointGenerator::project_onto_mesh(std::vector& po Vec3f& p = points[idx].pos; // Project the point upward and downward and choose the closer intersection with the mesh. - sla::IndexedMesh::hit_result hit_up = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., 1.)); - sla::IndexedMesh::hit_result hit_down = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., -1.)); + AABBMesh::hit_result hit_up = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., 1.)); + AABBMesh::hit_result hit_down = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., -1.)); bool up = hit_up.is_hit(); bool down = hit_down.is_hit(); @@ -106,7 +106,7 @@ void SupportPointGenerator::project_onto_mesh(std::vector& po if (!up && !down) return; - sla::IndexedMesh::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down; + AABBMesh::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down; p = p + (hit.distance() * hit.direction()).cast(); }, gransize); } @@ -127,7 +127,7 @@ static std::vector make_layers( //const float pixel_area = pow(wxGetApp().preset_bundle->project_config.option("display_width") / wxGetApp().preset_bundle->project_config.option("display_pixels_x"), 2.f); // const float pixel_area = pow(0.047f, 2.f); - ccr_par::for_each(size_t(0), layers.size(), + execution::for_each(ex_tbb, size_t(0), layers.size(), [&layers, &slices, &heights, pixel_area, throw_on_cancel](size_t layer_id) { if ((layer_id % 8) == 0) @@ -152,7 +152,7 @@ static std::vector make_layers( }, 32 /*gransize*/); // Calculate overlap of successive layers. Link overlapping islands. - ccr_par::for_each(size_t(1), layers.size(), + execution::for_each(ex_tbb, size_t(1), layers.size(), [&layers, &heights, throw_on_cancel] (size_t layer_id) { if ((layer_id % 2) == 0) diff --git a/src/libslic3r/SLA/SupportPointGenerator.hpp b/src/libslic3r/SLA/SupportPointGenerator.hpp index 441b82de1..5cb7b7069 100644 --- a/src/libslic3r/SLA/SupportPointGenerator.hpp +++ b/src/libslic3r/SLA/SupportPointGenerator.hpp @@ -3,9 +3,9 @@ #include -#include -#include +#include +#include #include #include #include @@ -28,10 +28,10 @@ public: inline float tear_pressure() const { return 1.f; } // pressure that the display exerts (the force unit per mm2) }; - SupportPointGenerator(const IndexedMesh& emesh, const std::vector& slices, + SupportPointGenerator(const AABBMesh& emesh, const std::vector& slices, const std::vector& heights, const Config& config, std::function throw_on_cancel, std::function statusfn); - SupportPointGenerator(const IndexedMesh& emesh, const Config& config, std::function throw_on_cancel, std::function statusfn); + SupportPointGenerator(const AABBMesh& emesh, const Config& config, std::function throw_on_cancel, std::function statusfn); const std::vector& output() const { return m_output; } std::vector& output() { return m_output; } @@ -217,7 +217,7 @@ private: static void output_structures(const std::vector &structures); #endif // SLA_SUPPORTPOINTGEN_DEBUG - const IndexedMesh& m_emesh; + const AABBMesh& m_emesh; std::function m_throw_on_cancel; std::function m_statusfn; diff --git a/src/libslic3r/SLA/SupportTree.cpp b/src/libslic3r/SLA/SupportTree.cpp index a0f88ed60..466fade88 100644 --- a/src/libslic3r/SLA/SupportTree.cpp +++ b/src/libslic3r/SLA/SupportTree.cpp @@ -7,15 +7,13 @@ #include #include #include -#include +#include #include #include #include #include -#include -#include #include #include @@ -23,43 +21,96 @@ //! return same string #define L(s) Slic3r::I18N::translate(s) -namespace Slic3r { -namespace sla { +namespace Slic3r { namespace sla { -void SupportTree::retrieve_full_mesh(indexed_triangle_set &outmesh) const { - its_merge(outmesh, retrieve_mesh(MeshType::Support)); - its_merge(outmesh, retrieve_mesh(MeshType::Pad)); +indexed_triangle_set create_support_tree(const SupportableMesh &sm, + const JobController &ctl) +{ + auto builder = make_unique(ctl); + + if (sm.cfg.enabled) { + switch (sm.cfg.tree_type) { + case SupportTreeType::Default: { + create_default_tree(*builder, sm); + break; + } + default:; + } + + builder->merge_and_cleanup(); // clean metadata, leave only the meshes. + } + + indexed_triangle_set out = builder->retrieve_mesh(MeshType::Support); + + return out; } -std::vector SupportTree::slice(const std::vector &grid, - float cr) const +indexed_triangle_set create_pad(const SupportableMesh &sm, + const indexed_triangle_set &support_mesh, + const JobController &ctl) { - const indexed_triangle_set &sup_mesh = retrieve_mesh(MeshType::Support); - const indexed_triangle_set &pad_mesh = retrieve_mesh(MeshType::Pad); + ExPolygons model_contours; // This will store the base plate of the pad. + double pad_h = sm.pad_cfg.full_height(); + + float zstart = ground_level(sm); + float zend = zstart + float(pad_h + EPSILON); + auto heights = grid(zstart, zend, 0.1f); + + if (!sm.cfg.enabled || sm.pad_cfg.embed_object) { + // No support (thus no elevation) or zero elevation mode + // we sometimes call it "builtin pad" is enabled so we will + // get a sample from the bottom of the mesh and use it for pad + // creation. + sla::pad_blueprint(*sm.emesh.get_triangle_mesh(), model_contours, + heights, ctl.cancelfn); + } + + ExPolygons sup_contours; + pad_blueprint(support_mesh, sup_contours, heights, ctl.cancelfn); + + indexed_triangle_set out; + create_pad(sup_contours, model_contours, out, sm.pad_cfg); + + Vec3f offs{.0f, .0f, zstart}; + for (auto &p : out.vertices) p += offs; + + its_merge_vertices(out); + + return out; +} + +std::vector slice(const indexed_triangle_set &sup_mesh, + const indexed_triangle_set &pad_mesh, + const std::vector &grid, + float cr, + const JobController &ctl) +{ using Slices = std::vector; + auto slices = reserve_vector(2); if (!sup_mesh.empty()) { slices.emplace_back(); - slices.back() = slice_mesh_ex(sup_mesh, grid, cr, ctl().cancelfn); + slices.back() = slice_mesh_ex(sup_mesh, grid, cr, ctl.cancelfn); } if (!pad_mesh.empty()) { slices.emplace_back(); - auto bb = bounding_box(pad_mesh); + auto bb = bounding_box(pad_mesh); auto maxzit = std::upper_bound(grid.begin(), grid.end(), bb.max.z()); - - auto cap = grid.end() - maxzit; + + auto cap = grid.end() - maxzit; auto padgrid = reserve_vector(size_t(cap > 0 ? cap : 0)); std::copy(grid.begin(), maxzit, std::back_inserter(padgrid)); - slices.back() = slice_mesh_ex(pad_mesh, padgrid, cr, ctl().cancelfn); + slices.back() = slice_mesh_ex(pad_mesh, padgrid, cr, ctl.cancelfn); } size_t len = grid.size(); - for (const Slices &slv : slices) { len = std::min(len, slv.size()); } + for (const Slices &slv : slices) + len = std::min(len, slv.size()); // Either the support or the pad or both has to be non empty if (slices.empty()) return {}; @@ -77,22 +128,4 @@ std::vector SupportTree::slice(const std::vector &grid, return mrg; } -SupportTree::UPtr SupportTree::create(const SupportableMesh &sm, - const JobController & ctl) -{ - auto builder = make_unique(); - builder->m_ctl = ctl; - - if (sm.cfg.enabled) { - // Execute takes care about the ground_level - SupportTreeBuildsteps::execute(*builder, sm); - builder->merge_and_cleanup(); // clean metadata, leave only the meshes. - } else { - // If a pad gets added later, it will be in the right Z level - builder->ground_level = sm.emesh.ground_level(); - } - - return std::move(builder); -} - }} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTree.hpp b/src/libslic3r/SLA/SupportTree.hpp index 200b7cbde..0ac29ff7a 100644 --- a/src/libslic3r/SLA/SupportTree.hpp +++ b/src/libslic3r/SLA/SupportTree.hpp @@ -3,38 +3,29 @@ #include #include -#include + +#include +#include #include -#include #include #include +#include namespace Slic3r { -class TriangleMesh; -class Model; -class ModelInstance; -class ModelObject; -class Polygon; -class ExPolygon; - using Polygons = std::vector; using ExPolygons = std::vector; namespace sla { -enum class PillarConnectionMode -{ - zigzag, - cross, - dynamic -}; - struct SupportTreeConfig { bool enabled = true; - + + // Type of the support tree, for + SupportTreeType tree_type = SupportTreeType::Default; + // Radius in mm of the pointing side of the head. double head_front_radius_mm = 0.2; @@ -59,7 +50,7 @@ struct SupportTreeConfig // when bridges and pillars are merged. The resulting pillar should be a bit // thicker than the ones merging into it. How much thicker? I don't know // but it will be derived from this value. - double pillar_widening_factor = 0.5; + double pillar_widening_factor = .05; // Radius in mm of the pillar base. double base_radius_mm = 2.0; @@ -109,68 +100,51 @@ struct SupportTreeConfig }; -// TODO: Part of future refactor -//class SupportConfig { -// std::optional tree_cfg {std::in_place_t{}}; // fill up -// std::optional pad_cfg; -//}; - enum class MeshType { Support, Pad }; struct SupportableMesh { - IndexedMesh emesh; - SupportPoints pts; + AABBMesh emesh; + SupportPoints pts; SupportTreeConfig cfg; -// PadConfig pad_cfg; + PadConfig pad_cfg; - explicit SupportableMesh(const indexed_triangle_set & trmsh, - const SupportPoints &sp, - const SupportTreeConfig &c) + explicit SupportableMesh(const indexed_triangle_set &trmsh, + const SupportPoints &sp, + const SupportTreeConfig &c) : emesh{trmsh}, pts{sp}, cfg{c} {} - - explicit SupportableMesh(const IndexedMesh &em, - const SupportPoints &sp, + + explicit SupportableMesh(const AABBMesh &em, + const SupportPoints &sp, const SupportTreeConfig &c) : emesh{em}, pts{sp}, cfg{c} {} }; -/// The class containing mesh data for the generated supports. -class SupportTree +inline double ground_level(const SupportableMesh &sm) { - JobController m_ctl; -public: - using UPtr = std::unique_ptr; - - static UPtr create(const SupportableMesh &input, - const JobController &ctl = {}); - - virtual ~SupportTree() = default; - - virtual const indexed_triangle_set &retrieve_mesh(MeshType meshtype) const = 0; - - /// Adding the "pad" under the supports. - /// modelbase will be used according to the embed_object flag in PoolConfig. - /// If set, the plate will be interpreted as the model's intrinsic pad. - /// Otherwise, the modelbase will be unified with the base plate calculated - /// from the supports. - virtual const indexed_triangle_set &add_pad(const ExPolygons &modelbase, - const PadConfig & pcfg) = 0; - - virtual void remove_pad() = 0; - - std::vector slice(const std::vector &, - float closing_radius) const; - - void retrieve_full_mesh(indexed_triangle_set &outmesh) const; - - const JobController &ctl() const { return m_ctl; } -}; + double lvl = sm.emesh.ground_level() - + !bool(sm.pad_cfg.embed_object) * sm.cfg.enabled * sm.cfg.object_elevation_mm + + bool(sm.pad_cfg.embed_object) * sm.pad_cfg.wall_thickness_mm; + return lvl; } -} +indexed_triangle_set create_support_tree(const SupportableMesh &mesh, + const JobController &ctl); + +indexed_triangle_set create_pad(const SupportableMesh &model_mesh, + const indexed_triangle_set &support_mesh, + const JobController &ctl); + +std::vector slice(const indexed_triangle_set &support_mesh, + const indexed_triangle_set &pad_mesh, + const std::vector &grid, + float closing_radius, + const JobController &ctl); + +} // namespace sla +} // namespace Slic3r #endif // SLASUPPORTTREE_HPP diff --git a/src/libslic3r/SLA/SupportTreeBuilder.cpp b/src/libslic3r/SLA/SupportTreeBuilder.cpp index 86339d2ac..1786adb83 100644 --- a/src/libslic3r/SLA/SupportTreeBuilder.cpp +++ b/src/libslic3r/SLA/SupportTreeBuilder.cpp @@ -1,7 +1,7 @@ #define NOMINMAX #include -#include +#include #include //#include @@ -23,48 +23,15 @@ Head::Head(double r_big_mm, { } -Pad::Pad(const indexed_triangle_set &support_mesh, - const ExPolygons & model_contours, - double ground_level, - const PadConfig & pcfg, - ThrowOnCancel thr) - : cfg(pcfg) - , zlevel(ground_level + pcfg.full_height() - pcfg.required_elevation()) -{ - thr(); - - ExPolygons sup_contours; - - float zstart = float(zlevel); - float zend = zstart + float(pcfg.full_height() + EPSILON); - - pad_blueprint(support_mesh, sup_contours, grid(zstart, zend, 0.1f), thr); - create_pad(sup_contours, model_contours, tmesh, pcfg); - - Vec3f offs{.0f, .0f, float(zlevel)}; - for (auto &p : tmesh.vertices) p += offs; - - its_merge_vertices(tmesh); -} - -const indexed_triangle_set &SupportTreeBuilder::add_pad( - const ExPolygons &modelbase, const PadConfig &cfg) -{ - m_pad = Pad{merged_mesh(), modelbase, ground_level, cfg, ctl().cancelfn}; - return m_pad.tmesh; -} - SupportTreeBuilder::SupportTreeBuilder(SupportTreeBuilder &&o) : m_heads(std::move(o.m_heads)) , m_head_indices{std::move(o.m_head_indices)} , m_pillars{std::move(o.m_pillars)} , m_bridges{std::move(o.m_bridges)} , m_crossbridges{std::move(o.m_crossbridges)} - , m_pad{std::move(o.m_pad)} , m_meshcache{std::move(o.m_meshcache)} , m_meshcache_valid{o.m_meshcache_valid} , m_model_height{o.m_model_height} - , ground_level{o.ground_level} {} SupportTreeBuilder::SupportTreeBuilder(const SupportTreeBuilder &o) @@ -73,11 +40,9 @@ SupportTreeBuilder::SupportTreeBuilder(const SupportTreeBuilder &o) , m_pillars{o.m_pillars} , m_bridges{o.m_bridges} , m_crossbridges{o.m_crossbridges} - , m_pad{o.m_pad} , m_meshcache{o.m_meshcache} , m_meshcache_valid{o.m_meshcache_valid} , m_model_height{o.m_model_height} - , ground_level{o.ground_level} {} SupportTreeBuilder &SupportTreeBuilder::operator=(SupportTreeBuilder &&o) @@ -87,11 +52,9 @@ SupportTreeBuilder &SupportTreeBuilder::operator=(SupportTreeBuilder &&o) m_pillars = std::move(o.m_pillars); m_bridges = std::move(o.m_bridges); m_crossbridges = std::move(o.m_crossbridges); - m_pad = std::move(o.m_pad); m_meshcache = std::move(o.m_meshcache); m_meshcache_valid = o.m_meshcache_valid; m_model_height = o.m_model_height; - ground_level = o.ground_level; return *this; } @@ -102,11 +65,9 @@ SupportTreeBuilder &SupportTreeBuilder::operator=(const SupportTreeBuilder &o) m_pillars = o.m_pillars; m_bridges = o.m_bridges; m_crossbridges = o.m_crossbridges; - m_pad = o.m_pad; m_meshcache = o.m_meshcache; m_meshcache_valid = o.m_meshcache_valid; m_model_height = o.m_model_height; - ground_level = o.ground_level; return *this; } @@ -116,7 +77,7 @@ void SupportTreeBuilder::add_pillar_base(long pid, double baseheight, double rad assert(pid >= 0 && size_t(pid) < m_pillars.size()); Pillar& pll = m_pillars[size_t(pid)]; m_pedestals.emplace_back(pll.endpt, std::min(baseheight, pll.height), - std::max(radius, pll.r), pll.r); + std::max(radius, pll.r_start), pll.r_start); m_pedestals.back().id = m_pedestals.size() - 1; m_meshcache_valid = false; @@ -174,7 +135,7 @@ const indexed_triangle_set &SupportTreeBuilder::merged_mesh(size_t steps) const return m_meshcache; } - m_meshcache = merged; + m_meshcache = std::move(merged); // The mesh will be passed by const-pointer to TriangleMeshSlicer, // which will need this. @@ -187,36 +148,28 @@ const indexed_triangle_set &SupportTreeBuilder::merged_mesh(size_t steps) const return m_meshcache; } -double SupportTreeBuilder::full_height() const -{ - if (merged_mesh().indices.empty() && !pad().empty()) - return pad().cfg.full_height(); - - double h = mesh_height(); - if (!pad().empty()) h += pad().cfg.required_elevation(); - return h; -} - const indexed_triangle_set &SupportTreeBuilder::merge_and_cleanup() { // in case the mesh is not generated, it should be... auto &ret = merged_mesh(); // Doing clear() does not garantee to release the memory. - m_heads = {}; - m_head_indices = {}; - m_pillars = {}; - m_junctions = {}; - m_bridges = {}; - + clear_and_shrink(m_heads); + clear_and_shrink(m_head_indices); + clear_and_shrink(m_pillars); + clear_and_shrink(m_junctions); + clear_and_shrink(m_bridges); + return ret; } const indexed_triangle_set &SupportTreeBuilder::retrieve_mesh(MeshType meshtype) const { + static const indexed_triangle_set EMPTY_MESH; + switch(meshtype) { case MeshType::Support: return merged_mesh(); - case MeshType::Pad: return pad().tmesh; + case MeshType::Pad: return EMPTY_MESH; //pad().tmesh; } return m_meshcache; diff --git a/src/libslic3r/SLA/SupportTreeBuilder.hpp b/src/libslic3r/SLA/SupportTreeBuilder.hpp index df5484e3d..29d34ab8e 100644 --- a/src/libslic3r/SLA/SupportTreeBuilder.hpp +++ b/src/libslic3r/SLA/SupportTreeBuilder.hpp @@ -1,9 +1,8 @@ #ifndef SLA_SUPPORTTREEBUILDER_HPP #define SLA_SUPPORTTREEBUILDER_HPP -#include +#include #include -//#include #include #include #include @@ -50,13 +49,13 @@ namespace sla { * nearby pillar. */ -template double distance(const Vec& p) { - return std::sqrt(p.transpose() * p); +template T distance(const Vec& p) { + return p.norm(); } -template double distance(const Vec& pp1, const Vec& pp2) { - auto p = pp2 - pp1; - return distance(p); +template +T distance(const Vec& pp1, const Vec& pp2) { + return (pp1 - pp2).norm(); } const Vec3d DOWN = {0.0, 0.0, -1.0}; @@ -109,12 +108,6 @@ struct Head: public SupportTreeNode { { return pos + (fullwidth() - r_back_mm) * dir; } - - inline double request_pillar_radius(double radius) const - { - const double rmax = r_back_mm; - return radius > 0 && radius < rmax ? radius : rmax; - } }; // A junction connecting bridges and pillars @@ -125,8 +118,11 @@ struct Junction: public SupportTreeNode { Junction(const Vec3d &tr, double r_mm) : r(r_mm), pos(tr) {} }; +// A straight pillar. Only has an endpoint and a height. No explicit starting +// point is given, as it would allow the pillar to be angled. +// Some connection info with other primitives can also be tracked. struct Pillar: public SupportTreeNode { - double height, r; + double height, r_start, r_end; Vec3d endpt; // If the pillar connects to a head, this is the id of that head @@ -139,8 +135,17 @@ struct Pillar: public SupportTreeNode { // How many pillars are cascaded with this one unsigned links = 0; - Pillar(const Vec3d &endp, double h, double radius = 1.): - height{h}, r(radius), endpt(endp), starts_from_head(false) {} + Pillar(const Vec3d &endp, double h, double start_radius, double end_radius) + : height{h} + , r_start(start_radius) + , r_end(end_radius) + , endpt(endp) + , starts_from_head(false) + {} + + Pillar(const Vec3d &endp, double h, double start_radius) + : Pillar(endp, h, start_radius, start_radius) + {} Vec3d startpoint() const { @@ -186,38 +191,19 @@ struct DiffBridge: public Bridge { {} }; -// A wrapper struct around the pad -struct Pad { - indexed_triangle_set tmesh; - PadConfig cfg; - double zlevel = 0; - - Pad() = default; - - Pad(const indexed_triangle_set &support_mesh, - const ExPolygons & model_contours, - double ground_level, - const PadConfig & pcfg, - ThrowOnCancel thr); - - bool empty() const { return tmesh.indices.size() == 0; } -}; - -// This class will hold the support tree meshes with some additional -// bookkeeping as well. Various parts of the support geometry are stored -// separately and are merged when the caller queries the merged mesh. The -// merged result is cached for fast subsequent delivery of the merged mesh -// which can be quite complex. The support tree creation algorithm can use an -// instance of this class as a somewhat higher level tool for crafting the 3D -// support mesh. Parts can be added with the appropriate methods such as -// add_head or add_pillar which forwards the constructor arguments and fills -// the IDs of these substructures. The IDs are basically indices into the -// arrays of the appropriate type (heads, pillars, etc...). One can later query -// e.g. a pillar for a specific head... -// -// The support pad is considered an auxiliary geometry and is not part of the -// merged mesh. It can be retrieved using a dedicated method (pad()) -class SupportTreeBuilder: public SupportTree { +// This class will hold the support tree parts (not meshes, but logical parts) +// with some additional bookkeeping as well. Various parts of the support +// geometry are stored separately and are merged when the caller queries the +// merged mesh (for every part, there is a meshing routine, see +// SupportTreeMesher.hpp). The merged result is cached for fast subsequent +// delivery of the merged mesh which can be quite complex. The support tree +// creation algorithm can use an instance of this class as a somewhat higher +// level tool for crafting the 3D support mesh. Parts can be added with the +// appropriate methods such as add_head or add_pillar which forwards the +// constructor arguments and fills the IDs of these substructures. The IDs are +// basically indices into the arrays of the appropriate type (heads, pillars, +// etc...). One can later query e.g. a pillar for a specific head... +class SupportTreeBuilder { // For heads it is beneficial to use the same IDs as for the support points. std::vector m_heads; std::vector m_head_indices; @@ -229,9 +215,9 @@ class SupportTreeBuilder: public SupportTree { std::vector m_pedestals; std::vector m_anchors; - Pad m_pad; + JobController m_ctl; - using Mutex = ccr::SpinningMutex; + using Mutex = tbb::spin_mutex; mutable indexed_triangle_set m_meshcache; mutable Mutex m_mutex; @@ -249,14 +235,15 @@ class SupportTreeBuilder: public SupportTree { } public: - double ground_level = 0; - SupportTreeBuilder() = default; + explicit SupportTreeBuilder(const JobController &ctl = {}) : m_ctl{ctl} {} SupportTreeBuilder(SupportTreeBuilder &&o); SupportTreeBuilder(const SupportTreeBuilder &o); SupportTreeBuilder& operator=(SupportTreeBuilder &&o); SupportTreeBuilder& operator=(const SupportTreeBuilder &o); + const JobController &ctl() const { return m_ctl; } + template Head& add_head(unsigned id, Args&&... args) { std::lock_guard lk(m_mutex); @@ -270,7 +257,7 @@ public: return m_heads.back(); } - template long add_pillar(long headid, double length) + long add_pillar(long headid, double length) { std::lock_guard lk(m_mutex); if (m_pillars.capacity() < m_heads.size()) @@ -415,34 +402,21 @@ public: return m_pillars[size_t(id)]; } - - const Pad& pad() const { return m_pad; } - + // WITHOUT THE PAD!!! const indexed_triangle_set &merged_mesh(size_t steps = 45) const; - // WITH THE PAD - double full_height() const; - - // WITHOUT THE PAD!!! - inline double mesh_height() const - { - if (!m_meshcache_valid) merged_mesh(); - return m_model_height; - } - // Intended to be called after the generation is fully complete const indexed_triangle_set & merge_and_cleanup(); - - // Implement SupportTree interface: - const indexed_triangle_set &add_pad(const ExPolygons &modelbase, - const PadConfig & pcfg) override; + const indexed_triangle_set &retrieve_mesh( + MeshType meshtype = MeshType::Support) const; - void remove_pad() override { m_pad = Pad(); } - - virtual const indexed_triangle_set &retrieve_mesh( - MeshType meshtype = MeshType::Support) const override; + void retrieve_full_mesh(indexed_triangle_set &outmesh) const + { + its_merge(outmesh, retrieve_mesh(MeshType::Support)); + its_merge(outmesh, retrieve_mesh(MeshType::Pad)); + } }; }} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp deleted file mode 100644 index e4ca3f59c..000000000 --- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp +++ /dev/null @@ -1,1277 +0,0 @@ -#include - -#include -#include -#include - -namespace Slic3r { -namespace sla { - -using Slic3r::opt::initvals; -using Slic3r::opt::bounds; -using Slic3r::opt::StopCriteria; -using Slic3r::opt::Optimizer; -using Slic3r::opt::AlgNLoptSubplex; -using Slic3r::opt::AlgNLoptGenetic; - -StopCriteria get_criteria(const SupportTreeConfig &cfg) -{ - return StopCriteria{} - .rel_score_diff(cfg.optimizer_rel_score_diff) - .max_iterations(cfg.optimizer_max_iterations); -} - -template -static Hit min_hit(const C &hits) -{ - auto mit = std::min_element(hits.begin(), hits.end(), - [](const Hit &h1, const Hit &h2) { - return h1.distance() < h2.distance(); - }); - - return *mit; -} - -SupportTreeBuildsteps::SupportTreeBuildsteps(SupportTreeBuilder & builder, - const SupportableMesh &sm) - : m_cfg(sm.cfg) - , m_mesh(sm.emesh) - , m_support_pts(sm.pts) - , m_support_nmls(sm.pts.size(), 3) - , m_builder(builder) - , m_points(sm.pts.size(), 3) - , m_thr(builder.ctl().cancelfn) -{ - // Prepare the support points in Eigen/IGL format as well, we will use - // it mostly in this form. - - long i = 0; - for (const SupportPoint &sp : m_support_pts) { - m_points.row(i)(X) = double(sp.pos(X)); - m_points.row(i)(Y) = double(sp.pos(Y)); - m_points.row(i)(Z) = double(sp.pos(Z)); - ++i; - } -} - -bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, - const SupportableMesh &sm) -{ - if(sm.pts.empty()) return false; - - builder.ground_level = sm.emesh.ground_level() - sm.cfg.object_elevation_mm; - - SupportTreeBuildsteps alg(builder, sm); - - // Let's define the individual steps of the processing. We can experiment - // later with the ordering and the dependencies between them. - enum Steps { - BEGIN, - FILTER, - PINHEADS, - CLASSIFY, - ROUTING_GROUND, - ROUTING_NONGROUND, - CASCADE_PILLARS, - MERGE_RESULT, - DONE, - ABORT, - NUM_STEPS - //... - }; - - // Collect the algorithm steps into a nice sequence - std::array, NUM_STEPS> program = { - [] () { - // Begin... - // Potentially clear up the shared data (not needed for now) - }, - - std::bind(&SupportTreeBuildsteps::filter, &alg), - - std::bind(&SupportTreeBuildsteps::add_pinheads, &alg), - - std::bind(&SupportTreeBuildsteps::classify, &alg), - - std::bind(&SupportTreeBuildsteps::routing_to_ground, &alg), - - std::bind(&SupportTreeBuildsteps::routing_to_model, &alg), - - std::bind(&SupportTreeBuildsteps::interconnect_pillars, &alg), - - std::bind(&SupportTreeBuildsteps::merge_result, &alg), - - [] () { - // Done - }, - - [] () { - // Abort - } - }; - - Steps pc = BEGIN; - - if(sm.cfg.ground_facing_only) { - program[ROUTING_NONGROUND] = []() { - BOOST_LOG_TRIVIAL(info) - << "Skipping model-facing supports as requested."; - }; - } - - // Let's define a simple automaton that will run our program. - auto progress = [&builder, &pc] () { - static const std::array stepstr { - "Starting", - "Filtering", - "Generate pinheads", - "Classification", - "Routing to ground", - "Routing supports to model surface", - "Interconnecting pillars", - "Merging support mesh", - "Done", - "Abort" - }; - - static const std::array stepstate { - 0, - 10, - 30, - 50, - 60, - 70, - 80, - 99, - 100, - 0 - }; - - if(builder.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 = CASCADE_PILLARS; break; - case CASCADE_PILLARS: pc = MERGE_RESULT; break; - case MERGE_RESULT: pc = DONE; break; - case DONE: - case ABORT: break; - default: ; - } - - builder.ctl().statuscb(stepstate[pc], stepstr[pc]); - }; - - // Just here we run the computation... - while(pc < DONE) { - progress(); - program[pc](); - } - - return pc == ABORT; -} - -IndexedMesh::hit_result SupportTreeBuildsteps::pinhead_mesh_intersect( - const Vec3d &s, - const Vec3d &dir, - double r_pin, - double r_back, - double width, - double sd) -{ - static const size_t SAMPLES = 8; - - // Move away slightly from the touching point to avoid raycasting on the - // inner surface of the mesh. - - auto& m = m_mesh; - using HitResult = IndexedMesh::hit_result; - - // Hit results - std::array hits; - - struct Rings { - double rpin; - double rback; - Vec3d spin; - Vec3d sback; - PointRing ring; - - Vec3d backring(size_t idx) { return ring.get(idx, sback, rback); } - Vec3d pinring(size_t idx) { return ring.get(idx, spin, rpin); } - } rings {r_pin + sd, r_back + sd, s, s + width * dir, dir}; - - // We will shoot multiple rays from the head pinpoint in the direction - // of the pinhead robe (side) surface. The result will be the smallest - // hit distance. - - ccr::for_each(size_t(0), hits.size(), - [&m, &rings, sd, &hits](size_t i) { - - // Point on the circle on the pin sphere - Vec3d ps = rings.pinring(i); - // This is the point on the circle on the back sphere - Vec3d p = rings.backring(i); - - auto &hit = hits[i]; - - // Point ps is not on mesh but can be inside or - // outside as well. This would cause many problems - // with ray-casting. To detect the position we will - // use the ray-casting result (which has an is_inside - // predicate). - - Vec3d n = (p - ps).normalized(); - auto q = m.query_ray_hit(ps + sd * n, n); - - if (q.is_inside()) { // the hit is inside the model - if (q.distance() > rings.rpin) { - // If we are inside the model and the hit - // distance is bigger than our pin circle - // diameter, it probably indicates that the - // support point was already inside the - // model, or there is really no space - // around the point. We will assign a zero - // hit distance to these cases which will - // enforce the function return value to be - // an invalid ray with zero hit distance. - // (see min_element at the end) - hit = HitResult(0.0); - } else { - // re-cast the ray from the outside of the - // object. The starting point has an offset - // of 2*safety_distance because the - // original ray has also had an offset - auto q2 = m.query_ray_hit(ps + (q.distance() + 2 * sd) * n, n); - hit = q2; - } - } else - hit = q; - }); - - return min_hit(hits); -} - -IndexedMesh::hit_result SupportTreeBuildsteps::bridge_mesh_intersect( - const Vec3d &src, const Vec3d &dir, double r, double sd) -{ - static const size_t SAMPLES = 8; - PointRing ring{dir}; - - using Hit = IndexedMesh::hit_result; - - // Hit results - std::array hits; - - ccr::for_each(size_t(0), hits.size(), - [this, r, src, /*ins_check,*/ &ring, dir, sd, &hits] (size_t i) - { - Hit &hit = hits[i]; - - // Point on the circle on the pin sphere - Vec3d p = ring.get(i, src, r + sd); - - auto hr = m_mesh.query_ray_hit(p + r * dir, dir); - - if(/*ins_check && */hr.is_inside()) { - if(hr.distance() > 2 * r + sd) hit = Hit(0.0); - else { - // re-cast the ray from the outside of the object - hit = m_mesh.query_ray_hit(p + (hr.distance() + EPSILON) * dir, dir); - } - } else hit = hr; - }); - - return min_hit(hits); -} - -bool SupportTreeBuildsteps::interconnect(const Pillar &pillar, - const Pillar &nextpillar) -{ - // We need to get the starting point of the zig-zag pattern. We have to - // be aware that the two head junctions are at different heights. We - // may start from the lowest junction and call it a day but this - // strategy would leave unconnected a lot of pillar duos where the - // shorter pillar is too short to start a new bridge but the taller - // pillar could still be bridged with the shorter one. - bool was_connected = false; - - Vec3d supper = pillar.startpoint(); - Vec3d slower = nextpillar.startpoint(); - Vec3d eupper = pillar.endpoint(); - Vec3d elower = nextpillar.endpoint(); - - double zmin = m_builder.ground_level + m_cfg.base_height_mm; - eupper(Z) = std::max(eupper(Z), zmin); - elower(Z) = std::max(elower(Z), zmin); - - // The usable length of both pillars should be positive - if(slower(Z) - elower(Z) < 0) return false; - if(supper(Z) - eupper(Z) < 0) return false; - - double pillar_dist = distance(Vec2d{slower(X), slower(Y)}, - Vec2d{supper(X), supper(Y)}); - double bridge_distance = pillar_dist / std::cos(-m_cfg.bridge_slope); - double zstep = pillar_dist * std::tan(-m_cfg.bridge_slope); - - if(pillar_dist < 2 * m_cfg.head_back_radius_mm || - pillar_dist > m_cfg.max_pillar_link_distance_mm) return false; - - if(supper(Z) < slower(Z)) supper.swap(slower); - if(eupper(Z) < elower(Z)) eupper.swap(elower); - - double startz = 0, endz = 0; - - startz = slower(Z) - zstep < supper(Z) ? slower(Z) - zstep : slower(Z); - endz = eupper(Z) + zstep > elower(Z) ? eupper(Z) + zstep : eupper(Z); - - if(slower(Z) - eupper(Z) < std::abs(zstep)) { - // no space for even one cross - - // Get max available space - startz = std::min(supper(Z), slower(Z) - zstep); - endz = std::max(eupper(Z) + zstep, elower(Z)); - - // Align to center - double available_dist = (startz - endz); - double rounds = std::floor(available_dist / std::abs(zstep)); - startz -= 0.5 * (available_dist - rounds * std::abs(zstep)); - } - - auto pcm = m_cfg.pillar_connection_mode; - bool docrosses = - pcm == PillarConnectionMode::cross || - (pcm == PillarConnectionMode::dynamic && - pillar_dist > 2*m_cfg.base_radius_mm); - - // 'sj' means starting junction, 'ej' is the end junction of a bridge. - // They will be swapped in every iteration thus the zig-zag pattern. - // According to a config parameter, a second bridge may be added which - // results in a cross connection between the pillars. - Vec3d sj = supper, ej = slower; sj(Z) = startz; ej(Z) = sj(Z) + zstep; - - // TODO: This is a workaround to not have a faulty last bridge - while(ej(Z) >= eupper(Z) /*endz*/) { - if(bridge_mesh_distance(sj, dirv(sj, ej), pillar.r) >= bridge_distance) - { - m_builder.add_crossbridge(sj, ej, pillar.r); - was_connected = true; - } - - // double bridging: (crosses) - if(docrosses) { - Vec3d sjback(ej(X), ej(Y), sj(Z)); - Vec3d ejback(sj(X), sj(Y), ej(Z)); - if (sjback(Z) <= slower(Z) && ejback(Z) >= eupper(Z) && - bridge_mesh_distance(sjback, dirv(sjback, ejback), - pillar.r) >= bridge_distance) { - // need to check collision for the cross stick - m_builder.add_crossbridge(sjback, ejback, pillar.r); - was_connected = true; - } - } - - sj.swap(ej); - ej(Z) = sj(Z) + zstep; - } - - return was_connected; -} - -bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head, - long nearpillar_id) -{ - auto nearpillar = [this, nearpillar_id]() -> const Pillar& { - return m_builder.pillar(nearpillar_id); - }; - - if (m_builder.bridgecount(nearpillar()) > m_cfg.max_bridges_on_pillar) - return false; - - Vec3d headjp = head.junction_point(); - Vec3d nearjp_u = nearpillar().startpoint(); - Vec3d nearjp_l = nearpillar().endpoint(); - - double r = head.r_back_mm; - double d2d = distance(to_2d(headjp), to_2d(nearjp_u)); - double d3d = distance(headjp, nearjp_u); - - double hdiff = nearjp_u(Z) - headjp(Z); - double slope = std::atan2(hdiff, d2d); - - Vec3d bridgestart = headjp; - Vec3d bridgeend = nearjp_u; - double max_len = r * m_cfg.max_bridge_length_mm / m_cfg.head_back_radius_mm; - double max_slope = m_cfg.bridge_slope; - double zdiff = 0.0; - - // check the default situation if feasible for a bridge - if(d3d > max_len || slope > -max_slope) { - // not feasible to connect the two head junctions. We have to search - // for a suitable touch point. - - double Zdown = headjp(Z) + d2d * std::tan(-max_slope); - Vec3d touchjp = bridgeend; touchjp(Z) = Zdown; - double D = distance(headjp, touchjp); - zdiff = Zdown - nearjp_u(Z); - - if(zdiff > 0) { - Zdown -= zdiff; - bridgestart(Z) -= zdiff; - touchjp(Z) = Zdown; - - double t = bridge_mesh_distance(headjp, DOWN, r); - - // We can't insert a pillar under the source head to connect - // with the nearby pillar's starting junction - if(t < zdiff) return false; - } - - if(Zdown <= nearjp_u(Z) && Zdown >= nearjp_l(Z) && D < max_len) - bridgeend(Z) = Zdown; - else - return false; - } - - // There will be a minimum distance from the ground where the - // bridge is allowed to connect. This is an empiric value. - double minz = m_builder.ground_level + 4 * head.r_back_mm; - if(bridgeend(Z) < minz) return false; - - double t = bridge_mesh_distance(bridgestart, dirv(bridgestart, bridgeend), r); - - // Cannot insert the bridge. (further search might not worth the hassle) - if(t < distance(bridgestart, bridgeend)) return false; - - std::lock_guard lk(m_bridge_mutex); - - if (m_builder.bridgecount(nearpillar()) < m_cfg.max_bridges_on_pillar) { - // A partial pillar is needed under the starting head. - if(zdiff > 0) { - m_builder.add_pillar(head.id, headjp.z() - bridgestart.z()); - m_builder.add_junction(bridgestart, r); - m_builder.add_bridge(bridgestart, bridgeend, r); - } else { - m_builder.add_bridge(head.id, bridgeend); - } - - m_builder.increment_bridges(nearpillar()); - } else return false; - - return true; -} - -bool SupportTreeBuildsteps::create_ground_pillar(const Vec3d &hjp, - const Vec3d &sourcedir, - double radius, - long head_id) -{ - Vec3d jp = hjp, endp = jp, dir = sourcedir; - long pillar_id = SupportTreeNode::ID_UNSET; - bool can_add_base = false, non_head = false; - - double gndlvl = 0.; // The Z level where pedestals should be - double jp_gnd = 0.; // The lowest Z where a junction center can be - double gap_dist = 0.; // The gap distance between the model and the pad - - auto to_floor = [&gndlvl](const Vec3d &p) { return Vec3d{p.x(), p.y(), gndlvl}; }; - - auto eval_limits = [this, &radius, &can_add_base, &gndlvl, &gap_dist, &jp_gnd] - (bool base_en = true) - { - can_add_base = base_en && radius >= m_cfg.head_back_radius_mm; - double base_r = can_add_base ? m_cfg.base_radius_mm : 0.; - gndlvl = m_builder.ground_level; - if (!can_add_base) gndlvl -= m_mesh.ground_level_offset(); - jp_gnd = gndlvl + (can_add_base ? 0. : m_cfg.head_back_radius_mm); - gap_dist = m_cfg.pillar_base_safety_distance_mm + base_r + EPSILON; - }; - - eval_limits(); - - // We are dealing with a mini pillar that's potentially too long - if (radius < m_cfg.head_back_radius_mm && jp.z() - gndlvl > 20 * radius) - { - std::optional diffbr = - search_widening_path(jp, dir, radius, m_cfg.head_back_radius_mm); - - if (diffbr && diffbr->endp.z() > jp_gnd) { - auto &br = m_builder.add_diffbridge(*diffbr); - if (head_id >= 0) m_builder.head(head_id).bridge_id = br.id; - endp = diffbr->endp; - radius = diffbr->end_r; - m_builder.add_junction(endp, radius); - non_head = true; - dir = diffbr->get_dir(); - eval_limits(); - } else return false; - } - - if (m_cfg.object_elevation_mm < EPSILON) - { - // get a suitable direction for the corrector bridge. It is the - // original sourcedir's azimuth but the polar angle is saturated to the - // configured bridge slope. - auto [polar, azimuth] = dir_to_spheric(dir); - polar = PI - m_cfg.bridge_slope; - Vec3d d = spheric_to_dir(polar, azimuth).normalized(); - double t = bridge_mesh_distance(endp, d, radius); - double tmax = std::min(m_cfg.max_bridge_length_mm, t); - t = 0.; - - double zd = endp.z() - jp_gnd; - double tmax2 = zd / std::sqrt(1 - m_cfg.bridge_slope * m_cfg.bridge_slope); - tmax = std::min(tmax, tmax2); - - Vec3d nexp = endp; - double dlast = 0.; - while (((dlast = std::sqrt(m_mesh.squared_distance(to_floor(nexp)))) < gap_dist || - !std::isinf(bridge_mesh_distance(nexp, DOWN, radius))) && t < tmax) { - t += radius; - nexp = endp + t * d; - } - - if (dlast < gap_dist && can_add_base) { - nexp = endp; - t = 0.; - can_add_base = false; - eval_limits(can_add_base); - - zd = endp.z() - jp_gnd; - tmax2 = zd / std::sqrt(1 - m_cfg.bridge_slope * m_cfg.bridge_slope); - tmax = std::min(tmax, tmax2); - - while (((dlast = std::sqrt(m_mesh.squared_distance(to_floor(nexp)))) < gap_dist || - !std::isinf(bridge_mesh_distance(nexp, DOWN, radius))) && t < tmax) { - t += radius; - nexp = endp + t * d; - } - } - - // Could not find a path to avoid the pad gap - if (dlast < gap_dist) return false; - - if (t > 0.) { // Need to make additional bridge - const Bridge& br = m_builder.add_bridge(endp, nexp, radius); - if (head_id >= 0) m_builder.head(head_id).bridge_id = br.id; - - m_builder.add_junction(nexp, radius); - endp = nexp; - non_head = true; - } - } - - Vec3d gp = to_floor(endp); - double h = endp.z() - gp.z(); - - pillar_id = head_id >= 0 && !non_head ? m_builder.add_pillar(head_id, h) : - m_builder.add_pillar(gp, h, radius); - - if (can_add_base) - add_pillar_base(pillar_id); - - if(pillar_id >= 0) // Save the pillar endpoint in the spatial index - m_pillar_index.guarded_insert(m_builder.pillar(pillar_id).endpt, - unsigned(pillar_id)); - - return true; -} - -std::optional SupportTreeBuildsteps::search_widening_path( - const Vec3d &jp, const Vec3d &dir, double radius, double new_radius) -{ - double w = radius + 2 * m_cfg.head_back_radius_mm; - double stopval = w + jp.z() - m_builder.ground_level; - Optimizer solver(get_criteria(m_cfg).stop_score(stopval)); - - auto [polar, azimuth] = dir_to_spheric(dir); - - double fallback_ratio = radius / m_cfg.head_back_radius_mm; - - auto oresult = solver.to_max().optimize( - [this, jp, radius, new_radius](const opt::Input<3> &input) { - auto &[plr, azm, t] = input; - - auto d = spheric_to_dir(plr, azm).normalized(); - double ret = pinhead_mesh_intersect(jp, d, radius, new_radius, t) - .distance(); - double down = bridge_mesh_distance(jp + t * d, d, new_radius); - - if (ret > t && std::isinf(down)) - ret += jp.z() - m_builder.ground_level; - - return ret; - }, - initvals({polar, azimuth, w}), // start with what we have - bounds({ - {PI - m_cfg.bridge_slope, PI}, // Must not exceed the slope limit - {-PI, PI}, // azimuth can be a full search - {radius + m_cfg.head_back_radius_mm, - fallback_ratio * m_cfg.max_bridge_length_mm} - })); - - if (oresult.score >= stopval) { - polar = std::get<0>(oresult.optimum); - azimuth = std::get<1>(oresult.optimum); - double t = std::get<2>(oresult.optimum); - Vec3d endp = jp + t * spheric_to_dir(polar, azimuth); - - return DiffBridge(jp, endp, radius, m_cfg.head_back_radius_mm); - } - - return {}; -} - -void SupportTreeBuildsteps::filter() -{ - // Get the points that are too close to each other and keep only the - // first one - auto aliases = cluster(m_points, D_SP, 2); - - PtIndices filtered_indices; - filtered_indices.reserve(aliases.size()); - m_iheads.reserve(aliases.size()); - m_iheadless.reserve(aliases.size()); - for(auto& a : aliases) { - // Here we keep only the front point of the cluster. - filtered_indices.emplace_back(a.front()); - } - - // calculate the normals to the triangles for filtered points - auto nmls = sla::normals(m_points, m_mesh, m_cfg.head_front_radius_mm, - m_thr, filtered_indices); - - // 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. - - std::vector heads; heads.reserve(m_support_pts.size()); - for (const SupportPoint &sp : m_support_pts) { - m_thr(); - heads.emplace_back( - NaNd, - sp.head_front_radius, - 0., - m_cfg.head_penetration_mm, - Vec3d::Zero(), // dir - sp.pos.cast() // displacement - ); - } - - std::function filterfn; - filterfn = [this, &nmls, &heads, &filterfn](unsigned fidx, size_t i, double back_r) { - m_thr(); - - auto n = nmls.row(Eigen::Index(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. - - auto [polar, azimuth] = dir_to_spheric(n); - - // skip if the tilt is not sane - if (polar < PI - m_cfg.normal_cutoff_angle) return; - - // We saturate the polar angle to 3pi/4 - polar = std::max(polar, PI - m_cfg.bridge_slope); - - // save the head (pinpoint) position - Vec3d hp = m_points.row(fidx); - - double lmin = m_cfg.head_width_mm, lmax = lmin; - - if (back_r < m_cfg.head_back_radius_mm) { - lmin = 0., lmax = m_cfg.head_penetration_mm; - } - - // The distance needed for a pinhead to not collide with model. - double w = lmin + 2 * back_r + 2 * m_cfg.head_front_radius_mm - - m_cfg.head_penetration_mm; - - double pin_r = double(m_support_pts[fidx].head_front_radius); - - // Reassemble the now corrected normal - auto nn = spheric_to_dir(polar, azimuth).normalized(); - - // check available distance - IndexedMesh::hit_result t = pinhead_mesh_intersect(hp, nn, pin_r, - back_r, w); - - if (t.distance() < w) { - // Let's try to optimize this angle, there might be a - // viable normal that doesn't collide with the model - // geometry and its very close to the default. - - Optimizer solver(get_criteria(m_cfg)); - solver.seed(0); // we want deterministic behavior - - auto oresult = solver.to_max().optimize( - [this, pin_r, back_r, hp](const opt::Input<3> &input) - { - auto &[plr, azm, l] = input; - - auto dir = spheric_to_dir(plr, azm).normalized(); - - return pinhead_mesh_intersect( - hp, dir, pin_r, back_r, l).distance(); - }, - initvals({polar, azimuth, (lmin + lmax) / 2.}), // start with what we have - bounds({ - {PI - m_cfg.bridge_slope, PI}, // Must not exceed the slope limit - {-PI, PI}, // azimuth can be a full search - {lmin, lmax} - })); - - if(oresult.score > w) { - polar = std::get<0>(oresult.optimum); - azimuth = std::get<1>(oresult.optimum); - nn = spheric_to_dir(polar, azimuth).normalized(); - lmin = std::get<2>(oresult.optimum); - t = IndexedMesh::hit_result(oresult.score); - } - } - - if (t.distance() > w && hp(Z) + w * nn(Z) >= m_builder.ground_level) { - Head &h = heads[fidx]; - h.id = fidx; h.dir = nn; h.width_mm = lmin; h.r_back_mm = back_r; - } else if (back_r > m_cfg.head_fallback_radius_mm) { - filterfn(fidx, i, m_cfg.head_fallback_radius_mm); - } - }; - - ccr::for_each(size_t(0), filtered_indices.size(), - [this, &filterfn, &filtered_indices] (size_t i) { - filterfn(filtered_indices[i], i, m_cfg.head_back_radius_mm); - }); - - for (size_t i = 0; i < heads.size(); ++i) - if (heads[i].is_valid()) { - m_builder.add_head(i, heads[i]); - m_iheads.emplace_back(i); - } - - m_thr(); -} - -void SupportTreeBuildsteps::add_pinheads() -{ -} - -void SupportTreeBuildsteps::classify() -{ - // We should first get the heads that reach the ground directly - PtIndices ground_head_indices; - ground_head_indices.reserve(m_iheads.size()); - m_iheads_onmodel.reserve(m_iheads.size()); - - // First we decide which heads reach the ground and can be full - // pillars and which shall be connected to the model surface (or - // search a suitable path around the surface that leads to the - // ground -- TODO) - for(unsigned i : m_iheads) { - m_thr(); - - Head &head = m_builder.head(i); - double r = head.r_back_mm; - Vec3d headjp = head.junction_point(); - - // collision check - auto hit = bridge_mesh_intersect(headjp, DOWN, r); - - if(std::isinf(hit.distance())) ground_head_indices.emplace_back(i); - else if(m_cfg.ground_facing_only) head.invalidate(); - else m_iheads_onmodel.emplace_back(i); - - m_head_to_ground_scans[i] = hit; - } - - // 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 pointfn = [this](unsigned i) { - return m_builder.head(i).junction_point(); - }; - - auto predicate = [this](const PointIndexEl &e1, - const PointIndexEl &e2) { - double d2d = distance(to_2d(e1.first), to_2d(e2.first)); - double d3d = distance(e1.first, e2.first); - return d2d < 2 * m_cfg.base_radius_mm - && d3d < m_cfg.max_bridge_length_mm; - }; - - m_pillar_clusters = cluster(ground_head_indices, pointfn, predicate, - m_cfg.max_bridges_on_pillar); -} - -void SupportTreeBuildsteps::routing_to_ground() -{ - ClusterEl cl_centroids; - cl_centroids.reserve(m_pillar_clusters.size()); - - for (auto &cl : m_pillar_clusters) { - m_thr(); - - // 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. - - if (cl.empty()) continue; - - // get the current cluster centroid - auto & thr = m_thr; - const auto &points = m_points; - - long lcid = cluster_centroid( - cl, [&points](size_t idx) { return points.row(long(idx)); }, - [thr](const Vec3d &p1, const Vec3d &p2) { - thr(); - return distance(Vec2d(p1(X), p1(Y)), Vec2d(p2(X), p2(Y))); - }); - - assert(lcid >= 0); - unsigned hid = cl[size_t(lcid)]; // Head ID - - cl_centroids.emplace_back(hid); - - Head &h = m_builder.head(hid); - - if (!create_ground_pillar(h.junction_point(), h.dir, h.r_back_mm, h.id)) { - BOOST_LOG_TRIVIAL(warning) - << "Pillar cannot be created for support point id: " << hid; - m_iheads_onmodel.emplace_back(h.id); - continue; - } - } - - // 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. - size_t ci = 0; - for (auto cl : m_pillar_clusters) { - m_thr(); - - auto cidx = cl_centroids[ci++]; - - auto q = m_pillar_index.query(m_builder.head(cidx).junction_point(), 1); - if (!q.empty()) { - long centerpillarID = q.front().second; - for (auto c : cl) { - m_thr(); - if (c == cidx) continue; - - auto &sidehead = m_builder.head(c); - - if (!connect_to_nearpillar(sidehead, centerpillarID) && - !search_pillar_and_connect(sidehead)) { - Vec3d pstart = sidehead.junction_point(); - // Vec3d pend = Vec3d{pstart(X), pstart(Y), gndlvl}; - // Could not find a pillar, create one - create_ground_pillar(pstart, sidehead.dir, sidehead.r_back_mm, sidehead.id); - } - } - } - } -} - -bool SupportTreeBuildsteps::connect_to_ground(Head &head, const Vec3d &dir) -{ - auto hjp = head.junction_point(); - double r = head.r_back_mm; - double t = bridge_mesh_distance(hjp, dir, head.r_back_mm); - double d = 0, tdown = 0; - t = std::min(t, m_cfg.max_bridge_length_mm * r / m_cfg.head_back_radius_mm); - - while (d < t && !std::isinf(tdown = bridge_mesh_distance(hjp + d * dir, DOWN, r))) - d += r; - - if(!std::isinf(tdown)) return false; - - Vec3d endp = hjp + d * dir; - bool ret = false; - - if ((ret = create_ground_pillar(endp, dir, head.r_back_mm))) { - m_builder.add_bridge(head.id, endp); - m_builder.add_junction(endp, head.r_back_mm); - } - - return ret; -} - -bool SupportTreeBuildsteps::connect_to_ground(Head &head) -{ - if (connect_to_ground(head, head.dir)) return true; - - // Optimize bridge direction: - // Straight path failed so we will try to search for a suitable - // direction out of the cavity. - auto [polar, azimuth] = dir_to_spheric(head.dir); - - Optimizer solver(get_criteria(m_cfg).stop_score(1e6)); - solver.seed(0); // we want deterministic behavior - - double r_back = head.r_back_mm; - Vec3d hjp = head.junction_point(); - auto oresult = solver.to_max().optimize( - [this, hjp, r_back](const opt::Input<2> &input) { - auto &[plr, azm] = input; - Vec3d n = spheric_to_dir(plr, azm).normalized(); - return bridge_mesh_distance(hjp, n, r_back); - }, - initvals({polar, azimuth}), // let's start with what we have - bounds({ {PI - m_cfg.bridge_slope, PI}, {-PI, PI} }) - ); - - Vec3d bridgedir = spheric_to_dir(oresult.optimum).normalized(); - return connect_to_ground(head, bridgedir); -} - -bool SupportTreeBuildsteps::connect_to_model_body(Head &head) -{ - if (head.id <= SupportTreeNode::ID_UNSET) return false; - - auto it = m_head_to_ground_scans.find(unsigned(head.id)); - if (it == m_head_to_ground_scans.end()) return false; - - auto &hit = it->second; - - if (!hit.is_hit()) { - // TODO scan for potential anchor points on model surface - return false; - } - - Vec3d hjp = head.junction_point(); - double zangle = std::asin(hit.direction()(Z)); - zangle = std::max(zangle, PI/4); - double h = std::sin(zangle) * head.fullwidth(); - - // The width of the tail head that we would like to have... - h = std::min(hit.distance() - head.r_back_mm, h); - - // If this is a mini pillar dont bother with the tail width, can be 0. - if (head.r_back_mm < m_cfg.head_back_radius_mm) h = std::max(h, 0.); - else if (h <= 0.) return false; - - Vec3d endp{hjp(X), hjp(Y), hjp(Z) - hit.distance() + h}; - auto center_hit = m_mesh.query_ray_hit(hjp, DOWN); - - double hitdiff = center_hit.distance() - hit.distance(); - Vec3d hitp = std::abs(hitdiff) < 2*head.r_back_mm? - center_hit.position() : hit.position(); - - long pillar_id = m_builder.add_pillar(head.id, hjp.z() - endp.z()); - Pillar &pill = m_builder.pillar(pillar_id); - - Vec3d taildir = endp - hitp; - double dist = (hitp - endp).norm() + m_cfg.head_penetration_mm; - double w = dist - 2 * head.r_pin_mm - head.r_back_mm; - - if (w < 0.) { - BOOST_LOG_TRIVIAL(warning) << "Pinhead width is negative!"; - w = 0.; - } - - m_builder.add_anchor(head.r_back_mm, head.r_pin_mm, w, - m_cfg.head_penetration_mm, taildir, hitp); - - m_pillar_index.guarded_insert(pill.endpoint(), pill.id); - - return true; -} - -bool SupportTreeBuildsteps::search_pillar_and_connect(const Head &source) -{ - // Hope that a local copy takes less time than the whole search loop. - // We also need to remove elements progressively from the copied index. - PointIndex spindex = m_pillar_index.guarded_clone(); - - long nearest_id = SupportTreeNode::ID_UNSET; - - Vec3d querypt = source.junction_point(); - - while(nearest_id < 0 && !spindex.empty()) { m_thr(); - // 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(querypt(X), querypt(Y), m_builder.ground_level); - auto qres = spindex.nearest(qp, 1); - if(qres.empty()) break; - - auto ne = qres.front(); - nearest_id = ne.second; - - if(nearest_id >= 0) { - if (size_t(nearest_id) < m_builder.pillarcount()) { - if(!connect_to_nearpillar(source, nearest_id) || - m_builder.pillar(nearest_id).r < source.r_back_mm) { - nearest_id = SupportTreeNode::ID_UNSET; // continue searching - spindex.remove(ne); // without the current pillar - } - } - } - } - - return nearest_id >= 0; -} - -void SupportTreeBuildsteps::routing_to_model() -{ - // We need to check if there is an easy way out to the bed surface. - // If it can be routed there with a bridge shorter than - // min_bridge_distance. - - ccr::for_each(m_iheads_onmodel.begin(), m_iheads_onmodel.end(), - [this] (const unsigned idx) { - m_thr(); - - auto& head = m_builder.head(idx); - - // Search nearby pillar - if (search_pillar_and_connect(head)) { return; } - - // Cannot connect to nearby pillar. We will try to search for - // a route to the ground. - if (connect_to_ground(head)) { return; } - - // No route to the ground, so connect to the model body as a last resort - if (connect_to_model_body(head)) { return; } - - // We have failed to route this head. - BOOST_LOG_TRIVIAL(warning) - << "Failed to route model facing support point. ID: " << idx; - - head.invalidate(); - }); -} - -void SupportTreeBuildsteps::interconnect_pillars() -{ - // Now comes the algorithm that connects pillars with each other. - // Ideally every pillar should be connected with at least one of its - // neighbors if that neighbor is within max_pillar_link_distance - - // Pillars with height exceeding H1 will require at least one neighbor - // to connect with. Height exceeding H2 require two neighbors. - double H1 = m_cfg.max_solo_pillar_height_mm; - double H2 = m_cfg.max_dual_pillar_height_mm; - double d = m_cfg.max_pillar_link_distance_mm; - - //A connection between two pillars only counts if the height ratio is - // bigger than 50% - double min_height_ratio = 0.5; - - std::set pairs; - - // A function to connect one pillar with its neighbors. THe number of - // neighbors is given in the configuration. This function if called - // for every pillar in the pillar index. A pair of pillar will not - // be connected multiple times this is ensured by the 'pairs' set which - // remembers the processed pillar pairs - auto cascadefn = - [this, d, &pairs, min_height_ratio, H1] (const PointIndexEl& el) - { - Vec3d qp = el.first; // endpoint of the pillar - - const Pillar& pillar = m_builder.pillar(el.second); // actual pillar - - // Get the max number of neighbors a pillar should connect to - unsigned neighbors = m_cfg.pillar_cascade_neighbors; - - // connections are already enough for the pillar - if(pillar.links >= neighbors) return; - - double max_d = d * pillar.r / m_cfg.head_back_radius_mm; - // Query all remaining points within reach - auto qres = m_pillar_index.query([qp, max_d](const PointIndexEl& e){ - return distance(e.first, qp) < max_d; - }); - - // sort the result by distance (have to check if this is needed) - std::sort(qres.begin(), qres.end(), - [qp](const PointIndexEl& e1, const PointIndexEl& e2){ - return distance(e1.first, qp) < distance(e2.first, qp); - }); - - for(auto& re : qres) { // process the queried neighbors - - if(re.second == el.second) continue; // Skip self - - auto a = el.second, b = re.second; - - // Get unique hash for the given pair (order doesn't matter) - auto hashval = pairhash(a, b); - - // Search for the pair amongst the remembered pairs - if(pairs.find(hashval) != pairs.end()) continue; - - const Pillar& neighborpillar = m_builder.pillar(re.second); - - // this neighbor is occupied, skip - if (neighborpillar.links >= neighbors) continue; - if (neighborpillar.r < pillar.r) continue; - - if(interconnect(pillar, neighborpillar)) { - pairs.insert(hashval); - - // If the interconnection length between the two pillars is - // less than 50% of the longer pillar's height, don't count - if(pillar.height < H1 || - neighborpillar.height / pillar.height > min_height_ratio) - m_builder.increment_links(pillar); - - if(neighborpillar.height < H1 || - pillar.height / neighborpillar.height > min_height_ratio) - m_builder.increment_links(neighborpillar); - - } - - // connections are enough for one pillar - if(pillar.links >= neighbors) break; - } - }; - - // Run the cascade for the pillars in the index - m_pillar_index.foreach(cascadefn); - - // We would be done here if we could allow some pillars to not be - // connected with any neighbors. But this might leave the support tree - // unprintable. - // - // The current solution is to insert additional pillars next to these - // lonely pillars. One or even two additional pillar might get inserted - // depending on the length of the lonely pillar. - - size_t pillarcount = m_builder.pillarcount(); - - // Again, go through all pillars, this time in the whole support tree - // not just the index. - for(size_t pid = 0; pid < pillarcount; pid++) { - auto pillar = [this, pid]() { return m_builder.pillar(pid); }; - - // Decide how many additional pillars will be needed: - - unsigned needpillars = 0; - if (pillar().bridges > m_cfg.max_bridges_on_pillar) - needpillars = 3; - else if (pillar().links < 2 && pillar().height > H2) { - // Not enough neighbors to support this pillar - needpillars = 2; - } else if (pillar().links < 1 && pillar().height > H1) { - // No neighbors could be found and the pillar is too long. - needpillars = 1; - } - - needpillars = std::max(pillar().links, needpillars) - pillar().links; - if (needpillars == 0) continue; - - // Search for new pillar locations: - - bool found = false; - double alpha = 0; // goes to 2Pi - double r = 2 * m_cfg.base_radius_mm; - Vec3d pillarsp = pillar().startpoint(); - - // temp value for starting point detection - Vec3d sp(pillarsp(X), pillarsp(Y), pillarsp(Z) - r); - - // A vector of bool for placement feasbility - std::vector canplace(needpillars, false); - std::vector spts(needpillars); // vector of starting points - - double gnd = m_builder.ground_level; - double min_dist = m_cfg.pillar_base_safety_distance_mm + - m_cfg.base_radius_mm + EPSILON; - - while(!found && alpha < 2*PI) { - for (unsigned n = 0; - n < needpillars && (!n || canplace[n - 1]); - n++) - { - double a = alpha + n * PI / 3; - Vec3d s = sp; - s(X) += std::cos(a) * r; - s(Y) += std::sin(a) * r; - spts[n] = s; - - // Check the path vertically down - Vec3d check_from = s + Vec3d{0., 0., pillar().r}; - auto hr = bridge_mesh_intersect(check_from, DOWN, pillar().r); - Vec3d gndsp{s(X), s(Y), gnd}; - - // If the path is clear, check for pillar base collisions - canplace[n] = std::isinf(hr.distance()) && - std::sqrt(m_mesh.squared_distance(gndsp)) > - min_dist; - } - - found = std::all_of(canplace.begin(), canplace.end(), - [](bool v) { return v; }); - - // 20 angles will be tried... - alpha += 0.1 * PI; - } - - std::vector newpills; - newpills.reserve(needpillars); - - if (found) - for (unsigned n = 0; n < needpillars; n++) { - Vec3d s = spts[n]; - Pillar p(Vec3d{s.x(), s.y(), gnd}, s.z() - gnd, pillar().r); - - if (interconnect(pillar(), p)) { - Pillar &pp = m_builder.pillar(m_builder.add_pillar(p)); - - add_pillar_base(pp.id); - - m_pillar_index.insert(pp.endpoint(), unsigned(pp.id)); - - m_builder.add_junction(s, pillar().r); - double t = bridge_mesh_distance(pillarsp, dirv(pillarsp, s), - pillar().r); - if (distance(pillarsp, s) < t) - m_builder.add_bridge(pillarsp, s, pillar().r); - - if (pillar().endpoint()(Z) > m_builder.ground_level + pillar().r) - m_builder.add_junction(pillar().endpoint(), pillar().r); - - newpills.emplace_back(pp.id); - m_builder.increment_links(pillar()); - m_builder.increment_links(pp); - } - } - - if(!newpills.empty()) { - for(auto it = newpills.begin(), nx = std::next(it); - nx != newpills.end(); ++it, ++nx) { - const Pillar& itpll = m_builder.pillar(*it); - const Pillar& nxpll = m_builder.pillar(*nx); - if(interconnect(itpll, nxpll)) { - m_builder.increment_links(itpll); - m_builder.increment_links(nxpll); - } - } - - m_pillar_index.foreach(cascadefn); - } - } -} - -}} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTreeMesher.cpp b/src/libslic3r/SLA/SupportTreeMesher.cpp index 17f57f457..3f0b6c841 100644 --- a/src/libslic3r/SLA/SupportTreeMesher.cpp +++ b/src/libslic3r/SLA/SupportTreeMesher.cpp @@ -93,68 +93,6 @@ indexed_triangle_set sphere(double rho, Portion portion, double fa) { return ret; } -indexed_triangle_set cylinder(double r, double h, size_t ssteps, const Vec3d &sp) -{ - assert(ssteps > 0); - - indexed_triangle_set ret; - - auto steps = int(ssteps); - auto& points = ret.vertices; - auto& indices = ret.indices; - points.reserve(2*ssteps); - double a = 2*PI/steps; - - Vec3d jp = sp; - Vec3d endp = {sp(X), sp(Y), sp(Z) + h}; - - // Upper circle points - for(int i = 0; i < steps; ++i) { - double phi = i*a; - auto ex = float(endp(X) + r*std::cos(phi)); - auto ey = float(endp(Y) + r*std::sin(phi)); - points.emplace_back(ex, ey, float(endp(Z))); - } - - // Lower circle points - for(int i = 0; i < steps; ++i) { - double phi = i*a; - auto x = float(jp(X) + r*std::cos(phi)); - auto y = float(jp(Y) + r*std::sin(phi)); - points.emplace_back(x, y, float(jp(Z))); - } - - // Now create long triangles connecting upper and lower circles - indices.reserve(2*ssteps); - 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); - } - - // Last triangle connecting the first and last vertices - auto last = steps - 1; - indices.emplace_back(0, last, offs); - indices.emplace_back(last, offs + last, offs); - - // According to the slicing algorithms, we need to aid them with generating - // a watertight body. So we create a triangle fan for the upper and lower - // ending of the cylinder to close the geometry. - points.emplace_back(jp.cast()); int ci = int(points.size() - 1); - for(int i = 0; i < steps - 1; ++i) - indices.emplace_back(i + offs + 1, i + offs, ci); - - indices.emplace_back(offs, steps + offs - 1, ci); - - points.emplace_back(endp.cast()); ci = int(points.size() - 1); - for(int i = 0; i < steps - 1; ++i) - indices.emplace_back(ci, i, i + 1); - - indices.emplace_back(steps - 1, 0, ci); - - return ret; -} - indexed_triangle_set pinhead(double r_pin, double r_back, double length, @@ -183,13 +121,16 @@ indexed_triangle_set pinhead(double r_pin, const double h = r_back + r_pin + length; double phi = PI / 2. - std::acos((r_back - r_pin) / h); + if (std::isnan(phi)) + return mesh; + // 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_back, make_portion(0, PI / 2 + phi), detail); - auto &&s2 = sphere(r_pin, make_portion(PI / 2 + phi, PI), detail); + auto s1 = sphere(r_back, make_portion(0, PI / 2 + phi), detail); + auto s2 = sphere(r_pin, make_portion(PI / 2 + phi, PI), detail); for (auto &p : s2.vertices) p.z() += h; @@ -267,4 +208,63 @@ indexed_triangle_set halfcone(double baseheight, return base; } +indexed_triangle_set get_mesh(const Head &h, size_t steps) +{ + indexed_triangle_set mesh = pinhead(h.r_pin_mm, h.r_back_mm, h.width_mm, steps); + + for (auto& p : mesh.vertices) p.z() -= (h.fullwidth() - h.r_back_mm); + + 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(Vec3f{0.f, 0.f, -1.f}, + h.dir.cast()); + + Vec3f pos = h.pos.cast(); + for (auto& p : mesh.vertices) p = quatern * p + pos; + + return mesh; +} + +indexed_triangle_set get_mesh(const Bridge &br, size_t steps) +{ + using Quaternion = Eigen::Quaternion; + Vec3d v = (br.endp - br.startp); + Vec3d dir = v.normalized(); + double d = v.norm(); + + indexed_triangle_set mesh = cylinder(br.r, d, steps); + + auto quater = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f}, + dir.cast()); + + Vec3f startp = br.startp.cast(); + for(auto& p : mesh.vertices) p = quater * p + startp; + + return mesh; +} + +indexed_triangle_set get_mesh(const DiffBridge &br, size_t steps) +{ + double h = br.get_length(); + indexed_triangle_set mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps); + + 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(Vec3f{0.f, 0.f, 1.f}, + br.get_dir().cast()); + + Vec3f startp = br.startp.cast(); + for(auto& p : mesh.vertices) p = quatern * p + startp; + + return mesh; +} + }} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTreeMesher.hpp b/src/libslic3r/SLA/SupportTreeMesher.hpp index 63ccc5fce..c25972923 100644 --- a/src/libslic3r/SLA/SupportTreeMesher.hpp +++ b/src/libslic3r/SLA/SupportTreeMesher.hpp @@ -17,18 +17,20 @@ inline Portion make_portion(double a, double b) } indexed_triangle_set sphere(double rho, - Portion portion = make_portion(0., 2. * PI), + Portion portion = make_portion(0., PI), double fa = (2. * PI / 360.)); // Down facing cylinder in Z direction with arguments: // r: radius -// h: Height +// h: height // ssteps: how many edges will create the base circle // sp: starting point -indexed_triangle_set cylinder(double r, +inline indexed_triangle_set cylinder(double r, double h, - size_t steps = 45, - const Vec3d &sp = Vec3d::Zero()); + size_t steps = 45) +{ + return its_make_cylinder(r, h, 2 * PI / steps); +} indexed_triangle_set pinhead(double r_pin, double r_back, @@ -41,33 +43,15 @@ indexed_triangle_set halfcone(double baseheight, const Vec3d &pt = Vec3d::Zero(), size_t steps = 45); -inline indexed_triangle_set get_mesh(const Head &h, size_t steps) -{ - indexed_triangle_set mesh = pinhead(h.r_pin_mm, h.r_back_mm, h.width_mm, steps); - - for (auto& p : mesh.vertices) p.z() -= (h.fullwidth() - h.r_back_mm); - - 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(Vec3f{0.f, 0.f, -1.f}, - h.dir.cast()); - - Vec3f pos = h.pos.cast(); - for (auto& p : mesh.vertices) p = quatern * p + pos; - - return mesh; -} +indexed_triangle_set get_mesh(const Head &h, size_t steps); inline indexed_triangle_set get_mesh(const Pillar &p, size_t steps) { if(p.height > EPSILON) { // Endpoint is below the starting point // We just create a bridge geometry with the pillar parameters and // move the data. - return cylinder(p.r, p.height, steps, p.endpoint()); + //return cylinder(p.r_start, p.height, steps, p.endpoint()); + return halfcone(p.height, p.r_end, p.r_start, p.endpt, steps); } return {}; @@ -80,49 +64,15 @@ inline indexed_triangle_set get_mesh(const Pedestal &p, size_t steps) inline indexed_triangle_set get_mesh(const Junction &j, size_t steps) { - indexed_triangle_set mesh = sphere(j.r, make_portion(0, PI), 2 *PI / steps); + indexed_triangle_set mesh = sphere(j.r, make_portion(0, PI), 2 * PI / steps); Vec3f pos = j.pos.cast(); for(auto& p : mesh.vertices) p += pos; return mesh; } -inline indexed_triangle_set get_mesh(const Bridge &br, size_t steps) -{ - using Quaternion = Eigen::Quaternion; - Vec3d v = (br.endp - br.startp); - Vec3d dir = v.normalized(); - double d = v.norm(); +indexed_triangle_set get_mesh(const Bridge &br, size_t steps); - indexed_triangle_set mesh = cylinder(br.r, d, steps); - - auto quater = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f}, - dir.cast()); - - Vec3f startp = br.startp.cast(); - for(auto& p : mesh.vertices) p = quater * p + startp; - - return mesh; -} - -inline indexed_triangle_set get_mesh(const DiffBridge &br, size_t steps) -{ - double h = br.get_length(); - indexed_triangle_set mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps); - - 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(Vec3f{0.f, 0.f, 1.f}, - br.get_dir().cast()); - - Vec3f startp = br.startp.cast(); - for(auto& p : mesh.vertices) p = quatern * p + startp; - - return mesh; -} +indexed_triangle_set get_mesh(const DiffBridge &br, size_t steps); }} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/SupportTreeStrategies.hpp b/src/libslic3r/SLA/SupportTreeStrategies.hpp new file mode 100644 index 000000000..d04ceaa2d --- /dev/null +++ b/src/libslic3r/SLA/SupportTreeStrategies.hpp @@ -0,0 +1,13 @@ +#ifndef SUPPORTTREESTRATEGIES_HPP +#define SUPPORTTREESTRATEGIES_HPP + +#include + +namespace Slic3r { namespace sla { + +enum class SupportTreeType { Default, Clever }; +enum class PillarConnectionMode { zigzag, cross, dynamic }; + +}} // namespace Slic3r::sla + +#endif // SUPPORTTREESTRATEGIES_HPP diff --git a/src/libslic3r/SLA/SupportTreeUtils.hpp b/src/libslic3r/SLA/SupportTreeUtils.hpp new file mode 100644 index 000000000..9c3c1f584 --- /dev/null +++ b/src/libslic3r/SLA/SupportTreeUtils.hpp @@ -0,0 +1,799 @@ +#ifndef SLASUPPORTTREEUTILS_H +#define SLASUPPORTTREEUTILS_H + +#include +#include + +#include +#include +#include +#include +#include + +#include + +namespace Slic3r { namespace sla { + +using Slic3r::opt::initvals; +using Slic3r::opt::bounds; +using Slic3r::opt::StopCriteria; +using Slic3r::opt::Optimizer; +using Slic3r::opt::AlgNLoptSubplex; +using Slic3r::opt::AlgNLoptGenetic; +using Slic3r::Geometry::dir_to_spheric; +using Slic3r::Geometry::spheric_to_dir; + +// Helper function for pillar interconnection where pairs of already connected +// pillars should be checked for not to be processed again. This can be done +// in constant time with a set of hash values uniquely representing a pair of +// integers. The order of numbers within the pair should not matter, it has +// the same unique hash. The hash value has to have twice as many bits as the +// arguments need. If the same integral type is used for args and return val, +// make sure the arguments use only the half of the type's bit depth. +template> +IntegerOnly pairhash(I a, I b) +{ + using std::ceil; using std::log2; using std::max; using std::min; + static const auto constexpr Ibits = int(sizeof(I) * CHAR_BIT); + static const auto constexpr DoubleIbits = int(sizeof(DoubleI) * CHAR_BIT); + static const auto constexpr shift = DoubleIbits / 2 < Ibits ? Ibits / 2 : Ibits; + + I g = min(a, b), l = max(a, b); + + // Assume the hash will fit into the output variable + assert((g ? (ceil(log2(g))) : 0) <= shift); + assert((l ? (ceil(log2(l))) : 0) <= shift); + + return (DoubleI(g) << shift) + l; +} + +// Give points on a 3D ring with given center, radius and orientation +// method based on: +// https://math.stackexchange.com/questions/73237/parametric-equation-of-a-circle-in-3d-space +template +class PointRing { + std::array m_phis; + + // Two vectors that will be perpendicular to each other and to the + // axis. Values for a(X) and a(Y) are now arbitrary, a(Z) is just a + // placeholder. + // a and b vectors are perpendicular to the ring direction and to each other. + // Together they define the plane where we have to iterate with the + // given angles in the 'm_phis' vector + Vec3d a = {0, 1, 0}, b; + double m_radius = 0.; + + static inline bool constexpr is_one(double val) + { + constexpr double eps = 1e-20; + + return std::abs(std::abs(val) - 1) < eps; + } + +public: + + PointRing(const Vec3d &n) : m_phis{linspace_array(0., 2 * PI)} + { + // We have to address the case when the direction vector v (same as + // dir) is coincident with one of the world axes. In this case two of + // its components will be completely zero and one is 1.0. Our method + // becomes dangerous here due to division with zero. Instead, vector + // 'a' can be an element-wise rotated version of 'v' + if(is_one(n(X)) || is_one(n(Y)) || is_one(n(Z))) { + a = {n(Z), n(X), n(Y)}; + b = {n(Y), n(Z), n(X)}; + } + else { + a(Z) = -(n(Y)*a(Y)) / n(Z); a.normalize(); + b = a.cross(n); + } + } + + Vec3d get(size_t idx, const Vec3d &src, double r) const + { + double phi = m_phis[idx]; + double sinphi = std::sin(phi); + double cosphi = std::cos(phi); + + double rpscos = r * cosphi; + double rpssin = r * sinphi; + + // Point on the sphere + return {src(X) + rpscos * a(X) + rpssin * b(X), + src(Y) + rpscos * a(Y) + rpssin * b(Y), + src(Z) + rpscos * a(Z) + rpssin * b(Z)}; + } +}; + +template +Vec dirv(const Vec& startp, const Vec& endp) { + return (endp - startp).normalized(); +} + +using Hit = AABBMesh::hit_result; + +template Hit min_hit(It from, It to) +{ + auto mit = std::min_element(from, to, [](const Hit &h1, const Hit &h2) { + return h1.distance() < h2.distance(); + }); + + return *mit; +} + +inline StopCriteria get_criteria(const SupportTreeConfig &cfg) +{ + return StopCriteria{} + .rel_score_diff(cfg.optimizer_rel_score_diff) + .max_iterations(cfg.optimizer_max_iterations); +} + +// A simple sphere with a center and a radius +struct Ball { Vec3d p; double R; }; + +struct Beam { // Defines a set of rays displaced along a cone's surface + static constexpr size_t SAMPLES = 8; + + Vec3d src; + Vec3d dir; + double r1; + double r2; // radius of the beam 1 unit further from src in dir direction + + Beam(const Vec3d &s, const Vec3d &d, double R1, double R2): + src{s}, dir{d}, r1{R1}, r2{R2} {}; + + Beam(const Ball &src_ball, const Ball &dst_ball): + src{src_ball.p}, dir{dirv(src_ball.p, dst_ball.p)}, r1{src_ball.R} + { + r2 = src_ball.R + + (dst_ball.R - src_ball.R) / distance(src_ball.p, dst_ball.p); + } + + Beam(const Vec3d &s, const Vec3d &d, double R) + : src{s}, dir{d}, r1{R}, r2{R} + {} +}; + +template +Hit beam_mesh_hit(Ex ex, const AABBMesh &mesh, const Beam &beam, double sd) +{ + Vec3d src = beam.src; + Vec3d dst = src + beam.dir; + double r_src = beam.r1; + double r_dst = beam.r2; + + Vec3d D = (dst - src); + Vec3d dir = D.normalized(); + PointRing ring{dir}; + + using Hit = AABBMesh::hit_result; + + // Hit results + std::array hits; + + execution::for_each( + ex, size_t(0), hits.size(), + [&mesh, r_src, r_dst, src, dst, &ring, dir, sd, &hits](size_t i) { + Hit &hit = hits[i]; + + // Point on the circle on the pin sphere + Vec3d p_src = ring.get(i, src, r_src + sd); + Vec3d p_dst = ring.get(i, dst, r_dst + sd); + Vec3d raydir = (p_dst - p_src).normalized(); + + auto hr = mesh.query_ray_hit(p_src + r_src * raydir, raydir); + + if (hr.is_inside()) { + if (hr.distance() > 2 * r_src + sd) + hit = Hit(0.0); + else { + // re-cast the ray from the outside of the object + auto q = p_src + (hr.distance() + EPSILON) * raydir; + hit = mesh.query_ray_hit(q, raydir); + } + } else + hit = hr; + }, std::min(execution::max_concurrency(ex), Beam::SAMPLES)); + + return min_hit(hits.begin(), hits.end()); +} + +template +Hit pinhead_mesh_hit(Ex ex, + const AABBMesh &mesh, + const Vec3d &s, + const Vec3d &dir, + double r_pin, + double r_back, + double width, + double sd) +{ + static const size_t SAMPLES = 8; + + // Move away slightly from the touching point to avoid raycasting on the + // inner surface of the mesh. + + auto &m = mesh; + using HitResult = AABBMesh::hit_result; + + // Hit results + std::array hits; + + struct Rings + { + double rpin; + double rback; + Vec3d spin; + Vec3d sback; + PointRing ring; + + Vec3d backring(size_t idx) { return ring.get(idx, sback, rback); } + Vec3d pinring(size_t idx) { return ring.get(idx, spin, rpin); } + } rings{r_pin + sd, r_back + sd, s, s + (r_pin + width + r_back) * dir, dir}; + + // We will shoot multiple rays from the head pinpoint in the direction + // of the pinhead robe (side) surface. The result will be the smallest + // hit distance. + + execution::for_each( + ex, size_t(0), hits.size(), [&m, &rings, sd, &hits](size_t i) { + // Point on the circle on the pin sphere + Vec3d ps = rings.pinring(i); + // This is the point on the circle on the back sphere + Vec3d p = rings.backring(i); + + auto &hit = hits[i]; + + // Point ps is not on mesh but can be inside or + // outside as well. This would cause many problems + // with ray-casting. To detect the position we will + // use the ray-casting result (which has an is_inside + // predicate). + + Vec3d n = (p - ps).normalized(); + auto q = m.query_ray_hit(ps + sd * n, n); + + if (q.is_inside()) { // the hit is inside the model + if (q.distance() > rings.rpin) { + // If we are inside the model and the hit + // distance is bigger than our pin circle + // diameter, it probably indicates that the + // support point was already inside the + // model, or there is really no space + // around the point. We will assign a zero + // hit distance to these cases which will + // enforce the function return value to be + // an invalid ray with zero hit distance. + // (see min_element at the end) + hit = HitResult(0.0); + } else { + // re-cast the ray from the outside of the + // object. The starting point has an offset + // of 2*safety_distance because the + // original ray has also had an offset + auto q2 = m.query_ray_hit(ps + (q.distance() + 2 * sd) * n, n); + hit = q2; + } + } else + hit = q; + }, std::min(execution::max_concurrency(ex), SAMPLES)); + + return min_hit(hits.begin(), hits.end()); +} + +template +Hit pinhead_mesh_hit(Ex ex, + const AABBMesh &mesh, + const Head &head, + double safety_d) +{ + return pinhead_mesh_hit(ex, mesh, head.pos, head.dir, head.r_pin_mm, + head.r_back_mm, head.width_mm, safety_d); +} + +template +std::optional search_widening_path(Ex policy, + const SupportableMesh &sm, + const Vec3d &jp, + const Vec3d &dir, + double radius, + double new_radius) +{ + double w = radius + 2 * sm.cfg.head_back_radius_mm; + double stopval = w + jp.z() - ground_level(sm); + Optimizer solver(get_criteria(sm.cfg).stop_score(stopval)); + + auto [polar, azimuth] = dir_to_spheric(dir); + + double fallback_ratio = radius / sm.cfg.head_back_radius_mm; + + auto oresult = solver.to_max().optimize( + [&policy, &sm, jp, radius, new_radius](const opt::Input<3> &input) { + auto &[plr, azm, t] = input; + + auto d = spheric_to_dir(plr, azm).normalized(); + + auto sd = new_radius * sm.cfg.safety_distance_mm / + sm.cfg.head_back_radius_mm; + + double ret = pinhead_mesh_hit(policy, sm.emesh, jp, d, radius, + new_radius, t, sd) + .distance(); + + Beam beam{jp + t * d, d, new_radius}; + double down = beam_mesh_hit(policy, sm.emesh, beam, sd).distance(); + + if (ret > t && std::isinf(down)) + ret += jp.z() - ground_level(sm); + + return ret; + }, + initvals({polar, azimuth, w}), // start with what we have + bounds({ + {PI - sm.cfg.bridge_slope, PI}, // Must not exceed the slope limit + {-PI, PI}, // azimuth can be a full search + {radius + sm.cfg.head_back_radius_mm, + fallback_ratio * sm.cfg.max_bridge_length_mm} + })); + + if (oresult.score >= stopval) { + polar = std::get<0>(oresult.optimum); + azimuth = std::get<1>(oresult.optimum); + double t = std::get<2>(oresult.optimum); + Vec3d endp = jp + t * spheric_to_dir(polar, azimuth); + + return DiffBridge(jp, endp, radius, sm.cfg.head_back_radius_mm); + } + + return {}; +} + +// This is a proxy function for pillar creation which will mind the gap +// between the pad and the model bottom in zero elevation mode. +// 'pinhead_junctionpt' is the starting junction point which needs to be +// routed down. sourcedir is the allowed direction of an optional bridge +// between the jp junction and the final pillar. +template +std::pair create_ground_pillar( + Ex policy, + SupportTreeBuilder &builder, + const SupportableMesh &sm, + const Vec3d &pinhead_junctionpt, + const Vec3d &sourcedir, + double radius, + double end_radius, + long head_id = SupportTreeNode::ID_UNSET) +{ + Vec3d jp = pinhead_junctionpt, endp = jp, dir = sourcedir; + long pillar_id = SupportTreeNode::ID_UNSET; + bool can_add_base = false, non_head = false; + + double gndlvl = 0.; // The Z level where pedestals should be + double jp_gnd = 0.; // The lowest Z where a junction center can be + double gap_dist = 0.; // The gap distance between the model and the pad + + double r2 = radius + (end_radius - radius) / (jp.z() - ground_level(sm)); + + auto to_floor = [&gndlvl](const Vec3d &p) { return Vec3d{p.x(), p.y(), gndlvl}; }; + + auto eval_limits = [&sm, &radius, &can_add_base, &gndlvl, &gap_dist, &jp_gnd] + (bool base_en = true) + { + can_add_base = base_en && radius >= sm.cfg.head_back_radius_mm; + double base_r = can_add_base ? sm.cfg.base_radius_mm : 0.; + gndlvl = ground_level(sm); + if (!can_add_base) gndlvl -= sm.pad_cfg.wall_thickness_mm; + jp_gnd = gndlvl + (can_add_base ? 0. : sm.cfg.head_back_radius_mm); + gap_dist = sm.cfg.pillar_base_safety_distance_mm + base_r + EPSILON; + }; + + eval_limits(); + + // We are dealing with a mini pillar that's potentially too long + if (radius < sm.cfg.head_back_radius_mm && jp.z() - gndlvl > 20 * radius) + { + std::optional diffbr = + search_widening_path(policy, sm, jp, dir, radius, + sm.cfg.head_back_radius_mm); + + if (diffbr && diffbr->endp.z() > jp_gnd) { + auto &br = builder.add_diffbridge(*diffbr); + if (head_id >= 0) builder.head(head_id).bridge_id = br.id; + endp = diffbr->endp; + radius = diffbr->end_r; + builder.add_junction(endp, radius); + non_head = true; + dir = diffbr->get_dir(); + eval_limits(); + } else return {false, pillar_id}; + } + + if (sm.cfg.object_elevation_mm < EPSILON) + { + // get a suitable direction for the corrector bridge. It is the + // original sourcedir's azimuth but the polar angle is saturated to the + // configured bridge slope. + auto [polar, azimuth] = dir_to_spheric(dir); + polar = PI - sm.cfg.bridge_slope; + Vec3d d = spheric_to_dir(polar, azimuth).normalized(); + auto sd = radius * sm.cfg.safety_distance_mm / sm.cfg.head_back_radius_mm; + double t = beam_mesh_hit(policy, sm.emesh, Beam{endp, d, radius, r2}, sd).distance(); + double tmax = std::min(sm.cfg.max_bridge_length_mm, t); + t = 0.; + + double zd = endp.z() - jp_gnd; + double tmax2 = zd / std::sqrt(1 - sm.cfg.bridge_slope * sm.cfg.bridge_slope); + tmax = std::min(tmax, tmax2); + + Vec3d nexp = endp; + double dlast = 0.; + while (((dlast = std::sqrt(sm.emesh.squared_distance(to_floor(nexp)))) < gap_dist || + !std::isinf(beam_mesh_hit(policy, sm.emesh, Beam{nexp, DOWN, radius, r2}, sd).distance())) && + t < tmax) + { + t += radius; + nexp = endp + t * d; + } + + if (dlast < gap_dist && can_add_base) { + nexp = endp; + t = 0.; + can_add_base = false; + eval_limits(can_add_base); + + zd = endp.z() - jp_gnd; + tmax2 = zd / std::sqrt(1 - sm.cfg.bridge_slope * sm.cfg.bridge_slope); + tmax = std::min(tmax, tmax2); + + while (((dlast = std::sqrt(sm.emesh.squared_distance(to_floor(nexp)))) < gap_dist || + !std::isinf(beam_mesh_hit(policy, sm.emesh, Beam{nexp, DOWN, radius}, sd).distance())) && t < tmax) { + t += radius; + nexp = endp + t * d; + } + } + + // Could not find a path to avoid the pad gap + if (dlast < gap_dist) return {false, pillar_id}; + + if (t > 0.) { // Need to make additional bridge + const Bridge& br = builder.add_bridge(endp, nexp, radius); + if (head_id >= 0) builder.head(head_id).bridge_id = br.id; + + builder.add_junction(nexp, radius); + endp = nexp; + non_head = true; + } + } + + Vec3d gp = to_floor(endp); + double h = endp.z() - gp.z(); + + pillar_id = head_id >= 0 && !non_head ? builder.add_pillar(head_id, h) : + builder.add_pillar(gp, h, radius, end_radius); + + if (can_add_base) + builder.add_pillar_base(pillar_id, sm.cfg.base_height_mm, + sm.cfg.base_radius_mm); + + return {true, pillar_id}; +} + +inline double distance(const SupportPoint &a, const SupportPoint &b) +{ + return (a.pos - b.pos).norm(); +} + +template +std::vector non_duplicate_suppt_indices(const PtIndex &index, + const SupportPoints &suppts, + double eps) +{ + std::vector to_remove(suppts.size(), false); + + for (size_t i = 0; i < suppts.size(); ++i) { + size_t closest_idx = + find_closest_point(index, suppts[i].pos, + [&i, &to_remove](size_t i_closest) { + return i_closest != i && + !to_remove[i_closest]; + }); + + if ((suppts[i].pos - suppts[closest_idx].pos).norm() < eps) + to_remove[i] = true; + } + + auto ret = reserve_vector(suppts.size()); + for (size_t i = 0; i < to_remove.size(); i++) + if (!to_remove[i]) + ret.emplace_back(i); + + return ret; +} + +template +bool optimize_pinhead_placement(Ex policy, + const SupportableMesh &m, + Head &head) +{ + Vec3d n = get_normal(m.emesh, head.pos); + + // for all normals the spherical coordinates are generated and + // the polar angle is saturated to 45 degrees from the bottom then + // converted back to standard coordinates to get the new normal. + // Then a simple quaternion is created from the two normals + // (Quaternion::FromTwoVectors) and the rotation is applied to the + // pinhead. + + auto [polar, azimuth] = dir_to_spheric(n); + + double back_r = head.r_back_mm; + + // skip if the tilt is not sane + if (polar < PI - m.cfg.normal_cutoff_angle) return false; + + // We saturate the polar angle to 3pi/4 + polar = std::max(polar, PI - m.cfg.bridge_slope); + + // save the head (pinpoint) position + Vec3d hp = head.pos; + + double lmin = m.cfg.head_width_mm, lmax = lmin; + + if (back_r < m.cfg.head_back_radius_mm) { + lmin = 0., lmax = m.cfg.head_penetration_mm; + } + + // The distance needed for a pinhead to not collide with model. + double w = lmin + 2 * back_r + 2 * m.cfg.head_front_radius_mm - + m.cfg.head_penetration_mm; + + double pin_r = head.r_pin_mm; + + // Reassemble the now corrected normal + auto nn = spheric_to_dir(polar, azimuth).normalized(); + + double sd = back_r * m.cfg.safety_distance_mm / + m.cfg.head_back_radius_mm; + + // check available distance + Hit t = pinhead_mesh_hit(policy, m.emesh, hp, nn, pin_r, back_r, w, + sd); + + if (t.distance() < w) { + // Let's try to optimize this angle, there might be a + // viable normal that doesn't collide with the model + // geometry and its very close to the default. + + Optimizer solver(get_criteria(m.cfg).stop_score(w).max_iterations(100)); + solver.seed(0); // we want deterministic behavior + + auto oresult = solver.to_max().optimize( + [&m, pin_r, back_r, hp, sd, policy](const opt::Input<3> &input) { + auto &[plr, azm, l] = input; + + auto dir = spheric_to_dir(plr, azm).normalized(); + + return pinhead_mesh_hit(policy, m.emesh, hp, dir, pin_r, + back_r, l, sd).distance(); + }, + initvals({polar, azimuth, + (lmin + lmax) / 2.}), // start with what we have + bounds({{PI - m.cfg.bridge_slope, PI}, // Must not exceed the slope limit + {-PI, PI}, // azimuth can be a full search + {lmin, lmax}})); + + if(oresult.score > w) { + polar = std::get<0>(oresult.optimum); + azimuth = std::get<1>(oresult.optimum); + nn = spheric_to_dir(polar, azimuth).normalized(); + lmin = std::get<2>(oresult.optimum); + t = AABBMesh::hit_result(oresult.score); + } + } + + bool ret = false; + if (t.distance() > w && hp.z() + w * nn.z() >= ground_level(m)) { + head.dir = nn; + head.width_mm = lmin; + head.r_back_mm = back_r; + + ret = true; + } else if (back_r > m.cfg.head_fallback_radius_mm) { + head.r_back_mm = m.cfg.head_fallback_radius_mm; + ret = optimize_pinhead_placement(policy, m, head); + } + + return ret; +} + +template +std::optional calculate_pinhead_placement(Ex policy, + const SupportableMesh &sm, + size_t suppt_idx) +{ + if (suppt_idx >= sm.pts.size()) + return {}; + + const SupportPoint &sp = sm.pts[suppt_idx]; + Head head{ + sm.cfg.head_back_radius_mm, + sp.head_front_radius, + 0., + sm.cfg.head_penetration_mm, + Vec3d::Zero(), // dir + sp.pos.cast() // displacement + }; + + if (optimize_pinhead_placement(policy, sm, head)) { + head.id = suppt_idx; + + return head; + } + + return {}; +} + +template +std::pair connect_to_ground(Ex policy, + SupportTreeBuilder &builder, + const SupportableMesh &sm, + const Junction &j, + const Vec3d &dir, + double end_r) +{ + auto hjp = j.pos; + double r = j.r; + auto sd = r * sm.cfg.safety_distance_mm / sm.cfg.head_back_radius_mm; + double r2 = j.r + (end_r - j.r) / (j.pos.z() - ground_level(sm)); + + double t = beam_mesh_hit(policy, sm.emesh, Beam{hjp, dir, r, r2}, sd).distance(); + double d = 0, tdown = 0; + t = std::min(t, sm.cfg.max_bridge_length_mm * r / sm.cfg.head_back_radius_mm); + + while (d < t && + !std::isinf(tdown = beam_mesh_hit(policy, sm.emesh, + Beam{hjp + d * dir, DOWN, r, r2}, sd) + .distance())) { + d += r; + } + + if(!std::isinf(tdown)) + return {false, SupportTreeNode::ID_UNSET}; + + Vec3d endp = hjp + d * dir; + auto ret = create_ground_pillar(policy, builder, sm, endp, dir, r, end_r); + + if (ret.second >= 0) { + builder.add_bridge(hjp, endp, r); + builder.add_junction(endp, r); + } + + return ret; +} + +template +std::pair search_ground_route(Ex policy, + SupportTreeBuilder &builder, + const SupportableMesh &sm, + const Junction &j, + double end_radius, + const Vec3d &init_dir = DOWN) +{ + double downdst = j.pos.z() - ground_level(sm); + + auto res = connect_to_ground(policy, builder, sm, j, init_dir, end_radius); + if (res.first) + return res; + + // Optimize bridge direction: + // Straight path failed so we will try to search for a suitable + // direction out of the cavity. + auto [polar, azimuth] = dir_to_spheric(init_dir); + + Optimizer solver(get_criteria(sm.cfg).stop_score(1e6)); + solver.seed(0); // we want deterministic behavior + + auto sd = j.r * sm.cfg.safety_distance_mm / sm.cfg.head_back_radius_mm; + auto oresult = solver.to_max().optimize( + [&j, sd, &policy, &sm, &downdst, &end_radius](const opt::Input<2> &input) { + auto &[plr, azm] = input; + Vec3d n = spheric_to_dir(plr, azm).normalized(); + Beam beam{Ball{j.pos, j.r}, Ball{j.pos + downdst * n, end_radius}}; + return beam_mesh_hit(policy, sm.emesh, beam, sd).distance(); + }, + initvals({polar, azimuth}), // let's start with what we have + bounds({ {PI - sm.cfg.bridge_slope, PI}, {-PI, PI} }) + ); + + Vec3d bridgedir = spheric_to_dir(oresult.optimum).normalized(); + + return connect_to_ground(policy, builder, sm, j, bridgedir, end_radius); +} + +template +bool optimize_anchor_placement(Ex policy, + const SupportableMesh &sm, + const Junction &from, + Anchor &anchor) +{ + Vec3d n = get_normal(sm.emesh, anchor.pos); + + auto [polar, azimuth] = dir_to_spheric(n); + + // Saturate the polar angle to 3pi/4 + polar = std::min(polar, sm.cfg.bridge_slope); + + double lmin = 0; + double lmax = std::min(sm.cfg.head_width_mm, + distance(from.pos, anchor.pos) - 2 * from.r); + + double sd = anchor.r_back_mm * sm.cfg.safety_distance_mm / + sm.cfg.head_back_radius_mm; + + Optimizer solver(get_criteria(sm.cfg) + .stop_score(anchor.fullwidth()) + .max_iterations(100)); + + solver.seed(0); // deterministic behavior + + auto oresult = solver.to_max().optimize( + [&sm, &anchor, sd, policy](const opt::Input<3> &input) { + auto &[plr, azm, l] = input; + + auto dir = spheric_to_dir(plr, azm).normalized(); + + anchor.width_mm = l; + anchor.dir = dir; + + return pinhead_mesh_hit(policy, sm.emesh, anchor, sd) + .distance(); + }, + initvals({polar, azimuth, (lmin + lmax) / 2.}), + bounds({{0., sm.cfg.bridge_slope}, // Must not exceed the slope limit + {-PI, PI}, // azimuth can be a full search + {lmin, lmax}})); + + polar = std::get<0>(oresult.optimum); + azimuth = std::get<1>(oresult.optimum); + anchor.dir = spheric_to_dir(polar, azimuth).normalized(); + anchor.width_mm = std::get<2>(oresult.optimum); + + if (oresult.score < anchor.fullwidth()) { + // Unsuccesful search, the anchor does not fit into its intended space. + return false; + } + + return true; +} + +template +std::optional calculate_anchor_placement(Ex policy, + const SupportableMesh &sm, + const Junction &from, + const Vec3d &to_hint) +{ + double back_r = from.r; + double pin_r = sm.cfg.head_front_radius_mm; + double penetr = sm.cfg.head_penetration_mm; + double hwidth = sm.cfg.head_width_mm; + Vec3d bridgedir = dirv(from.pos, to_hint); + Vec3d anchordir = -bridgedir; + + std::optional ret; + + Anchor anchor(back_r, pin_r, hwidth, penetr, anchordir, to_hint); + + if (optimize_anchor_placement(policy, sm, from, anchor)) { + ret = anchor; + } else if (anchor.r_back_mm = sm.cfg.head_fallback_radius_mm; + optimize_anchor_placement(policy, sm, from, anchor)) { + // Retrying with the fallback strut radius as a last resort. + ret = anchor; + } + + return anchor; +} + +}} // namespace Slic3r::sla + +#endif // SLASUPPORTTREEUTILS_H diff --git a/src/libslic3r/SLAPrint.cpp b/src/libslic3r/SLAPrint.cpp index c4494c3c2..87da7af12 100644 --- a/src/libslic3r/SLAPrint.cpp +++ b/src/libslic3r/SLAPrint.cpp @@ -55,14 +55,7 @@ sla::SupportTreeConfig make_support_cfg(const SLAPrintObjectConfig& c) scfg.bridge_slope = c.support_critical_angle.getFloat() * PI / 180.0 ; scfg.max_bridge_length_mm = c.support_max_bridge_length.getFloat(); scfg.max_pillar_link_distance_mm = c.support_max_pillar_link_distance.getFloat(); - switch(c.support_pillar_connection_mode.getInt()) { - case slapcmZigZag: - scfg.pillar_connection_mode = sla::PillarConnectionMode::zigzag; break; - case slapcmCross: - scfg.pillar_connection_mode = sla::PillarConnectionMode::cross; break; - case slapcmDynamic: - scfg.pillar_connection_mode = sla::PillarConnectionMode::dynamic; break; - } + scfg.pillar_connection_mode = c.support_pillar_connection_mode.value; scfg.ground_facing_only = c.support_buildplate_only.getBool(); scfg.pillar_widening_factor = c.support_pillar_widening_factor.getFloat(); scfg.base_radius_mm = 0.5*c.support_base_diameter.getFloat(); @@ -844,6 +837,7 @@ bool SLAPrintObject::invalidate_state_by_config_options(const std::vector& SLAPrintObject::get_support_points() const { - return m_supportdata? m_supportdata->pts : EMPTY_SUPPORT_POINTS; + return m_supportdata? m_supportdata->input.pts : EMPTY_SUPPORT_POINTS; } const std::vector &SLAPrintObject::get_support_slices() const diff --git a/src/libslic3r/SLAPrint.hpp b/src/libslic3r/SLAPrint.hpp index 56309b30b..528d3c28b 100644 --- a/src/libslic3r/SLAPrint.hpp +++ b/src/libslic3r/SLAPrint.hpp @@ -309,31 +309,24 @@ private: // Caching the transformed (m_trafo) raw mesh of the object mutable CachedObject m_transformed_rmesh; - class SupportData : public sla::SupportableMesh + struct SupportData { - public: - sla::SupportTree::UPtr support_tree_ptr; // the supports + sla::SupportableMesh input; // the input std::vector support_slices; // sliced supports - TriangleMesh tree_mesh, pad_mesh, full_mesh; + TriangleMesh tree_mesh, pad_mesh, full_mesh; // cached artifacts inline SupportData(const TriangleMesh &t) - : sla::SupportableMesh{t.its, {}, {}} + : input{t.its, {}, {}} {} - sla::SupportTree::UPtr &create_support_tree(const sla::JobController &ctl) + void create_support_tree(const sla::JobController &ctl) { - support_tree_ptr = sla::SupportTree::create(*this, ctl); - tree_mesh = TriangleMesh{support_tree_ptr->retrieve_mesh(sla::MeshType::Support)}; - return support_tree_ptr; + tree_mesh = TriangleMesh{sla::create_support_tree(input, ctl)}; } - void create_pad(const ExPolygons &blueprint, const sla::PadConfig &pcfg) + void create_pad(const sla::JobController &ctl) { - if (!support_tree_ptr) - return; - - support_tree_ptr->add_pad(blueprint, pcfg); - pad_mesh = TriangleMesh{support_tree_ptr->retrieve_mesh(sla::MeshType::Pad)}; + pad_mesh = TriangleMesh{sla::create_pad(input, tree_mesh.its, ctl)}; } }; diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp index c3ec3e399..25779a662 100644 --- a/src/libslic3r/SLAPrintSteps.cpp +++ b/src/libslic3r/SLAPrintSteps.cpp @@ -8,7 +8,7 @@ // Need the cylinder method for the the drainholes in hollowing step #include -#include +#include #include #include @@ -546,12 +546,13 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po) std::vector interior_slices = slice_mesh_ex(interiormesh, slice_grid, params, thr); - sla::ccr::for_each(size_t(0), interior_slices.size(), - [&po, &interior_slices] (size_t i) { - const ExPolygons &slice = interior_slices[i]; - po.m_model_slices[i] = - diff_ex(po.m_model_slices[i], slice); - }); + execution::for_each( + ex_tbb, size_t(0), interior_slices.size(), + [&po, &interior_slices](size_t i) { + const ExPolygons &slice = interior_slices[i]; + po.m_model_slices[i] = diff_ex(po.m_model_slices[i], slice); + }, + execution::max_concurrency(ex_tbb)); } auto mit = slindex_it; @@ -621,16 +622,17 @@ void SLAPrint::Steps::support_points(SLAPrintObject &po) // Construction of this object does the calculation. throw_if_canceled(); sla::SupportPointGenerator auto_supports( - po.m_supportdata->emesh, po.get_model_slices(), heights, config, - [this]() { throw_if_canceled(); }, statuscb); + po.m_supportdata->input.emesh, po.get_model_slices(), + heights, config, [this]() { throw_if_canceled(); }, statuscb); // Now let's extract the result. const std::vector& points = auto_supports.output(); throw_if_canceled(); - po.m_supportdata->pts = points; + po.m_supportdata->input.pts = points; - BOOST_LOG_TRIVIAL(debug) << "Automatic support points: " - << po.m_supportdata->pts.size(); + BOOST_LOG_TRIVIAL(debug) + << "Automatic support points: " + << po.m_supportdata->input.pts.size(); // Using RELOAD_SLA_SUPPORT_POINTS to tell the Plater to pass // the update status to GLGizmoSlaSupports @@ -639,7 +641,7 @@ void SLAPrint::Steps::support_points(SLAPrintObject &po) } else { // There are either some points on the front-end, or the user // removed them on purpose. No calculation will be done. - po.m_supportdata->pts = po.transformed_support_points(); + po.m_supportdata->input.pts = po.transformed_support_points(); } } @@ -647,19 +649,22 @@ void SLAPrint::Steps::support_tree(SLAPrintObject &po) { if(!po.m_supportdata) return; - sla::PadConfig pcfg = make_pad_cfg(po.m_config); +// sla::PadConfig pcfg = make_pad_cfg(po.m_config); - if (pcfg.embed_object) - po.m_supportdata->emesh.ground_level_offset(pcfg.wall_thickness_mm); +// if (pcfg.embed_object) +// po.m_supportdata->emesh.ground_level_offset(pcfg.wall_thickness_mm); // If the zero elevation mode is engaged, we have to filter out all the // points that are on the bottom of the object if (is_zero_elevation(po.config())) { - remove_bottom_points(po.m_supportdata->pts, - float(po.m_supportdata->emesh.ground_level() + EPSILON)); + remove_bottom_points(po.m_supportdata->input.pts, + float( + po.m_supportdata->input.emesh.ground_level() + + EPSILON)); } - po.m_supportdata->cfg = make_support_cfg(po.m_config); + po.m_supportdata->input.cfg = make_support_cfg(po.m_config); + po.m_supportdata->input.pad_cfg = make_pad_cfg(po.m_config); // po.m_supportdata->emesh.load_holes(po.transformed_drainhole_points()); // scaling for the sub operations @@ -689,7 +694,7 @@ void SLAPrint::Steps::support_tree(SLAPrintObject &po) report_status(-1, L("Visualizing supports")); BOOST_LOG_TRIVIAL(debug) << "Processed support point count " - << po.m_supportdata->pts.size(); + << po.m_supportdata->input.pts.size(); // Check the mesh for later troubleshooting. if(po.support_mesh().empty()) @@ -705,31 +710,34 @@ void SLAPrint::Steps::generate_pad(SLAPrintObject &po) { if(po.m_config.pad_enable.getBool()) { // Get the distilled pad configuration from the config + // (Again, despite it was retrieved in the previous step. Note that + // on a param change event, the previous step might not be executed + // depending on the specific parameter that has changed). sla::PadConfig pcfg = make_pad_cfg(po.m_config); + po.m_supportdata->input.pad_cfg = pcfg; - ExPolygons bp; // This will store the base plate of the pad. - double pad_h = pcfg.full_height(); - const TriangleMesh &trmesh = po.transformed_mesh(); +// if (!po.m_config.supports_enable.getBool() || pcfg.embed_object) { +// // No support (thus no elevation) or zero elevation mode +// // we sometimes call it "builtin pad" is enabled so we will +// // get a sample from the bottom of the mesh and use it for pad +// // creation. +// sla::pad_blueprint(trmesh.its, bp, float(pad_h), +// float(po.m_config.layer_height.getFloat()), +// [this](){ throw_if_canceled(); }); +// } - if (!po.m_config.supports_enable.getBool() || pcfg.embed_object) { - // No support (thus no elevation) or zero elevation mode - // we sometimes call it "builtin pad" is enabled so we will - // get a sample from the bottom of the mesh and use it for pad - // creation. - sla::pad_blueprint(trmesh.its, bp, float(pad_h), - float(po.m_config.layer_height.getFloat()), - [this](){ throw_if_canceled(); }); - } + sla::JobController ctl; + ctl.stopcondition = [this]() { return canceled(); }; + ctl.cancelfn = [this]() { throw_if_canceled(); }; + po.m_supportdata->create_pad(ctl); - po.m_supportdata->create_pad(bp, pcfg); - - if (!validate_pad(po.m_supportdata->support_tree_ptr->retrieve_mesh(sla::MeshType::Pad), pcfg)) + if (!validate_pad(po.m_supportdata->pad_mesh.its, pcfg)) throw Slic3r::SlicingError( L("No pad can be generated for this model with the " "current configuration")); - } else if(po.m_supportdata && po.m_supportdata->support_tree_ptr) { - po.m_supportdata->support_tree_ptr->remove_pad(); + } else if(po.m_supportdata) { + po.m_supportdata->pad_mesh = {}; } throw_if_canceled(); @@ -748,13 +756,18 @@ void SLAPrint::Steps::slice_supports(SLAPrintObject &po) { if (!po.m_config.supports_enable.getBool() && !po.m_config.pad_enable.getBool()) return; - if(sd && sd->support_tree_ptr) { + if(sd) { auto heights = reserve_vector(po.m_slice_index.size()); for(auto& rec : po.m_slice_index) heights.emplace_back(rec.slice_level()); - sd->support_slices = sd->support_tree_ptr->slice( - heights, float(po.config().slice_closing_radius.value)); + sla::JobController ctl; + ctl.stopcondition = [this]() { return canceled(); }; + ctl.cancelfn = [this]() { throw_if_canceled(); }; + + sd->support_slices = + sla::slice(sd->tree_mesh.its, sd->pad_mesh.its, heights, + float(po.config().slice_closing_radius.value), ctl); } for (size_t i = 0; i < sd->support_slices.size() && i < po.m_slice_index.size(); ++i) @@ -910,8 +923,8 @@ void SLAPrint::Steps::merge_slices_and_eval_stats() { const double delta_fade_time = (init_exp_time - exp_time) / (fade_layers_cnt + 1); double fade_layer_time = init_exp_time; - sla::ccr::SpinningMutex mutex; - using Lock = std::lock_guard; + execution::SpinningMutex mutex; + using Lock = std::lock_guard; // Going to parallel: auto printlayerfn = [this, @@ -1027,7 +1040,8 @@ void SLAPrint::Steps::merge_slices_and_eval_stats() { // sequential version for debugging: // for(size_t i = 0; i < m_printer_input.size(); ++i) printlayerfn(i); - sla::ccr::for_each(size_t(0), printer_input.size(), printlayerfn); + execution::for_each(ex_tbb, size_t(0), printer_input.size(), printlayerfn, + execution::max_concurrency(ex_tbb)); auto SCALING2 = SCALING_FACTOR * SCALING_FACTOR; print_statistics.support_used_material = supports_volume * SCALING2; @@ -1066,8 +1080,7 @@ void SLAPrint::Steps::rasterize() double increment = (slot * sd) / m_print->m_printer_input.size(); double dstatus = current_status(); - sla::ccr::SpinningMutex slck; - using Lock = std::lock_guard; + execution::SpinningMutex slck; // procedure to process one height level. This will run in parallel auto lvlfn = @@ -1082,7 +1095,7 @@ void SLAPrint::Steps::rasterize() // Status indication guarded with the spinlock { - Lock lck(slck); + std::lock_guard lck(slck); dstatus += increment; double st = std::round(dstatus); if(st > pst) { diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 3ee8acfc3..6f71930ab 100644 --- a/src/libslic3r/libslic3r.h +++ b/src/libslic3r/libslic3r.h @@ -126,6 +126,15 @@ inline void append(std::vector& dest, std::vector&& src) src.shrink_to_fit(); } +template // Arbitrary allocator can be used +void clear_and_shrink(std::vector& vec) +{ + // shrink_to_fit does not garantee the release of memory nor does it clear() + std::vector tmp; + vec.swap(tmp); + assert(vec.capacity() == 0); +} + // Append the source in reverse. template inline void append_reversed(std::vector& dest, const std::vector& src) diff --git a/src/slic3r/GUI/MeshUtils.cpp b/src/slic3r/GUI/MeshUtils.cpp index 7782b8d75..3b26586e7 100644 --- a/src/slic3r/GUI/MeshUtils.cpp +++ b/src/slic3r/GUI/MeshUtils.cpp @@ -1,378 +1,380 @@ -#include "MeshUtils.hpp" - -#include "libslic3r/Tesselate.hpp" -#include "libslic3r/TriangleMesh.hpp" -#include "libslic3r/TriangleMeshSlicer.hpp" -#include "libslic3r/ClipperUtils.hpp" -#include "libslic3r/Model.hpp" - -#if ENABLE_LEGACY_OPENGL_REMOVAL -#include "slic3r/GUI/GUI_App.hpp" -#endif // ENABLE_LEGACY_OPENGL_REMOVAL -#include "slic3r/GUI/Camera.hpp" -#if ENABLE_GL_SHADERS_ATTRIBUTES -#include "slic3r/GUI/Plater.hpp" -#endif // ENABLE_GL_SHADERS_ATTRIBUTES - -#include - -#include - - -namespace Slic3r { -namespace GUI { - -void MeshClipper::set_plane(const ClippingPlane& plane) -{ - if (m_plane != plane) { - m_plane = plane; - m_triangles_valid = false; - } -} - - -void MeshClipper::set_limiting_plane(const ClippingPlane& plane) -{ - if (m_limiting_plane != plane) { - m_limiting_plane = plane; - m_triangles_valid = false; - } -} - - - -void MeshClipper::set_mesh(const TriangleMesh& mesh) -{ - if (m_mesh != &mesh) { - m_mesh = &mesh; - m_triangles_valid = false; - m_triangles2d.resize(0); - } -} - -void MeshClipper::set_negative_mesh(const TriangleMesh& mesh) -{ - if (m_negative_mesh != &mesh) { - m_negative_mesh = &mesh; - m_triangles_valid = false; - m_triangles2d.resize(0); - } -} - - - -void MeshClipper::set_transformation(const Geometry::Transformation& trafo) -{ - if (! m_trafo.get_matrix().isApprox(trafo.get_matrix())) { - m_trafo = trafo; - m_triangles_valid = false; - m_triangles2d.resize(0); - } -} - -#if ENABLE_LEGACY_OPENGL_REMOVAL -void MeshClipper::render_cut(const ColorRGBA& color) -#else -void MeshClipper::render_cut() -#endif // ENABLE_LEGACY_OPENGL_REMOVAL -{ - if (! m_triangles_valid) - recalculate_triangles(); - -#if ENABLE_LEGACY_OPENGL_REMOVAL - if (m_model.vertices_count() == 0 || m_model.indices_count() == 0) - return; - - GLShaderProgram* curr_shader = wxGetApp().get_current_shader(); - if (curr_shader != nullptr) - curr_shader->stop_using(); - - GLShaderProgram* shader = wxGetApp().get_shader("flat"); - if (shader != nullptr) { - shader->start_using(); -#if ENABLE_GL_SHADERS_ATTRIBUTES - const Camera& camera = wxGetApp().plater()->get_camera(); - shader->set_uniform("view_model_matrix", camera.get_view_matrix()); - shader->set_uniform("projection_matrix", camera.get_projection_matrix()); -#endif // ENABLE_GL_SHADERS_ATTRIBUTES - m_model.set_color(color); - m_model.render(); - shader->stop_using(); - } - - if (curr_shader != nullptr) - curr_shader->start_using(); -#else - if (m_vertex_array.has_VBOs()) - m_vertex_array.render(); -#endif // ENABLE_LEGACY_OPENGL_REMOVAL -} - - - -void MeshClipper::recalculate_triangles() -{ -#if ENABLE_WORLD_COORDINATE - const Transform3f instance_matrix_no_translation_no_scaling = m_trafo.get_rotation_matrix().cast(); -#else - const Transform3f& instance_matrix_no_translation_no_scaling = m_trafo.get_matrix(true,false,true).cast(); -#endif // ENABLE_WORLD_COORDINATE - // Calculate clipping plane normal in mesh coordinates. - const Vec3f up_noscale = instance_matrix_no_translation_no_scaling.inverse() * m_plane.get_normal().cast(); - const Vec3d up = up_noscale.cast().cwiseProduct(m_trafo.get_scaling_factor()); - // Calculate distance from mesh origin to the clipping plane (in mesh coordinates). - const float height_mesh = m_plane.distance(m_trafo.get_offset()) * (up_noscale.norm()/up.norm()); - - // Now do the cutting - MeshSlicingParams slicing_params; - slicing_params.trafo.rotate(Eigen::Quaternion::FromTwoVectors(up, Vec3d::UnitZ())); - - ExPolygons expolys = union_ex(slice_mesh(m_mesh->its, height_mesh, slicing_params)); - - if (m_negative_mesh && !m_negative_mesh->empty()) { - const ExPolygons neg_expolys = union_ex(slice_mesh(m_negative_mesh->its, height_mesh, slicing_params)); - expolys = diff_ex(expolys, neg_expolys); - } - - // Triangulate and rotate the cut into world coords: - Eigen::Quaterniond q; - q.setFromTwoVectors(Vec3d::UnitZ(), up); - Transform3d tr = Transform3d::Identity(); - tr.rotate(q); - tr = m_trafo.get_matrix() * tr; - - if (m_limiting_plane != ClippingPlane::ClipsNothing()) - { - // Now remove whatever ended up below the limiting plane (e.g. sinking objects). - // First transform the limiting plane from world to mesh coords. - // Note that inverse of tr transforms the plane from world to horizontal. - const Vec3d normal_old = m_limiting_plane.get_normal().normalized(); - const Vec3d normal_new = (tr.matrix().block<3,3>(0,0).transpose() * normal_old).normalized(); - - // normal_new should now be the plane normal in mesh coords. To find the offset, - // transform a point and set offset so it belongs to the transformed plane. - Vec3d pt = Vec3d::Zero(); - const double plane_offset = m_limiting_plane.get_data()[3]; - if (std::abs(normal_old.z()) > 0.5) // normal is normalized, at least one of the coords if larger than sqrt(3)/3 = 0.57 - pt.z() = - plane_offset / normal_old.z(); - else if (std::abs(normal_old.y()) > 0.5) - pt.y() = - plane_offset / normal_old.y(); - else - pt.x() = - plane_offset / normal_old.x(); - pt = tr.inverse() * pt; - const double offset = -(normal_new.dot(pt)); - - if (std::abs(normal_old.dot(m_plane.get_normal().normalized())) > 0.99) { - // The cuts are parallel, show all or nothing. - if (normal_old.dot(m_plane.get_normal().normalized()) < 0.0 && offset < height_mesh) - expolys.clear(); - } else { - // The cut is a horizontal plane defined by z=height_mesh. - // ax+by+e=0 is the line of intersection with the limiting plane. - // Normalized so a^2 + b^2 = 1. - const double len = std::hypot(normal_new.x(), normal_new.y()); - if (len == 0.) - return; - const double a = normal_new.x() / len; - const double b = normal_new.y() / len; - const double e = (normal_new.z() * height_mesh + offset) / len; - - // We need a half-plane to limit the cut. Get angle of the intersecting line. - double angle = (b != 0.0) ? std::atan(-a / b) : ((a < 0.0) ? -0.5 * M_PI : 0.5 * M_PI); - if (b > 0) // select correct half-plane - angle += M_PI; - - // We'll take a big rectangle above x-axis and rotate and translate - // it so it lies on our line. This will be the figure to subtract - // from the cut. The coordinates must not overflow after the transform, - // make the rectangle a bit smaller. - const coord_t size = (std::numeric_limits::max() - scale_(std::max(std::abs(e*a), std::abs(e*b)))) / 4; - Polygons ep {Polygon({Point(-size, 0), Point(size, 0), Point(size, 2*size), Point(-size, 2*size)})}; - ep.front().rotate(angle); - ep.front().translate(scale_(-e * a), scale_(-e * b)); - expolys = diff_ex(expolys, ep); - } - } - - m_triangles2d = triangulate_expolygons_2f(expolys, m_trafo.get_matrix().matrix().determinant() < 0.); - - tr.pretranslate(0.001 * m_plane.get_normal().normalized()); // to avoid z-fighting - -#if ENABLE_LEGACY_OPENGL_REMOVAL - m_model.reset(); - - GLModel::Geometry init_data; - init_data.format = { GLModel::Geometry::EPrimitiveType::Triangles, GLModel::Geometry::EVertexLayout::P3N3 }; - init_data.reserve_vertices(m_triangles2d.size()); - init_data.reserve_indices(m_triangles2d.size()); - - // vertices + indices - for (auto it = m_triangles2d.cbegin(); it != m_triangles2d.cend(); it = it + 3) { - init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 0)).x(), (*(it + 0)).y(), height_mesh)).cast(), (Vec3f)up.cast()); - init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 1)).x(), (*(it + 1)).y(), height_mesh)).cast(), (Vec3f)up.cast()); - init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 2)).x(), (*(it + 2)).y(), height_mesh)).cast(), (Vec3f)up.cast()); - const size_t idx = it - m_triangles2d.cbegin(); - init_data.add_triangle((unsigned int)idx, (unsigned int)idx + 1, (unsigned int)idx + 2); - } - - if (!init_data.is_empty()) - m_model.init_from(std::move(init_data)); -#else - m_vertex_array.release_geometry(); - for (auto it=m_triangles2d.cbegin(); it != m_triangles2d.cend(); it=it+3) { - m_vertex_array.push_geometry(tr * Vec3d((*(it+0))(0), (*(it+0))(1), height_mesh), up); - m_vertex_array.push_geometry(tr * Vec3d((*(it+1))(0), (*(it+1))(1), height_mesh), up); - m_vertex_array.push_geometry(tr * Vec3d((*(it+2))(0), (*(it+2))(1), height_mesh), up); - const size_t idx = it - m_triangles2d.cbegin(); - m_vertex_array.push_triangle(idx, idx+1, idx+2); - } - m_vertex_array.finalize_geometry(true); -#endif // ENABLE_LEGACY_OPENGL_REMOVAL - - m_triangles_valid = true; -} - - -Vec3f MeshRaycaster::get_triangle_normal(size_t facet_idx) const -{ - return m_normals[facet_idx]; -} - -void MeshRaycaster::line_from_mouse_pos(const Vec2d& mouse_pos, const Transform3d& trafo, const Camera& camera, - Vec3d& point, Vec3d& direction) const -{ - Matrix4d modelview = camera.get_view_matrix().matrix(); - Matrix4d projection= camera.get_projection_matrix().matrix(); - Vec4i viewport(camera.get_viewport().data()); - - Vec3d pt1; - Vec3d pt2; - igl::unproject(Vec3d(mouse_pos(0), viewport[3] - mouse_pos(1), 0.), - modelview, projection, viewport, pt1); - igl::unproject(Vec3d(mouse_pos(0), viewport[3] - mouse_pos(1), 1.), - modelview, projection, viewport, pt2); - - Transform3d inv = trafo.inverse(); - pt1 = inv * pt1; - pt2 = inv * pt2; - - point = pt1; - direction = pt2-pt1; -} - - -bool MeshRaycaster::unproject_on_mesh(const Vec2d& mouse_pos, const Transform3d& trafo, const Camera& camera, - Vec3f& position, Vec3f& normal, const ClippingPlane* clipping_plane, - size_t* facet_idx) const -{ - Vec3d point; - Vec3d direction; - line_from_mouse_pos(mouse_pos, trafo, camera, point, direction); - - std::vector hits = m_emesh.query_ray_hits(point, direction); - - if (hits.empty()) - return false; // no intersection found - - unsigned i = 0; - - // Remove points that are obscured or cut by the clipping plane. - // Also, remove anything below the bed (sinking objects). - for (i=0; i= SINKING_Z_THRESHOLD && - (! clipping_plane || ! clipping_plane->is_point_clipped(transformed_hit))) - break; - } - - if (i==hits.size() || (hits.size()-i) % 2 != 0) { - // All hits are either clipped, or there is an odd number of unclipped - // hits - meaning the nearest must be from inside the mesh. - return false; - } - - // Now stuff the points in the provided vector and calculate normals if asked about them: - position = hits[i].position().cast(); - normal = hits[i].normal().cast(); - - if (facet_idx) - *facet_idx = hits[i].face(); - - return true; -} - - -std::vector MeshRaycaster::get_unobscured_idxs(const Geometry::Transformation& trafo, const Camera& camera, const std::vector& points, - const ClippingPlane* clipping_plane) const -{ - std::vector out; - -#if ENABLE_WORLD_COORDINATE - const Transform3d instance_matrix_no_translation_no_scaling = trafo.get_rotation_matrix(); -#else - const Transform3d& instance_matrix_no_translation_no_scaling = trafo.get_matrix(true,false,true); -#endif // ENABLE_WORLD_COORDINATE - Vec3d direction_to_camera = -camera.get_dir_forward(); - Vec3d direction_to_camera_mesh = (instance_matrix_no_translation_no_scaling.inverse() * direction_to_camera).normalized().eval(); - direction_to_camera_mesh = direction_to_camera_mesh.cwiseProduct(trafo.get_scaling_factor()); - const Transform3d inverse_trafo = trafo.get_matrix().inverse(); - - for (size_t i=0; iis_point_clipped(pt.cast())) - continue; - - bool is_obscured = false; - // Cast a ray in the direction of the camera and look for intersection with the mesh: - std::vector hits; - // Offset the start of the ray by EPSILON to account for numerical inaccuracies. - hits = m_emesh.query_ray_hits((inverse_trafo * pt.cast() + direction_to_camera_mesh * EPSILON), - direction_to_camera_mesh); - - if (! hits.empty()) { - // If the closest hit facet normal points in the same direction as the ray, - // we are looking through the mesh and should therefore discard the point: - if (hits.front().normal().dot(direction_to_camera_mesh.cast()) > 0) - is_obscured = true; - - // Eradicate all hits that the caller wants to ignore - for (unsigned j=0; jis_point_clipped(trafo.get_matrix() * hits[j].position())) { - hits.erase(hits.begin()+j); - --j; - } - } - - // FIXME: the intersection could in theory be behind the camera, but as of now we only have camera direction. - // Also, the threshold is in mesh coordinates, not in actual dimensions. - if (! hits.empty()) - is_obscured = true; - } - if (! is_obscured) - out.push_back(i); - } - return out; -} - - -Vec3f MeshRaycaster::get_closest_point(const Vec3f& point, Vec3f* normal) const -{ - int idx = 0; - Vec3d closest_point; - m_emesh.squared_distance(point.cast(), idx, closest_point); - if (normal) - *normal = m_normals[idx]; - - return closest_point.cast(); -} - -int MeshRaycaster::get_closest_facet(const Vec3f &point) const -{ - int facet_idx = 0; - Vec3d closest_point; - m_emesh.squared_distance(point.cast(), facet_idx, closest_point); - return facet_idx; -} - -} // namespace GUI -} // namespace Slic3r +#include "MeshUtils.hpp" + +#include "libslic3r/Tesselate.hpp" +#include "libslic3r/TriangleMesh.hpp" +#include "libslic3r/TriangleMeshSlicer.hpp" +#include "libslic3r/ClipperUtils.hpp" +#include "libslic3r/Model.hpp" + +#if ENABLE_LEGACY_OPENGL_REMOVAL +#include "slic3r/GUI/GUI_App.hpp" +#endif // ENABLE_LEGACY_OPENGL_REMOVAL +#include "slic3r/GUI/Camera.hpp" +#if ENABLE_GL_SHADERS_ATTRIBUTES +#include "slic3r/GUI/Plater.hpp" +#endif // ENABLE_GL_SHADERS_ATTRIBUTES + +#include + +#include + + +namespace Slic3r { +namespace GUI { + +void MeshClipper::set_plane(const ClippingPlane& plane) +{ + if (m_plane != plane) { + m_plane = plane; + m_triangles_valid = false; + } +} + + +void MeshClipper::set_limiting_plane(const ClippingPlane& plane) +{ + if (m_limiting_plane != plane) { + m_limiting_plane = plane; + m_triangles_valid = false; + } +} + + + +void MeshClipper::set_mesh(const TriangleMesh& mesh) +{ + if (m_mesh != &mesh) { + m_mesh = &mesh; + m_triangles_valid = false; + m_triangles2d.resize(0); + } +} + +void MeshClipper::set_negative_mesh(const TriangleMesh& mesh) +{ + if (m_negative_mesh != &mesh) { + m_negative_mesh = &mesh; + m_triangles_valid = false; + m_triangles2d.resize(0); + } +} + + + +void MeshClipper::set_transformation(const Geometry::Transformation& trafo) +{ + if (! m_trafo.get_matrix().isApprox(trafo.get_matrix())) { + m_trafo = trafo; + m_triangles_valid = false; + m_triangles2d.resize(0); + } +} + +#if ENABLE_LEGACY_OPENGL_REMOVAL +void MeshClipper::render_cut(const ColorRGBA& color) +#else +void MeshClipper::render_cut() +#endif // ENABLE_LEGACY_OPENGL_REMOVAL +{ + if (! m_triangles_valid) + recalculate_triangles(); + +#if ENABLE_LEGACY_OPENGL_REMOVAL + if (m_model.vertices_count() == 0 || m_model.indices_count() == 0) + return; + + GLShaderProgram* curr_shader = wxGetApp().get_current_shader(); + if (curr_shader != nullptr) + curr_shader->stop_using(); + + GLShaderProgram* shader = wxGetApp().get_shader("flat"); + if (shader != nullptr) { + shader->start_using(); +#if ENABLE_GL_SHADERS_ATTRIBUTES + const Camera& camera = wxGetApp().plater()->get_camera(); + shader->set_uniform("view_model_matrix", camera.get_view_matrix()); + shader->set_uniform("projection_matrix", camera.get_projection_matrix()); +#endif // ENABLE_GL_SHADERS_ATTRIBUTES + m_model.set_color(color); + m_model.render(); + shader->stop_using(); + } + + if (curr_shader != nullptr) + curr_shader->start_using(); +#else + if (m_vertex_array.has_VBOs()) + m_vertex_array.render(); +#endif // ENABLE_LEGACY_OPENGL_REMOVAL +} + + + +void MeshClipper::recalculate_triangles() +{ +#if ENABLE_WORLD_COORDINATE + const Transform3f instance_matrix_no_translation_no_scaling = m_trafo.get_rotation_matrix().cast(); +#else + const Transform3f& instance_matrix_no_translation_no_scaling = m_trafo.get_matrix(true,false,true).cast(); +#endif // ENABLE_WORLD_COORDINATE + // Calculate clipping plane normal in mesh coordinates. + const Vec3f up_noscale = instance_matrix_no_translation_no_scaling.inverse() * m_plane.get_normal().cast(); + const Vec3d up = up_noscale.cast().cwiseProduct(m_trafo.get_scaling_factor()); + // Calculate distance from mesh origin to the clipping plane (in mesh coordinates). + const float height_mesh = m_plane.distance(m_trafo.get_offset()) * (up_noscale.norm()/up.norm()); + + // Now do the cutting + MeshSlicingParams slicing_params; + slicing_params.trafo.rotate(Eigen::Quaternion::FromTwoVectors(up, Vec3d::UnitZ())); + + ExPolygons expolys = union_ex(slice_mesh(m_mesh->its, height_mesh, slicing_params)); + + if (m_negative_mesh && !m_negative_mesh->empty()) { + const ExPolygons neg_expolys = union_ex(slice_mesh(m_negative_mesh->its, height_mesh, slicing_params)); + expolys = diff_ex(expolys, neg_expolys); + } + + // Triangulate and rotate the cut into world coords: + Eigen::Quaterniond q; + q.setFromTwoVectors(Vec3d::UnitZ(), up); + Transform3d tr = Transform3d::Identity(); + tr.rotate(q); + tr = m_trafo.get_matrix() * tr; + + if (m_limiting_plane != ClippingPlane::ClipsNothing()) + { + // Now remove whatever ended up below the limiting plane (e.g. sinking objects). + // First transform the limiting plane from world to mesh coords. + // Note that inverse of tr transforms the plane from world to horizontal. + const Vec3d normal_old = m_limiting_plane.get_normal().normalized(); + const Vec3d normal_new = (tr.matrix().block<3,3>(0,0).transpose() * normal_old).normalized(); + + // normal_new should now be the plane normal in mesh coords. To find the offset, + // transform a point and set offset so it belongs to the transformed plane. + Vec3d pt = Vec3d::Zero(); + const double plane_offset = m_limiting_plane.get_data()[3]; + if (std::abs(normal_old.z()) > 0.5) // normal is normalized, at least one of the coords if larger than sqrt(3)/3 = 0.57 + pt.z() = - plane_offset / normal_old.z(); + else if (std::abs(normal_old.y()) > 0.5) + pt.y() = - plane_offset / normal_old.y(); + else + pt.x() = - plane_offset / normal_old.x(); + pt = tr.inverse() * pt; + const double offset = -(normal_new.dot(pt)); + + if (std::abs(normal_old.dot(m_plane.get_normal().normalized())) > 0.99) { + // The cuts are parallel, show all or nothing. + if (normal_old.dot(m_plane.get_normal().normalized()) < 0.0 && offset < height_mesh) + expolys.clear(); + } else { + // The cut is a horizontal plane defined by z=height_mesh. + // ax+by+e=0 is the line of intersection with the limiting plane. + // Normalized so a^2 + b^2 = 1. + const double len = std::hypot(normal_new.x(), normal_new.y()); + if (len == 0.) + return; + const double a = normal_new.x() / len; + const double b = normal_new.y() / len; + const double e = (normal_new.z() * height_mesh + offset) / len; + + // We need a half-plane to limit the cut. Get angle of the intersecting line. + double angle = (b != 0.0) ? std::atan(-a / b) : ((a < 0.0) ? -0.5 * M_PI : 0.5 * M_PI); + if (b > 0) // select correct half-plane + angle += M_PI; + + // We'll take a big rectangle above x-axis and rotate and translate + // it so it lies on our line. This will be the figure to subtract + // from the cut. The coordinates must not overflow after the transform, + // make the rectangle a bit smaller. + const coord_t size = (std::numeric_limits::max() - scale_(std::max(std::abs(e*a), std::abs(e*b)))) / 4; + Polygons ep {Polygon({Point(-size, 0), Point(size, 0), Point(size, 2*size), Point(-size, 2*size)})}; + ep.front().rotate(angle); + ep.front().translate(scale_(-e * a), scale_(-e * b)); + expolys = diff_ex(expolys, ep); + } + } + + m_triangles2d = triangulate_expolygons_2f(expolys, m_trafo.get_matrix().matrix().determinant() < 0.); + + tr.pretranslate(0.001 * m_plane.get_normal().normalized()); // to avoid z-fighting + +#if ENABLE_LEGACY_OPENGL_REMOVAL + m_model.reset(); + + GLModel::Geometry init_data; + init_data.format = { GLModel::Geometry::EPrimitiveType::Triangles, GLModel::Geometry::EVertexLayout::P3N3 }; + init_data.reserve_vertices(m_triangles2d.size()); + init_data.reserve_indices(m_triangles2d.size()); + + // vertices + indices + for (auto it = m_triangles2d.cbegin(); it != m_triangles2d.cend(); it = it + 3) { + init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 0)).x(), (*(it + 0)).y(), height_mesh)).cast(), (Vec3f)up.cast()); + init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 1)).x(), (*(it + 1)).y(), height_mesh)).cast(), (Vec3f)up.cast()); + init_data.add_vertex((Vec3f)(tr * Vec3d((*(it + 2)).x(), (*(it + 2)).y(), height_mesh)).cast(), (Vec3f)up.cast()); + const size_t idx = it - m_triangles2d.cbegin(); + init_data.add_triangle((unsigned int)idx, (unsigned int)idx + 1, (unsigned int)idx + 2); + } + + if (!init_data.is_empty()) + m_model.init_from(std::move(init_data)); +#else + m_vertex_array.release_geometry(); + for (auto it=m_triangles2d.cbegin(); it != m_triangles2d.cend(); it=it+3) { + m_vertex_array.push_geometry(tr * Vec3d((*(it+0))(0), (*(it+0))(1), height_mesh), up); + m_vertex_array.push_geometry(tr * Vec3d((*(it+1))(0), (*(it+1))(1), height_mesh), up); + m_vertex_array.push_geometry(tr * Vec3d((*(it+2))(0), (*(it+2))(1), height_mesh), up); + const size_t idx = it - m_triangles2d.cbegin(); + m_vertex_array.push_triangle(idx, idx+1, idx+2); + } + m_vertex_array.finalize_geometry(true); +#endif // ENABLE_LEGACY_OPENGL_REMOVAL + + m_triangles_valid = true; +} + + +Vec3f MeshRaycaster::get_triangle_normal(size_t facet_idx) const +{ + return m_normals[facet_idx]; +} + +void MeshRaycaster::line_from_mouse_pos(const Vec2d& mouse_pos, const Transform3d& trafo, const Camera& camera, + Vec3d& point, Vec3d& direction) const +{ + Matrix4d modelview = camera.get_view_matrix().matrix(); + Matrix4d projection= camera.get_projection_matrix().matrix(); + Vec4i viewport(camera.get_viewport().data()); + + Vec3d pt1; + Vec3d pt2; + igl::unproject(Vec3d(mouse_pos(0), viewport[3] - mouse_pos(1), 0.), + modelview, projection, viewport, pt1); + igl::unproject(Vec3d(mouse_pos(0), viewport[3] - mouse_pos(1), 1.), + modelview, projection, viewport, pt2); + + Transform3d inv = trafo.inverse(); + pt1 = inv * pt1; + pt2 = inv * pt2; + + point = pt1; + direction = pt2-pt1; +} + + +bool MeshRaycaster::unproject_on_mesh(const Vec2d& mouse_pos, const Transform3d& trafo, const Camera& camera, + Vec3f& position, Vec3f& normal, const ClippingPlane* clipping_plane, + size_t* facet_idx) const +{ + Vec3d point; + Vec3d direction; + line_from_mouse_pos(mouse_pos, trafo, camera, point, direction); + + std::vector hits = m_emesh.query_ray_hits(point, direction); + + if (hits.empty()) + return false; // no intersection found + + unsigned i = 0; + + // Remove points that are obscured or cut by the clipping plane. + // Also, remove anything below the bed (sinking objects). + for (i=0; i= SINKING_Z_THRESHOLD && + (! clipping_plane || ! clipping_plane->is_point_clipped(transformed_hit))) + break; + } + + if (i==hits.size() || (hits.size()-i) % 2 != 0) { + // All hits are either clipped, or there is an odd number of unclipped + // hits - meaning the nearest must be from inside the mesh. + return false; + } + + // Now stuff the points in the provided vector and calculate normals if asked about them: + position = hits[i].position().cast(); + normal = hits[i].normal().cast(); + + if (facet_idx) + *facet_idx = hits[i].face(); + + return true; +} + + +std::vector MeshRaycaster::get_unobscured_idxs(const Geometry::Transformation& trafo, const Camera& camera, const std::vector& points, + const ClippingPlane* clipping_plane) const +{ + std::vector out; + +#if ENABLE_WORLD_COORDINATE + const Transform3d instance_matrix_no_translation_no_scaling = trafo.get_rotation_matrix(); +#else + const Transform3d& instance_matrix_no_translation_no_scaling = trafo.get_matrix(true,false,true); +#endif // ENABLE_WORLD_COORDINATE + Vec3d direction_to_camera = -camera.get_dir_forward(); + Vec3d direction_to_camera_mesh = (instance_matrix_no_translation_no_scaling.inverse() * direction_to_camera).normalized().eval(); + direction_to_camera_mesh = direction_to_camera_mesh.cwiseProduct(trafo.get_scaling_factor()); + const Transform3d inverse_trafo = trafo.get_matrix().inverse(); + + for (size_t i=0; iis_point_clipped(pt.cast())) + continue; + + bool is_obscured = false; + // Cast a ray in the direction of the camera and look for intersection with the mesh: + std::vector hits; + // Offset the start of the ray by EPSILON to account for numerical inaccuracies. + hits = m_emesh.query_ray_hits((inverse_trafo * pt.cast() + direction_to_camera_mesh * EPSILON), + direction_to_camera_mesh); + + if (! hits.empty()) { + // If the closest hit facet normal points in the same direction as the ray, + // we are looking through the mesh and should therefore discard the point: + if (hits.front().normal().dot(direction_to_camera_mesh.cast()) > 0) + is_obscured = true; + + // Eradicate all hits that the caller wants to ignore + for (unsigned j=0; jis_point_clipped(trafo.get_matrix() * hits[j].position())) { + hits.erase(hits.begin()+j); + --j; + } + } + + // FIXME: the intersection could in theory be behind the camera, but as of now we only have camera direction. + // Also, the threshold is in mesh coordinates, not in actual dimensions. + if (! hits.empty()) + is_obscured = true; + } + if (! is_obscured) + out.push_back(i); + } + return out; +} + + +Vec3f MeshRaycaster::get_closest_point(const Vec3f& point, Vec3f* normal) const +{ + int idx = 0; + Vec3d closest_point; + Vec3d pointd = point.cast(); + m_emesh.squared_distance(pointd, idx, closest_point); + if (normal) + // TODO: consider: get_normal(m_emesh, pointd).cast(); + *normal = m_normals[idx]; + + return closest_point.cast(); +} + +int MeshRaycaster::get_closest_facet(const Vec3f &point) const +{ + int facet_idx = 0; + Vec3d closest_point; + m_emesh.squared_distance(point.cast(), facet_idx, closest_point); + return facet_idx; +} + +} // namespace GUI +} // namespace Slic3r diff --git a/src/slic3r/GUI/MeshUtils.hpp b/src/slic3r/GUI/MeshUtils.hpp index 3b8d247c1..df5c3ce43 100644 --- a/src/slic3r/GUI/MeshUtils.hpp +++ b/src/slic3r/GUI/MeshUtils.hpp @@ -4,7 +4,7 @@ #include "libslic3r/Point.hpp" #include "libslic3r/Geometry.hpp" #include "libslic3r/TriangleMesh.hpp" -#include "libslic3r/SLA/IndexedMesh.hpp" +#include "libslic3r/AABBMesh.hpp" #include "admesh/stl.h" #if ENABLE_LEGACY_OPENGL_REMOVAL @@ -167,7 +167,7 @@ public: Vec3f get_triangle_normal(size_t facet_idx) const; private: - sla::IndexedMesh m_emesh; + AABBMesh m_emesh; std::vector m_normals; }; diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp index f7b0df339..a5ed0d029 100644 --- a/tests/sla_print/sla_print_tests.cpp +++ b/tests/sla_print/sla_print_tests.cpp @@ -8,7 +8,6 @@ #include #include -#include namespace { @@ -43,7 +42,7 @@ TEST_CASE("Support point generator should be deterministic if seeded", "[SLASupportGeneration], [SLAPointGen]") { TriangleMesh mesh = load_model("A_upsidedown.obj"); - sla::IndexedMesh emesh{mesh}; + AABBMesh emesh{mesh}; sla::SupportTreeConfig supportcfg; sla::SupportPointGenerator::Config autogencfg; @@ -126,33 +125,69 @@ TEST_CASE("WingedPadAroundObjectIsValid", "[SLASupportGeneration]") { for (auto &fname : AROUND_PAD_TEST_OBJECTS) test_pad(fname, padcfg); } -TEST_CASE("ElevatedSupportGeometryIsValid", "[SLASupportGeneration]") { +TEST_CASE("DefaultSupports::ElevatedSupportGeometryIsValid", "[SLASupportGeneration]") { sla::SupportTreeConfig supportcfg; supportcfg.object_elevation_mm = 10.; for (auto fname : SUPPORT_TEST_MODELS) test_supports(fname, supportcfg); } -TEST_CASE("FloorSupportGeometryIsValid", "[SLASupportGeneration]") { +TEST_CASE("DefaultSupports::FloorSupportGeometryIsValid", "[SLASupportGeneration]") { sla::SupportTreeConfig supportcfg; supportcfg.object_elevation_mm = 0; for (auto &fname: SUPPORT_TEST_MODELS) test_supports(fname, supportcfg); } -TEST_CASE("ElevatedSupportsDoNotPierceModel", "[SLASupportGeneration]") { +TEST_CASE("DefaultSupports::ElevatedSupportsDoNotPierceModel", "[SLASupportGeneration]") { + sla::SupportTreeConfig supportcfg; + supportcfg.object_elevation_mm = 10.; + + for (auto fname : SUPPORT_TEST_MODELS) + test_support_model_collision(fname, supportcfg); +} + +TEST_CASE("DefaultSupports::FloorSupportsDoNotPierceModel", "[SLASupportGeneration]") { sla::SupportTreeConfig supportcfg; + supportcfg.object_elevation_mm = 0; for (auto fname : SUPPORT_TEST_MODELS) test_support_model_collision(fname, supportcfg); } -TEST_CASE("FloorSupportsDoNotPierceModel", "[SLASupportGeneration]") { - +//TEST_CASE("CleverSupports::ElevatedSupportGeometryIsValid", "[SLASupportGeneration][Clever]") { +// sla::SupportTreeConfig supportcfg; +// supportcfg.object_elevation_mm = 10.; +// supportcfg.tree_type = sla::SupportTreeType::Clever; + +// for (auto fname : SUPPORT_TEST_MODELS) test_supports(fname, supportcfg); +//} + +//TEST_CASE("CleverSupports::FloorSupportGeometryIsValid", "[SLASupportGeneration][Clever]") { +// sla::SupportTreeConfig supportcfg; +// supportcfg.object_elevation_mm = 0; +// supportcfg.tree_type = sla::SupportTreeType::Clever; + +// for (auto &fname: SUPPORT_TEST_MODELS) test_supports(fname, supportcfg); +//} + +TEST_CASE("CleverSupports::ElevatedSupportsDoNotPierceModel", "[SLASupportGeneration][Clever]") { + + sla::SupportTreeConfig supportcfg; + supportcfg.object_elevation_mm = 10.; + supportcfg.tree_type = sla::SupportTreeType::Clever; + + for (auto fname : SUPPORT_TEST_MODELS) + test_support_model_collision(fname, supportcfg); +} + +TEST_CASE("CleverSupports::FloorSupportsDoNotPierceModel", "[SLASupportGeneration][Clever]") { + sla::SupportTreeConfig supportcfg; supportcfg.object_elevation_mm = 0; - + supportcfg.tree_type = sla::SupportTreeType::Clever; + for (auto fname : SUPPORT_TEST_MODELS) test_support_model_collision(fname, supportcfg); } diff --git a/tests/sla_print/sla_raycast_tests.cpp b/tests/sla_print/sla_raycast_tests.cpp index a69e09412..1f2962bec 100644 --- a/tests/sla_print/sla_raycast_tests.cpp +++ b/tests/sla_print/sla_raycast_tests.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include #include "sla_test_utils.hpp" diff --git a/tests/sla_print/sla_supptgen_tests.cpp b/tests/sla_print/sla_supptgen_tests.cpp index f2cf7e834..e5e602438 100644 --- a/tests/sla_print/sla_supptgen_tests.cpp +++ b/tests/sla_print/sla_supptgen_tests.cpp @@ -3,6 +3,7 @@ #include #include +#include #include "sla_test_utils.hpp" diff --git a/tests/sla_print/sla_test_utils.cpp b/tests/sla_print/sla_test_utils.cpp index a58b33ca0..07f19300c 100644 --- a/tests/sla_print/sla_test_utils.cpp +++ b/tests/sla_print/sla_test_utils.cpp @@ -1,6 +1,7 @@ #include "sla_test_utils.hpp" #include "libslic3r/TriangleMeshSlicer.hpp" #include "libslic3r/SLA/AGGRaster.hpp" +#include "libslic3r/SLA/DefaultSupportTree.hpp" #include @@ -15,14 +16,16 @@ void test_support_model_collision(const std::string &obj_filename, // Set head penetration to a small negative value which should ensure that // the supports will not touch the model body. - supportcfg.head_penetration_mm = -0.15; + supportcfg.head_penetration_mm = -0.2; test_supports(obj_filename, supportcfg, hollowingcfg, drainholes, byproducts); // Slice the support mesh given the slice grid of the model. std::vector support_slices = - byproducts.supporttree.slice(byproducts.slicegrid, CLOSING_RADIUS); - + sla::slice(byproducts.supporttree.retrieve_mesh(sla::MeshType::Support), + byproducts.supporttree.retrieve_mesh(sla::MeshType::Pad), + byproducts.slicegrid, CLOSING_RADIUS, {}); + // The slices originate from the same slice grid so the numbers must match bool support_mesh_is_empty = @@ -47,13 +50,15 @@ void test_support_model_collision(const std::string &obj_filename, notouch = notouch && area(intersections) < PI * pinhead_r * pinhead_r; } - /*if (!notouch) */export_failed_case(support_slices, byproducts); + if (!notouch) + export_failed_case(support_slices, byproducts); REQUIRE(notouch); } void export_failed_case(const std::vector &support_slices, const SupportByproducts &byproducts) { + bool do_export_stl = false; for (size_t n = 0; n < support_slices.size(); ++n) { const ExPolygons &sup_slice = support_slices[n]; const ExPolygons &mod_slice = byproducts.model_slices[n]; @@ -68,14 +73,18 @@ void export_failed_case(const std::vector &support_slices, const Sup svg.draw(intersections, "red"); svg.Close(); } + + do_export_stl = do_export_stl || !intersections.empty(); } - indexed_triangle_set its; - byproducts.supporttree.retrieve_full_mesh(its); - TriangleMesh m{its}; - m.merge(byproducts.input_mesh); - m.WriteOBJFile((Catch::getResultCapture().getCurrentTestName() + "_" + - byproducts.obj_fname).c_str()); + if (do_export_stl) { + indexed_triangle_set its; + byproducts.supporttree.retrieve_full_mesh(its); + TriangleMesh m{its}; + m.merge(byproducts.input_mesh); + m.WriteOBJFile((Catch::getResultCapture().getCurrentTestName() + "_" + + byproducts.obj_fname).c_str()); + } } void test_supports(const std::string &obj_filename, @@ -107,7 +116,7 @@ void test_supports(const std::string &obj_filename, // Create the special index-triangle mesh with spatial indexing which // is the input of the support point and support mesh generators - sla::IndexedMesh emesh{mesh}; + AABBMesh emesh{mesh}; #ifdef SLIC3R_HOLE_RAYCASTER if (hollowingcfg.enabled) @@ -144,10 +153,16 @@ void test_supports(const std::string &obj_filename, // Generate the actual support tree sla::SupportTreeBuilder treebuilder; sla::SupportableMesh sm{emesh, support_points, supportcfg}; - sla::SupportTreeBuildsteps::execute(treebuilder, sm); - - check_support_tree_integrity(treebuilder, supportcfg); - + + switch (sm.cfg.tree_type) { + case sla::SupportTreeType::Default: { + sla::DefaultSupportTree::execute(treebuilder, sm); + check_support_tree_integrity(treebuilder, supportcfg, sla::ground_level(sm)); + break; + } + default:; + } + TriangleMesh output_mesh{treebuilder.retrieve_mesh(sla::MeshType::Support)}; check_validity(output_mesh, validityflags); @@ -171,9 +186,9 @@ void test_supports(const std::string &obj_filename, } void check_support_tree_integrity(const sla::SupportTreeBuilder &stree, - const sla::SupportTreeConfig &cfg) + const sla::SupportTreeConfig &cfg, + double gnd) { - double gnd = stree.ground_level; double H1 = cfg.max_solo_pillar_height_mm; double H2 = cfg.max_dual_pillar_height_mm; @@ -426,7 +441,7 @@ sla::SupportPoints calc_support_pts( std::vector slices = slice_mesh_ex(mesh.its, heights, CLOSING_RADIUS); // Prepare the support point calculator - sla::IndexedMesh emesh{mesh}; + AABBMesh emesh{mesh}; sla::SupportPointGenerator spgen{emesh, cfg, []{}, [](int){}}; // Calculate the support points diff --git a/tests/sla_print/sla_test_utils.hpp b/tests/sla_print/sla_test_utils.hpp index 7171914d4..c5943bd5c 100644 --- a/tests/sla_print/sla_test_utils.hpp +++ b/tests/sla_print/sla_test_utils.hpp @@ -14,7 +14,7 @@ #include "libslic3r/TriangleMesh.hpp" #include "libslic3r/SLA/Pad.hpp" #include "libslic3r/SLA/SupportTreeBuilder.hpp" -#include "libslic3r/SLA/SupportTreeBuildsteps.hpp" +#include "libslic3r/SLA/SupportTreeUtils.hpp" #include "libslic3r/SLA/SupportPointGenerator.hpp" #include "libslic3r/SLA/AGGRaster.hpp" #include "libslic3r/SLA/ConcaveHull.hpp" @@ -67,7 +67,8 @@ struct SupportByproducts const constexpr float CLOSING_RADIUS = 0.005f; void check_support_tree_integrity(const sla::SupportTreeBuilder &stree, - const sla::SupportTreeConfig &cfg); + const sla::SupportTreeConfig &cfg, + double gnd); void test_supports(const std::string &obj_filename, const sla::SupportTreeConfig &supportcfg,