Started work on extending EigenMesh3D to account for possible drain holes when raycasting

src/libslic3r/SLA/Common.cpp

@@ -6,6 +6,7 @@
 #include <libslic3r/SLA/EigenMesh3D.hpp>
 #include <libslic3r/SLA/Contour3D.hpp>
 #include <libslic3r/SLA/Clustering.hpp>
+#include <libslic3r/SLA/Hollowing.hpp>
 
 
 // Workaround: IGL signed_distance.h will define PI in the igl namespace.
@@ -181,6 +182,19 @@ void BoxIndex::foreach(std::function<void (const BoxIndexEl &)> fn)
     for(auto& el : m_impl->m_store) fn(el);
 }
 
+
+namespace {
+// Iterates over hits and holes and returns the true hit, possibly
+// on the inside of a hole. Free function so it can return igl::Hit
+// without including igl in a header.
+igl::Hit filter_hits(const std::vector<EigenMesh3D::hit_result>& hits,
+                     const std::vector<DrainHole>& holes)
+{
+    return igl::Hit();
+}
+
+} // namespace
+
 /* ****************************************************************************
  * EigenMesh3D implementation
  * ****************************************************************************/
@@ -261,11 +275,18 @@ EigenMesh3D &EigenMesh3D::operator=(const EigenMesh3D &other)
 }
 
 EigenMesh3D::hit_result
-EigenMesh3D::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
+EigenMesh3D::query_ray_hit(const Vec3d &s,
+                           const Vec3d &dir,
+                           const std::vector<DrainHole>* holes
+                          ) const
 {
     igl::Hit hit;
     hit.t = std::numeric_limits<float>::infinity();
+
+    if (! holes)
         m_aabb->intersect_ray(m_V, m_F, s, dir, hit);
+    else
+        hit = filter_hits(query_ray_hits(s, dir), *holes);
     
     hit_result ret(*this);
     ret.m_t = double(hit.t);

src/libslic3r/SLA/EigenMesh3D.hpp

@@ -10,10 +10,11 @@ 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.
-//  Implemented in SLASupportTreeIGL.cpp
+//  Implemented in libslic3r/SLA/Common.cpp
 class EigenMesh3D {
     class AABBImpl;
     
@@ -83,7 +84,9 @@ public:
     };
     
     // 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;
+    hit_result query_ray_hit(const Vec3d &s,
+                             const Vec3d &dir,
+                             const std::vector<DrainHole>* holes = nullptr) 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;

src/libslic3r/SLA/Hollowing.cpp

@@ -126,4 +126,100 @@ bool DrainHole::operator==(const DrainHole &sp) const
             is_approx(height, sp.height);
 }
 
+bool DrainHole::is_inside(const Vec3f& pt) const
+{
+    Eigen::Hyperplane<float, 3> plane(normal, pos);
+    float dist = plane.signedDistance(pt);
+    if (dist < EPSILON || dist > height)
+        return false;
+
+    Eigen::ParametrizedLine<float, 3> axis(pos, normal);
+    if ( axis.squaredDistance(pt) < pow(radius, 2.f))
+        return true;
+
+    return false;
+}
+
+
+// Given a line s+dir*t, return parameter t of intersections  with the hole.
+// If there is no intersection, returns nan.
+std::pair<float, float> DrainHole::get_intersections(const Vec3f& s,
+                                                     const Vec3f& dir) const
+{
+    assert(is_approx(normal.norm(), 1.f));
+    const Eigen::ParametrizedLine<float, 3> ray(s, dir.normalized());
+
+    std::pair<float, float> out(std::nan(""), std::nan(""));
+    const float sqr_radius = pow(radius, 2.f);
+
+    // first check a bounding sphere of the hole:
+    Vec3f center = pos+normal*height/2.f;
+    float sqr_dist_limit = pow(height/2.f, 2.f) + sqr_radius ;
+    if (ray.squaredDistance(center) > sqr_dist_limit)
+        return out;
+
+    // The line intersects the bounding sphere, look for intersections with
+    // bases of the cylinder.
+
+    size_t found = 0; // counts how many intersections were found
+    Eigen::Hyperplane<float, 3> base;
+    if (! is_approx(ray.direction().dot(normal), 0.f)) {
+        for (size_t i=1; i<=1; --i) {
+            Vec3f cylinder_center = pos+i*height*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.
+            if ((cylinder_center-intersection).squaredNorm() < sqr_radius) {
+                (found ? out.second : out.first) = ray.intersectionParameter(base);
+                ++found;
+            }
+        }
+    }
+    else
+    {
+        // 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);
+    }
+
+    // In case there is still an intersection to be found, check the wall
+    if (found != 2 && ! is_approx(std::abs(ray.direction().dot(normal)), 1.f)) {
+        // Project the ray onto the base plane
+        Vec3f proj_origin = base.projection(ray.origin());
+        Vec3f proj_dir = base.projection(ray.origin()+ray.direction())-proj_origin;
+        // save how the parameter scales and normalize the projected direction
+        float par_scale = proj_dir.norm();
+        proj_dir = proj_dir/par_scale;
+        Eigen::ParametrizedLine<float, 3> projected_ray(proj_origin, proj_dir);
+        // Calculate point on the secant that's closest to the center
+        // and its distance to the circle along the projected line
+        Vec3f closest = projected_ray.projection(pos);
+        float dist = sqrt((sqr_radius - (closest-pos).squaredNorm()));
+        // Unproject both intersections on the original line and check
+        // they are on the cylinder and not past it:
+        for (int i=-1; i<=1 && found !=2; i+=2) {
+            Vec3f isect = closest + i*dist * projected_ray.direction();
+            float par = (isect-proj_origin).norm() / par_scale;
+            isect = ray.pointAt(par);
+            // check that the intersection is between the base planes:
+            float vert_dist = base.signedDistance(isect);
+            if (vert_dist > 0.f && vert_dist < height) {
+                (found ? out.second : out.first) = par;
+                ++found;
+            }
+        }
+    }
+
+    // If only one intersection was found, it is some corner case,
+    // no intersection will be returned:
+    if (found != 0)
+        return std::pair<float, float>(std::nan(""), std::nan(""));
+
+    // Sort the intersections:
+    if (out.first > out.second)
+        std::swap(out.first, out.second);
+
+    return out;
+}
+
 }} // namespace Slic3r::sla

