Create smaller supports in problematic areas with established strategies

Completely remove the concept of CompactBridge.

Replace it with Heads having the same back radius as front radius. 

Try to apply the same rules for mini supports as in the route_to_model step.

Increased accuracy of bridge_mesh_intersect shot from support points


Refining mini support integration

1 files changed, 205 insertions, 217 deletions

diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
index 29ad6057f..df9de3555 100644
--- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
+++ b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
@@ -42,6 +42,8 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder &   builder,
 {
     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
@@ -166,64 +168,6 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder &   builder,
     return pc == ABORT;
 }
 
-// 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<size_t N>
-class PointRing {
-    std::array<double, N> 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<N>(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<class C, class Hit = EigenMesh3D::hit_result> 
 static Hit min_hit(const C &hits)
 {
@@ -312,7 +256,7 @@ EigenMesh3D::hit_result SupportTreeBuildsteps::pinhead_mesh_intersect(
 }
 
 EigenMesh3D::hit_result SupportTreeBuildsteps::bridge_mesh_intersect(
-    const Vec3d &src, const Vec3d &dir, double r, bool ins_check)
+    const Vec3d &src, const Vec3d &dir, double r, double safety_d)
 {
     static const size_t SAMPLES = 8;
     PointRing<SAMPLES> ring{dir};
@@ -321,16 +265,19 @@ EigenMesh3D::hit_result SupportTreeBuildsteps::bridge_mesh_intersect(
     
     // Hit results
     std::array<Hit, SAMPLES> hits;
+
+    double sd = std::isnan(safety_d) ? m_cfg.safety_distance_mm : safety_d;
+    sd = sd * r / m_cfg.head_back_radius_mm;
+
+    bool ins_check = sd < m_cfg.safety_distance_mm;
     
     ccr::enumerate(hits.begin(), hits.end(), 
-                   [this, r, src, ins_check, &ring, dir] (Hit &hit, size_t i) {
-        
-        const double sd = m_cfg.safety_distance_mm;
-        
+                [this, r, src, ins_check, &ring, dir, sd] (Hit &hit, size_t 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 + sd * dir, dir);
+        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);
@@ -460,7 +407,7 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head,
     
     Vec3d bridgestart = headjp;
     Vec3d bridgeend = nearjp_u;
-    double max_len = m_cfg.max_bridge_length_mm;
+    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;
     
@@ -494,7 +441,7 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head,
     
     // 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 + 2 * m_cfg.head_width_mm;
+    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);
@@ -509,7 +456,7 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head,
         if(zdiff > 0) {
             m_builder.add_pillar(head.id, bridgestart, r);
             m_builder.add_junction(bridgestart, r);
-            m_builder.add_bridge(bridgestart, bridgeend, head.r_back_mm);
+            m_builder.add_bridge(bridgestart, bridgeend, r);
         } else {
             m_builder.add_bridge(head.id, bridgeend);
         }
@@ -520,40 +467,6 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head,
     return true;
 }
 
-bool SupportTreeBuildsteps::search_pillar_and_connect(const Head &head)
-{
-    PointIndex spindex = m_pillar_index.guarded_clone();
-    
-    long nearest_id = ID_UNSET;
-    
-    Vec3d querypoint = head.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(querypoint(X), querypoint(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(head, nearest_id)) {
-                    nearest_id = ID_UNSET;    // continue searching
-                    spindex.remove(ne);       // without the current pillar
-                }
-            }
-        }
-    }
-    
-    return nearest_id >= 0;
-}
-
 void SupportTreeBuildsteps::create_ground_pillar(const Vec3d &jp,
                                                  const Vec3d &sourcedir,
                                                  double       radius,
@@ -565,9 +478,10 @@ void SupportTreeBuildsteps::create_ground_pillar(const Vec3d &jp,
     Vec3d  endp         = {jp(X), jp(Y), gndlvl};
     double sd           = m_cfg.pillar_base_safety_distance_mm;
     long   pillar_id    = ID_UNSET;
-    double min_dist     = sd + m_cfg.base_radius_mm + EPSILON;
+    bool   can_add_base = radius >= m_cfg.head_back_radius_mm;
+    double base_r       = can_add_base ? m_cfg.base_radius_mm : 0.;
+    double min_dist     = sd + base_r + EPSILON;
     double dist         = 0;
-    bool   can_add_base = true;
     bool   normal_mode  = true;
     
     // If in zero elevation mode and the pillar is too close to the model body,
@@ -612,7 +526,7 @@ void SupportTreeBuildsteps::create_ground_pillar(const Vec3d &jp,
         
         endp         = jp + std::get<0>(result.optimum) * dir;
         Vec3d pgnd   = {endp(X), endp(Y), gndlvl};
-        can_add_base = result.score > min_dist;
+        can_add_base = can_add_base && result.score > min_dist;
         
         double gnd_offs = m_mesh.ground_level_offset();
         auto abort_in_shame =
@@ -712,84 +626,85 @@ void SupportTreeBuildsteps::filter()
         auto [polar, azimuth] = dir_to_spheric(n);
         
         // skip if the tilt is not sane
-        if(polar >= PI - m_cfg.normal_cutoff_angle) {
-            
-            // We saturate the polar angle to 3pi/4
-            polar = std::max(polar, 3*PI / 4);
-            
-            // save the head (pinpoint) position
-            Vec3d hp = m_points.row(fidx);
-            
-            double w = m_cfg.head_width_mm +
-                       m_cfg.head_back_radius_mm +
-                       2*m_cfg.head_front_radius_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
-            EigenMesh3D::hit_result t
-                = pinhead_mesh_intersect(hp, // touching point
-                                         nn, // normal
-                                         pin_r,
-                                         m_cfg.head_back_radius_mm,
-                                         w);
+        if(polar < PI - m_cfg.normal_cutoff_angle) return;
             
-            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.
-                
-                StopCriteria stc;
-                stc.max_iterations = m_cfg.optimizer_max_iterations;
-                stc.relative_score_difference = m_cfg.optimizer_rel_score_diff;
-                stc.stop_score = w; // space greater than w is enough
-                GeneticOptimizer solver(stc);
-                solver.seed(0); // we want deterministic behavior
-                
-                auto oresult = solver.optimize_max(
-                    [this, pin_r, w, hp](double plr, double azm)
-                    {
-                        auto dir = spheric_to_dir(plr, azm).normalized();
-                        
-                        double score = pinhead_mesh_distance(
-                            hp, dir, pin_r, m_cfg.head_back_radius_mm, w);
-                        
-                        return score;
-                    },
-                    initvals(polar, azimuth), // start with what we have
-                    bound(3 * PI / 4, PI),    // Must not exceed the tilt limit
-                    bound(-PI, PI) // azimuth can be a full search
-                    );
-                
-                if(oresult.score > w) {
-                    polar = std::get<0>(oresult.optimum);
-                    azimuth = std::get<1>(oresult.optimum);
-                    nn = spheric_to_dir(polar, azimuth).normalized();
-                    t = EigenMesh3D::hit_result(oresult.score);
-                }
+        // We saturate the polar angle to 3pi/4
+        polar = std::max(polar, 3*PI / 4);
+
+        // save the head (pinpoint) position
+        Vec3d hp = m_points.row(fidx);
+
+        // The distance needed for a pinhead to not collide with model.
+        double w = m_cfg.head_width_mm +
+                   m_cfg.head_back_radius_mm +
+                   2*m_cfg.head_front_radius_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
+        EigenMesh3D::hit_result t
+            = pinhead_mesh_intersect(hp, // touching point
+                                     nn, // normal
+                                     pin_r,
+                                     m_cfg.head_back_radius_mm,
+                                     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.
+
+            StopCriteria stc;
+            stc.max_iterations = m_cfg.optimizer_max_iterations;
+            stc.relative_score_difference = m_cfg.optimizer_rel_score_diff;
+            stc.stop_score = w; // space greater than w is enough
+            GeneticOptimizer solver(stc);
+            solver.seed(0); // we want deterministic behavior
+
+            auto oresult = solver.optimize_max(
+                [this, pin_r, w, hp](double plr, double azm)
+                {
+                    auto dir = spheric_to_dir(plr, azm).normalized();
+
+                    double score = pinhead_mesh_intersect(
+                        hp, dir, pin_r, m_cfg.head_back_radius_mm, w).distance();
+
+                    return score;
+                },
+                initvals(polar, azimuth), // start with what we have
+                bound(3 * PI / 4, PI),    // Must not exceed the tilt limit
+                bound(-PI, PI) // azimuth can be a full search
+                );
+
+            if(oresult.score > w) {
+                polar = std::get<0>(oresult.optimum);
+                azimuth = std::get<1>(oresult.optimum);
+                nn = spheric_to_dir(polar, azimuth).normalized();
+                t = EigenMesh3D::hit_result(oresult.score);
             }
-            
-            // save the verified and corrected normal
-            m_support_nmls.row(fidx) = nn;
-            
-            if (t.distance() > w) {
-                // Check distance from ground, we might have zero elevation.
-                if (hp(Z) + w * nn(Z) < m_builder.ground_level) {
-                    addfn(m_iheadless, fidx);
-                } else {
-                    // mark the point for needing a head.
-                    addfn(m_iheads, fidx);
-                }
-            } else if (polar >= 3 * PI / 4) {
-                // Headless supports do not tilt like the headed ones
-                // so the normal should point almost to the ground.
+        }
+
+        // save the verified and corrected normal
+        m_support_nmls.row(fidx) = nn;
+
+        if (t.distance() > w) {
+            // Check distance from ground, we might have zero elevation.
+            if (hp(Z) + w * nn(Z) < m_builder.ground_level) {
                 addfn(m_iheadless, fidx);
+            } else {
+                // mark the point for needing a head.
+                addfn(m_iheads, fidx);
             }
+        } else if (polar >= 3 * PI / 4) {
+            // Headless supports do not tilt like the headed ones
+            // so the normal should point almost to the ground.
+            addfn(m_iheadless, fidx);
         }
+
     };
     
     ccr::enumerate(filtered_indices.begin(), filtered_indices.end(), filterfn);
@@ -811,6 +726,27 @@ void SupportTreeBuildsteps::add_pinheads()
             m_support_pts[i].pos.cast<double>() // displacement
             );
     }
+
+    for (unsigned i : m_iheadless) {
+        const auto R = double(m_support_pts[i].head_front_radius);
+
+        // The support point position on the mesh
+        Vec3d sph = m_support_pts[i].pos.cast<double>();
+
+        // Get an initial normal from the filtering step
+        Vec3d n = m_support_nmls.row(i);
+
+        // First we need to determine the available space for a mini pinhead.
+        // The goal is the move away from the model a little bit to make the
+        // contact point small as possible and avoid pearcing the model body.
+        double pin_space = std::min(2 * R, bridge_mesh_distance(sph, n, R, 0.));
+
+        if (pin_space <= 0) continue;
+
+        m_iheads.emplace_back(i);
+        m_builder.add_head(i, R, R, pin_space,
+                           m_cfg.head_penetration_mm, n, sph);
+    }
 }
 
 void SupportTreeBuildsteps::classify()
@@ -864,8 +800,6 @@ void SupportTreeBuildsteps::classify()
 
 void SupportTreeBuildsteps::routing_to_ground()
 {
-    const double pradius = m_cfg.head_back_radius_mm;
-    
     ClusterEl cl_centroids;
     cl_centroids.reserve(m_pillar_clusters.size());
     
@@ -931,7 +865,7 @@ void SupportTreeBuildsteps::routing_to_ground()
                 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, pradius, sidehead.id);
+                create_ground_pillar(pstart, sidehead.dir, sidehead.r_back_mm, sidehead.id);
             }
         }
     }
