From 73af7c64b8028a11d0123a95e236af014dc28008 Mon Sep 17 00:00:00 2001
From: Lukas Matena <lukasmatena@seznam.cz>
Date: Tue, 26 Nov 2019 11:36:09 +0100
Subject: [PATCH] SLATreeSupports generator now takes account for holes and can
 build supports through them

---
 src/libslic3r/SLA/Common.cpp                | 113 +++++++++++---------
 src/libslic3r/SLA/EigenMesh3D.hpp           |  44 ++++----
 src/libslic3r/SLA/Hollowing.cpp             |  10 +-
 src/libslic3r/SLA/SupportPointGenerator.cpp |   7 --
 src/libslic3r/SLA/SupportTreeBuildsteps.cpp |   2 +-
 src/libslic3r/SLAPrintSteps.cpp             |   6 ++
 6 files changed, 103 insertions(+), 79 deletions(-)

diff --git a/src/libslic3r/SLA/Common.cpp b/src/libslic3r/SLA/Common.cpp
index c54d78afa..3710bf3be 100644
--- a/src/libslic3r/SLA/Common.cpp
+++ b/src/libslic3r/SLA/Common.cpp
@@ -263,16 +263,13 @@ EigenMesh3D &EigenMesh3D::operator=(const EigenMesh3D &other)
 }
 
 EigenMesh3D::hit_result