src/libslic3r/SLA/Hollowing.hpp

@@ -37,6 +37,12 @@ struct DrainHole
     
     bool operator!=(const DrainHole &sp) const { return !(sp == (*this)); }
 
+    bool is_inside(const Vec3f& pt) const;
+
+    std::pair<float, float> get_intersections(const Vec3f& s,
+                                              const Vec3f& dir
+                                             ) const;
+    
     template<class Archive> inline void serialize(Archive &ar)
     {
         ar(pos, normal, radius, height);

src/libslic3r/SLAPrint.cpp

@@ -415,7 +415,8 @@ SLAPrint::ApplyStatus SLAPrint::apply(const Model &model, DynamicPrintConfig con
                 }
                 model_object.sla_points_status = model_object_new.sla_points_status;
                 
-                if (model_object.sla_drain_holes.size() != model_object_new.sla_drain_holes.size())
+                // Invalidate hollowing if drain holes have changed
+                if (model_object.sla_drain_holes != model_object_new.sla_drain_holes)
                 {
                     model_object.sla_drain_holes = model_object_new.sla_drain_holes;
                     update_apply_status(it_print_object_status->print_object->invalidate_step(slaposHollowing));

src/slic3r/GUI/Gizmos/GLGizmoSlaSupports.cpp

@@ -412,25 +412,6 @@ void GLGizmoSlaSupports::update_mesh()
 }
 
 
-bool GLGizmoSlaSupports::is_point_in_hole(const Vec3f& pt) const
-{
-    auto squared_distance_from_line = [](const Vec3f pt, const Vec3f& line_pt, const Vec3f& dir) -> float {
-        Vec3f diff = line_pt - pt;
-        return (diff - diff.dot(dir) * dir).squaredNorm();
-    };
-
-
-    for (const sla::DrainHole& hole : m_model_object->sla_drain_holes) {
-        if ( hole.normal.dot(pt-hole.pos) < EPSILON
-         || hole.normal.dot(pt-(hole.pos+hole.height * hole.normal)) > 0.f)
-            continue;
-        if ( squared_distance_from_line(pt, hole.pos, hole.normal) < pow(hole.radius, 2.f))
-            return true;
-    }
-
-    return false;
-}
-
 // Unprojects the mouse position on the mesh and saves hit point and normal of the facet into pos_and_normal
 // Return false if no intersection was found, true otherwise.
 bool GLGizmoSlaSupports::unproject_on_mesh(const Vec2d& mouse_pos, std::pair<Vec3f, Vec3f>& pos_and_normal)
@@ -448,13 +429,22 @@ bool GLGizmoSlaSupports::unproject_on_mesh(const Vec2d& mouse_pos, std::pair<Vec
     // The raycaster query
     Vec3f hit;
     Vec3f normal;
-    if (m_mesh_raycaster->unproject_on_mesh(mouse_pos, trafo.get_matrix(), camera, hit, normal, m_clipping_plane.get())
+    if (m_mesh_raycaster->unproject_on_mesh(mouse_pos, trafo.get_matrix(), camera, hit, normal, m_clipping_plane.get())) {
-     && ! is_point_in_hole(hit)) {
+        // Check whether the hit is in a hole
+        bool in_hole = false;
+        for (const sla::DrainHole& hole : m_model_object->sla_drain_holes) {
+            if (hole.is_inside(hit)) {
+                in_hole = true;
+                break;
+            }
+        }
+        if (! in_hole) {
             // Return both the point and the facet normal.
             pos_and_normal = std::make_pair(hit, normal);
             return true;
         }
-    else
+    }
+
     return false;
 }