@@ -943,7 +877,7 @@ bool SupportTreeBuildsteps::connect_to_ground(Head &head, const Vec3d &dir)
     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);
+    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;
@@ -1041,6 +975,42 @@ bool SupportTreeBuildsteps::connect_to_model_body(Head &head)
     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 = 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)) {
+                    nearest_id = 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.
@@ -1054,18 +1024,18 @@ void SupportTreeBuildsteps::routing_to_model()
         auto& head = m_builder.head(idx);
         
         // Search nearby pillar
-        if(search_pillar_and_connect(head)) { head.transform(); return; }
+        if (search_pillar_and_connect(head)) { head.transform(); return; }
         
         // Cannot connect to nearby pillar. We will try to search for
         // a route to the ground.
-        if(connect_to_ground(head)) { head.transform(); return; }
+        if (connect_to_ground(head)) { head.transform(); 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;
+                << "Failed to route model facing support point. ID: " << idx;
         
         head.invalidate();
     });
@@ -1107,9 +1077,10 @@ void SupportTreeBuildsteps::interconnect_pillars()
         // 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, d](const PointIndexEl& e){
-            return distance(e.first, qp) < d;
+        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)
@@ -1288,37 +1259,54 @@ void SupportTreeBuildsteps::routing_headless()
     
     // We will sink the pins into the model surface for a distance of 1/3 of
     // the pin radius