-EigenMesh3D::query_ray_hit(const Vec3d &s,
-                           const Vec3d &dir,
-                           const std::vector<DrainHole>* holes
-                          ) const
+EigenMesh3D::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
 {
     assert(is_approx(dir.norm(), 1.));
     igl::Hit hit;
     hit.t = std::numeric_limits<float>::infinity();
 
-    if (! holes) {
+    if (m_holes.empty()) {
         m_aabb->intersect_ray(m_V, m_F, s, dir, hit);
         hit_result ret(*this);
         ret.m_t = double(hit.t);
@@ -286,7 +283,7 @@ EigenMesh3D::query_ray_hit(const Vec3d &s,
     else {
         // If there are holes, the hit_results will be made by
         // query_ray_hits (object) and filter_hits (holes):
-        return filter_hits(query_ray_hits(s, dir), *holes);
+        return filter_hits(query_ray_hits(s, dir));
     }
 }
 
@@ -316,46 +313,59 @@ EigenMesh3D::query_ray_hits(const Vec3d &s, const Vec3d &dir) const
 }
 
 EigenMesh3D::hit_result EigenMesh3D::filter_hits(
-                     const std::vector<EigenMesh3D::hit_result>& object_hits,
-                     const std::vector<DrainHole>& holes) const
+                     const std::vector<EigenMesh3D::hit_result>& object_hits) const
 {
+    assert(! m_holes.empty());
     hit_result out(*this);
-    out.m_t = std::nan("");
 
-    if (! holes.empty() && ! object_hits.empty()) {
-        Vec3d s = object_hits.front().source();
-        Vec3d dir = object_hits.front().direction();
+    if (object_hits.empty())
+        return out;
 
-        struct HoleHit {
-            HoleHit(float t_p, const Vec3d& normal_p, bool entry_p) :
-                t(t_p), normal(normal_p), entry(entry_p) {}
-            float t;
-            Vec3d normal;
-            bool entry;
-        };
-        std::vector<HoleHit> hole_isects;
+    const Vec3d& s = object_hits.front().source();
+    const Vec3d& dir = object_hits.front().direction();
 
-        // Collect hits on all holes, preserve information about entry/exit
-        for (const sla::DrainHole& hole : holes) {
-            std::array<std::pair<float, Vec3d>, 2> isects;
-            if (hole.get_intersections(s.cast<float>(),
-                                       dir.cast<float>(), isects)) {
-                hole_isects.emplace_back(isects[0].first, isects[0].second, true);
-                hole_isects.emplace_back(isects[1].first, isects[1].second, false);
-            }
+    // A helper struct to save an intersetion with a hole
+    struct HoleHit {
+        HoleHit(float t_p, const Vec3d& normal_p, bool entry_p) :
+            t(t_p), normal(normal_p), entry(entry_p) {}
+        float t;
+        Vec3d normal;
+        bool entry;
+    };
+    std::vector<HoleHit> hole_isects;
+
+    // Collect hits on all holes, preserve information about entry/exit
+    for (const sla::DrainHole& hole : m_holes) {
+        std::array<std::pair<float, Vec3d>, 2> isects;
+        if (hole.get_intersections(s.cast<float>(),
+                                   dir.cast<float>(), isects)) {
+            hole_isects.emplace_back(isects[0].first, isects[0].second, true);
+            hole_isects.emplace_back(isects[1].first, isects[1].second, false);
         }
-        // Holes can intersect each other, sort the hits by t
-        std::sort(hole_isects.begin(), hole_isects.end(),
-                  [](const HoleHit& a, const HoleHit& b) { return a.t < b.t; });
+    }
+    // Remove hole hits behind the source
+    for (int i=0; i<int(hole_isects.size()); ++i)
+        if (hole_isects[i].t < 0.f)
+            hole_isects.erase(hole_isects.begin() + (i--));
 
-        // Now inspect the intersections with object and holes, keep track how
-        // deep are we nested in mesh/holes and pick the correct intersection
-        int hole_nested = 0;
-        int object_nested = 0;
+    // Holes can intersect each other, sort the hits by t
+    std::sort(hole_isects.begin(), hole_isects.end(),
+              [](const HoleHit& a, const HoleHit& b) { return a.t < b.t; });
+
+    // Now inspect the intersections with object and holes, in the order of
+    // increasing distance. Keep track how deep are we nested in mesh/holes and
+    // pick the correct intersection.
+    // This needs to be done twice - first to find out how deep in the structure
+    // the source is, then to pick the correct intersection.
+    int hole_nested = 0;
+    int object_nested = 0;
+    for (int dry_run=1; dry_run>=0; --dry_run) {
+        hole_nested = -hole_nested;
+        object_nested = -object_nested;
 
         bool is_hole = false;
         bool is_entry = false;
-        const HoleHit* next_hole_hit = &hole_isects.front();
+        const HoleHit* next_hole_hit = hole_isects.empty() ? nullptr : &hole_isects.front();
         const hit_result* next_mesh_hit = &object_hits.front();
 
         while (next_hole_hit || next_mesh_hit) {
@@ -367,23 +377,25 @@ EigenMesh3D::hit_result EigenMesh3D::filter_hits(
             // Is this entry or exit hit?
             is_entry = is_hole ? next_hole_hit->entry : ! next_mesh_hit->is_inside();
 
-            if (! is_hole && is_entry && hole_nested == 0) {
-                // This mesh point is the one we seek
-                return *next_mesh_hit;
-            }
-            if (is_hole && ! is_entry && object_nested != 0) {
-                // This holehit is the one we seek
-                out.m_t = next_hole_hit->t;
-                out.m_normal = next_hole_hit->normal;
-                out.m_source = s;
-                out.m_dir = dir;
-                return out;
+            if (! dry_run) {
+                if (! is_hole && hole_nested == 0) {
+                    // This is a valid object hit
+                    return *next_mesh_hit;
+                }
+                if (is_hole && ! is_entry && object_nested != 0) {
+                    // This holehit is the one we seek
+                    out.m_t = next_hole_hit->t;
+                    out.m_normal = next_hole_hit->normal;
+                    out.m_source = s;
+                    out.m_dir = dir;
+                    return out;
+                }
             }
 
-            hole_nested += (is_hole ? (is_entry ? 1 : -1) : 0);
-            object_nested += (! is_hole ? (is_entry ? 1 : -1) : 0);
+            // Increase/decrease the counter
+            (is_hole ? hole_nested : object_nested) += (is_entry ? 1 : -1);
 
-            // Advance the pointer
+            // Advance the respective pointer
             if (is_hole && next_hole_hit++ == &hole_isects.back())
                 next_hole_hit = nullptr;
             if (! is_hole && next_mesh_hit++ == &object_hits.back())
@@ -391,6 +403,7 @@ EigenMesh3D::hit_result EigenMesh3D::filter_hits(
         }
     }
 
+    // if we got here, the ray ended up in infinity
     return out;
 }
 
diff --git a/src/libslic3r/SLA/EigenMesh3D.hpp b/src/libslic3r/SLA/EigenMesh3D.hpp
index c92094756..15176f0cb 100644
--- a/src/libslic3r/SLA/EigenMesh3D.hpp
+++ b/src/libslic3r/SLA/EigenMesh3D.hpp
@@ -2,6 +2,7 @@
 #define SLA_EIGENMESH3D_H
 
 #include <libslic3r/SLA/Common.hpp>
+#include "libslic3r/SLA/Hollowing.hpp"
 
 namespace Slic3r {
 
@@ -10,7 +11,6 @@ class TriangleMesh;
 namespace sla {
 
 struct Contour3D;
-struct DrainHole;
 
 /// An index-triangle structure for libIGL functions. Also serves as an
 /// alternative (raw) input format for the SLASupportTree.
@@ -23,6 +23,11 @@ class EigenMesh3D {
     double m_ground_level = 0, m_gnd_offset = 0;
     
     std::unique_ptr<AABBImpl> m_aabb;
+
+    // This holds a copy of holes in the mesh. Initialized externally
+    // by load_mesh setter.
+    std::vector<DrainHole> m_holes;
+
 public:
     
     EigenMesh3D(const TriangleMesh&);
@@ -41,57 +46,58 @@ public:
     
     // Result of a raycast
     class hit_result {
-        double m_t = std::nan("");
+        // m_t holds a distance from m_source to the intersection.
+        double m_t = infty();
         const EigenMesh3D *m_mesh = nullptr;
         Vec3d m_dir;
         Vec3d m_source;
-        Vec3d m_normal = Vec3d::Zero();
+        Vec3d m_normal;
         friend class EigenMesh3D;
         
         // A valid object of this class can only be obtained from
         // EigenMesh3D::query_ray_hit method.
         explicit inline hit_result(const EigenMesh3D& em): m_mesh(&em) {}
     public:
+        // This denotes no hit on the mesh.
+        static inline constexpr double infty() { return std::numeric_limits<double>::infinity(); }
         
-        // This can create a placeholder object which is invalid (not created
-        // by a query_ray_hit call) but the distance can be preset to
-        // a specific value for distinguishing the placeholder.
-        inline hit_result(double val = std::nan("")): m_t(val) {}
+        explicit inline hit_result(double val = infty()) : m_t(val) {}
         
         inline double distance() const { return m_t; }
         inline const Vec3d& direction() const { return m_dir; }
         inline const Vec3d& source() const { return m_source; }
         inline Vec3d position() const { return m_source + m_dir * m_t; }
         inline bool is_valid() const { return m_mesh != nullptr; }
-        inline bool is_hit() const { return m_normal != Vec3d::Zero(); }
-        
+        inline bool is_hit() const { return m_t != infty(); }
+
         // Hit_result can decay into a double as the hit distance.
         inline operator double() const { return distance(); }
 
         inline const Vec3d& normal() const {
-            if(!is_valid())
-                throw std::runtime_error("EigenMesh3D::hit_result::normal() "
-                                         "called on invalid object.");
+            assert(is_valid());
             return m_normal;
         }
-        
+
         inline bool is_inside() const {
-            return normal().dot(m_dir) > 0;
+            return is_hit() && normal().dot(m_dir) > 0;
         }
     };
     
+    // Inform the object about location of holes
+    // creates internal copy of the vector
+    void load_holes(const std::vector<DrainHole>& holes) {
+        m_holes = holes;
+    }
+
     // Casting a ray on the mesh, returns the distance where the hit occures.
-    hit_result query_ray_hit(const Vec3d &s,
-                             const Vec3d &dir,
-                             const std::vector<DrainHole>* holes = nullptr) const;
+    hit_result query_ray_hit(const Vec3d &s, const Vec3d &dir) const;
     
     // Casts a ray on the mesh and returns all hits
     std::vector<hit_result> query_ray_hits(const Vec3d &s, const Vec3d &dir) const;
 
     // Iterates over hits and holes and returns the true hit, possibly
     // on the inside of a hole.
-    hit_result filter_hits(const std::vector<EigenMesh3D::hit_result>& obj_hits,
-                           const std::vector<DrainHole>& holes) const;
+    hit_result filter_hits(const std::vector<EigenMesh3D::hit_result>& obj_hits) const;
 
     class si_result {
         double m_value;
diff --git a/src/libslic3r/SLA/Hollowing.cpp b/src/libslic3r/SLA/Hollowing.cpp
index a510ab3da..46e60f91c 100644
--- a/src/libslic3r/SLA/Hollowing.cpp
+++ b/src/libslic3r/SLA/Hollowing.cpp
@@ -4,6 +4,7 @@
 #include <libslic3r/TriangleMesh.hpp>
 #include <libslic3r/SLA/Hollowing.hpp>
 #include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/SLA/EigenMesh3D.hpp>
 
 #include <boost/log/trivial.hpp>
 
@@ -152,7 +153,7 @@ bool DrainHole::get_intersections(const Vec3f& s, const Vec3f& dir,
     const Eigen::ParametrizedLine<float, 3> ray(s, dir.normalized());
 
     for (size_t i=0; i<2; ++i)
-        out[i] = std::make_pair(std::nan(""), Vec3d::Zero());
+        out[i] = std::make_pair(sla::EigenMesh3D::hit_result::infty(), Vec3d::Zero());
 
     const float sqr_radius = pow(radius, 2.f);
 
@@ -170,6 +171,11 @@ bool DrainHole::get_intersections(const Vec3f& s, const Vec3f& dir,
     if (! is_approx(ray.direction().dot(normal), 0.f)) {
         for (size_t i=1; i<=1; --i) {
             Vec3f cylinder_center = pos+i*height*normal;
+            if (i == 0) {
+                // The hole base can be identical to mesh surface if it is flat
+                // let's better move the base outward a bit
+                cylinder_center -= EPSILON*normal;
+            }
             base = Eigen::Hyperplane<float, 3>(normal, cylinder_center);
             Vec3f intersection = ray.intersectionPoint(base);
             // Only accept the point if it is inside the cylinder base.
@@ -184,7 +190,7 @@ bool DrainHole::get_intersections(const Vec3f& s, const Vec3f& dir,
     {
         // In case the line was perpendicular to the cylinder axis, previous
         // block was skipped, but base will later be assumed to be valid.
-        base = Eigen::Hyperplane<float, 3>(normal, pos);
+        base = Eigen::Hyperplane<float, 3>(normal, pos-EPSILON*normal);
     }
 
     // In case there is still an intersection to be found, check the wall
diff --git a/src/libslic3r/SLA/SupportPointGenerator.cpp b/src/libslic3r/SLA/SupportPointGenerator.cpp
index 36361f9ca..66b6ea3f4 100644
--- a/src/libslic3r/SLA/SupportPointGenerator.cpp
+++ b/src/libslic3r/SLA/SupportPointGenerator.cpp
@@ -77,9 +77,6 @@ void SupportPointGenerator::project_onto_mesh(std::vector<sla::SupportPoint>& po
                     m_throw_on_cancel();
                 Vec3f& p = points[point_id].pos;
                 // Project the point upward and downward and choose the closer intersection with the mesh.
-                //bool up   = igl::ray_mesh_intersect(p.cast<float>(), Vec3f(0., 0., 1.), m_V, m_F, hit_up);
-                //bool down = igl::ray_mesh_intersect(p.cast<float>(), Vec3f(0., 0., -1.), m_V, m_F, hit_down);
-
                 sla::EigenMesh3D::hit_result hit_up   = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., 1.));
                 sla::EigenMesh3D::hit_result hit_down = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., -1.));
 
@@ -90,10 +87,6 @@ void SupportPointGenerator::project_onto_mesh(std::vector<sla::SupportPoint>& po
                     continue;
 
                 sla::EigenMesh3D::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down;
-                //int fid = hit.face();
-                //Vec3f bc(1-hit.u-hit.v, hit.u, hit.v);
-                //p = (bc(0) * m_V.row(m_F(fid, 0)) + bc(1) * m_V.row(m_F(fid, 1)) + bc(2)*m_V.row(m_F(fid, 2))).cast<float>();
-
                 p = p + (hit.distance() * hit.direction()).cast<float>();
             }
         });
diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
index 46a6160c1..68afb7391 100644
--- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
+++ b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
@@ -778,7 +778,7 @@ void SupportTreeBuildsteps::filter()
                     nn = Vec3d(std::cos(azimuth) * std::sin(polar),
                                std::sin(azimuth) * std::sin(polar),
                                std::cos(polar)).normalized();
-                    t = oresult.score;
+                    t = EigenMesh3D::hit_result(oresult.score);
                 }
             }
             
diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp
index 55359cc5c..3135c33dd 100644
--- a/src/libslic3r/SLAPrintSteps.cpp
+++ b/src/libslic3r/SLAPrintSteps.cpp
@@ -253,6 +253,11 @@ void SLAPrint::Steps::support_points(SLAPrintObject &po)
         
         // calculate heights of slices (slices are calculated already)
         const std::vector<float>& heights = po.m_model_height_levels;
+
+        // Tell the mesh where drain holes are. Although the points are
+        // calculated on slices, the algorithm then raycasts the points
+        // so they actually lie on the mesh.
+        po.m_supportdata->emesh.load_holes(po.transformed_drainhole_points());
         
         throw_if_canceled();
         sla::SupportPointGenerator::Config config;
@@ -321,6 +326,7 @@ void SLAPrint::Steps::support_tree(SLAPrintObject &po)
         po.m_supportdata->emesh.ground_level_offset(pcfg.wall_thickness_mm);
     
     po.m_supportdata->cfg = make_support_cfg(po.m_config);
+    po.m_supportdata->emesh.load_holes(po.transformed_drainhole_points());
     
     // scaling for the sub operations
     double d = objectstep_scale * OBJ_STEP_LEVELS[slaposSupportTree] / 100.0;