-    for(unsigned i : m_iheadless) {
-        m_thr();
-        
-        const auto R = double(m_support_pts[i].head_front_radius);
-        const double HWIDTH_MM = std::min(R, m_cfg.head_penetration_mm);
-        
-        // Exact support position
-        Vec3d sph = m_support_pts[i].pos.cast<double>();
-        Vec3d n = m_support_nmls.row(i);   // mesh outward normal
-        Vec3d sp = sph - n * HWIDTH_MM;     // stick head start point
-        
-        Vec3d sj = sp + R * n;              // stick start point
-        
-        // This is only for checking
-        double idist = bridge_mesh_distance(sph, DOWN, R, true);
-        double realdist = ray_mesh_intersect(sj, DOWN).distance();
-        double dist = realdist;
-        
-        if (std::isinf(dist)) dist = sph(Z) - m_builder.ground_level;
-        
-        if(std::isnan(idist) || idist < 2*R || std::isnan(dist) || dist < 2*R) {
-            BOOST_LOG_TRIVIAL(warning) << "Can not find route for headless"
-                                       << " support stick at: "
-                                       << sj.transpose();
-            continue;
-        }
-        
-        bool use_endball = !std::isinf(realdist);
-        Vec3d ej = sj + (dist + HWIDTH_MM) * DOWN ;
-        m_builder.add_compact_bridge(sp, ej, n, R, use_endball);
-    }
+//    for(unsigned i : m_iheadless) {
+//        m_thr();
+
+//        const auto R = double(m_support_pts[i].head_front_radius);
+
+//        // The support point position on the mesh
+//        Vec3d sph = m_support_pts[i].pos.cast<double>();
+
+//        // Get an initial normal from the filtering step
+//        Vec3d n = m_support_nmls.row(i);
+
+//        // First we need to determine the available space for a mini pinhead.
+//        // The goal is the move away from the model a little bit to make the
+//        // contact point small as possible and avoid pearcing the model body.
+//        double pin_space = std::min(2 * R, bridge_mesh_distance(sph, n, R, 0.));
+
+//        if (pin_space <= 0) continue;
+
+//        auto &head = m_builder.add_head(i, R, R, pin_space,
+//                                        m_cfg.head_penetration_mm, n, sph);
+
+//        // collision check
+
+//            m_head_to_ground_scans[i] =
+//                bridge_mesh_intersect(head.junction_point(), DOWN, R);
+
+//        // Here the steps will be similar as in route_to_model step:
+//        // 1. Search for a nearby pillar, include other mini pillars
+
+//        // Search nearby pillar
+//        if (search_pillar_and_connect(head)) { head.transform(); continue; }
+
+//        if (std::isinf(m_head_to_ground_scans[i].distance())) {
+//            create_ground_pillar(head.junction_point(), head.dir, m_cfg.head_back_radius_mm, head.id);
+//        }
+
+//        // Cannot connect to nearby pillar. We will try to search for
+//        // a route to the ground.
+//        if (connect_to_ground(head)) { head.transform(); continue; }
+
+//        // No route to the ground, so connect to the model body as a last resort
+//        if (connect_to_model_body(head)) { continue; }
+
+//        BOOST_LOG_TRIVIAL(warning) << "Can not find route for headless"
+//                                   << " support stick at: "
+//                                   << sph.transpose();
+//        head.invalidate();
+//    }
 }
 
 }