SLA autosupports including islands

src/libslic3r/SLA/SLAAutoSupports.cpp

@@ -8,21 +8,41 @@
 #include "Point.hpp"
 
 #include <iostream>
-
+#include <random>
 
 namespace Slic3r {
-namespace SLAAutoSupports {
+
-SLAAutoSupports::SLAAutoSupports(ModelObject& mo, const SLAAutoSupports::Config& c)
+SLAAutoSupports::SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices, const std::vector<float>& heights, const Config& config)
-: m_model_object(mo), mesh(), m_config(c)
+: m_config(config), m_V(emesh.V), m_F(emesh.F)
-{}
+{
+    // FIXME: It might be safer to get rid of the rand() calls altogether, because it is probably
+    // not always thread-safe and can be slow if it is.
+    srand(time(NULL)); // rand() is used by igl::random_point_on_mesh
+
+    // Find all separate islands that will need support. The coord_t number denotes height
+    // of a point just below the mesh (so that we can later project the point precisely
+    // on the mesh by raycasting (done by igl) and not risking we will place the point inside).
+    std::vector<std::pair<ExPolygon, coord_t>> islands = find_islands(slices, heights);
+
+    // Uniformly cover each of the islands with support points.
+    for (const auto& island : islands) {
+        std::vector<Vec3d> points = uniformly_cover(island);
+        project_upward_onto_mesh(points);
+        m_output.insert(m_output.end(), points.begin(), points.end());
+    }
+
+    // We are done with the islands. Let's sprinkle the rest of the mesh.
+    // The function appends to m_output.
+    sprinkle_mesh(mesh);
+}
 
 
-float SLAAutoSupports::approximate_geodesic_distance(const Vec3f& p1, const Vec3f& p2, Vec3f& n1, Vec3f& n2)
+float SLAAutoSupports::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2)
 {
     n1.normalize();
     n2.normalize();
 
-    Vec3f v = (p2-p1);
+    Vec3d v = (p2-p1);
     v.normalize();
 
     float c1 = n1.dot(v);
@@ -35,15 +55,16 @@ float SLAAutoSupports::approximate_geodesic_distance(const Vec3f& p1, const Vec3
 }
 
 
-void SLAAutoSupports::generate()
+void SLAAutoSupports::sprinkle_mesh(const TriangleMesh& mesh)
 {
+    std::vector<Vec3d> points;
     // Loads the ModelObject raw_mesh and transforms it by first instance's transformation matrix (disregarding translation).
     // Instances only differ in z-rotation, so it does not matter which of them will be used for the calculation.
     // The supports point will be calculated on this mesh (so scaling ang vertical direction is correctly accounted for).
     // Results will be inverse-transformed to raw_mesh coordinates.
-    TriangleMesh mesh = m_model_object.raw_mesh();
+    //TriangleMesh mesh = m_model_object.raw_mesh();
-    Transform3d transformation_matrix = m_model_object.instances[0]->get_matrix(true/*dont_translate*/);
+    //Transform3d transformation_matrix = m_model_object.instances[0]->get_matrix(true/*dont_translate*/);
-    mesh.transform(transformation_matrix);
+    //mesh.transform(transformation_matrix);
 
     // Check that the object is thick enough to produce any support points
     BoundingBoxf3 bb = mesh.bounding_box();
@@ -51,30 +72,20 @@ void SLAAutoSupports::generate()
         return;
 
     // All points that we curretly have must be transformed too, so distance to them is correcly calculated.
-    for (Vec3f& point : m_model_object.sla_support_points)
+    //for (Vec3f& point : m_model_object.sla_support_points)
-        point = transformation_matrix.cast<float>() * point;
+    //    point = transformation_matrix.cast<float>() * point;
 
-    const stl_file& stl = mesh.stl;
-    Eigen::MatrixXf V;
-    Eigen::MatrixXi F;
-    V.resize(3 * stl.stats.number_of_facets, 3);
-    F.resize(stl.stats.number_of_facets, 3);
-    for (unsigned int i=0; i<stl.stats.number_of_facets; ++i) {
-        const stl_facet* facet = stl.facet_start+i;
-        V(3*i+0, 0) = facet->vertex[0](0); V(3*i+0, 1) = facet->vertex[0](1); V(3*i+0, 2) = facet->vertex[0](2);
-        V(3*i+1, 0) = facet->vertex[1](0); V(3*i+1, 1) = facet->vertex[1](1); V(3*i+1, 2) = facet->vertex[1](2);
-        V(3*i+2, 0) = facet->vertex[2](0); V(3*i+2, 1) = facet->vertex[2](1); V(3*i+2, 2) = facet->vertex[2](2);
-        F(i, 0) = 3*i+0;
-        F(i, 1) = 3*i+1;
-        F(i, 2) = 3*i+2;
-    }
 
     // In order to calculate distance to already placed points, we must keep know which facet the point lies on.
-    std::vector<Vec3f> facets_normals;
+    std::vector<Vec3d> facets_normals;
+
+    // Only points belonging to islands were added so far - they all lie on horizontal surfaces:
+    for (unsigned int i=0; i<m_output.size(); ++i)
+        facets_normals.push_back(Vec3d(0,0,-1));
 
     // The AABB hierarchy will be used to find normals of already placed points.
     // The points added automatically will just push_back the new normal on the fly.
-    igl::AABB<Eigen::MatrixXf,3> aabb;
+    /*igl::AABB<Eigen::MatrixXf,3> aabb;
     aabb.init(V, F);
     for (unsigned int i=0; i<m_model_object.sla_support_points.size(); ++i) {
         int facet_idx = 0;
@@ -86,61 +97,62 @@ void SLAAutoSupports::generate()
         Vec3f normal = a1.cross(a2);
         normal.normalize();
         facets_normals.push_back(normal);
-    }
+    }*/
 
     // New potential support point is randomly generated on the mesh and distance to all already placed points is calculated.
     // In case it is never smaller than certain limit (depends on the new point's facet normal), the point is accepted.
     // The process stops after certain number of points is refused in a row.
-    Vec3f point;
+    Vec3d point;
-    Vec3f normal;
+    Vec3d normal;
     int added_points = 0;
     int refused_points = 0;
     const int refused_limit = 30;
     // Angle at which the density reaches zero:
     const float threshold_angle = std::min(M_PI_2, M_PI_4 * acos(0.f/m_config.density_at_horizontal) / acos(m_config.density_at_45/m_config.density_at_horizontal));
 
-    srand(time(NULL)); // rand() is used by igl::random_point_on_mesh
     while (refused_points < refused_limit) {
         // Place a random point on the mesh and calculate corresponding facet's normal:
         Eigen::VectorXi FI;
-        Eigen::MatrixXf B;
+        Eigen::MatrixXd B;
-        igl::random_points_on_mesh(1, V, F, B, FI);
+        igl::random_points_on_mesh(1, m_V, m_F, B, FI);
-        point = B(0,0)*V.row(F(FI(0),0)) +
+        point = B(0,0)*m_V.row(m_F(FI(0),0)) +
-                B(0,1)*V.row(F(FI(0),1)) +
+                B(0,1)*m_V.row(m_F(FI(0),1)) +
-                B(0,2)*V.row(F(FI(0),2));
+                B(0,2)*m_V.row(m_F(FI(0),2));
         if (point(2) - bb.min(2) < m_config.minimal_z)
             continue;
 
-        Vec3f a1 = V.row(F(FI(0),1)) - V.row(F(FI(0),0));
+        Vec3d a1 = m_V.row(m_F(FI(0),1)) - m_V.row(m_F(FI(0),0));
-        Vec3f a2 = V.row(F(FI(0),2)) - V.row(F(FI(0),0));
+        Vec3d a2 = m_V.row(m_F(FI(0),2)) - m_V.row(m_F(FI(0),0));
         normal = a1.cross(a2);
         normal.normalize();
 
         // calculate angle between the normal and vertical:
-        float angle = angle_from_normal(normal);
+        float angle = angle_from_normal(normal.cast<float>());
         if (angle > threshold_angle)
             continue;
 
-        const float distance_limit = 1./(2.4*get_required_density(angle));
+        const float limit = distance_limit(angle);
         bool add_it = true;
-        for (unsigned int i=0; i<m_model_object.sla_support_points.size(); ++i) {
+        for (unsigned int i=0; i<points.size(); ++i) {
-            if (approximate_geodesic_distance(m_model_object.sla_support_points[i], point, facets_normals[i], normal) < distance_limit) {
+            if (approximate_geodesic_distance(points[i], point, facets_normals[i], normal) < limit) {
                 add_it = false;
                 ++refused_points;
                 break;
             }
         }
         if (add_it) {
-            m_model_object.sla_support_points.push_back(point);
+            points.push_back(point.cast<double>());
             facets_normals.push_back(normal);
             ++added_points;
             refused_points = 0;
         }
     }
 
+    m_output.insert(m_output.end(), points.begin(), points.end());
+
     // Now transform all support points to mesh coordinates:
-    for (Vec3f& point : m_model_object.sla_support_points)
+    //for (Vec3f& point : m_model_object.sla_support_points)
-        point = transformation_matrix.inverse().cast<float>() * point;
+    //    point = transformation_matrix.inverse().cast<float>() * point;
 }
 
 
@@ -153,13 +165,13 @@ float SLAAutoSupports::get_required_density(float angle) const
     return std::max(0.f, float(m_config.density_at_horizontal * cos(K*angle)));
 }
 
+float SLAAutoSupports::distance_limit(float angle) const
+{
+    return 1./(2.4*get_required_density(angle));
+}
 
 
-
+void SLAAutoSupports::output_expolygons(const ExPolygons& expolys, std::string filename) const
-
-
-
-void output_expolygons(const ExPolygons& expolys, std::string filename)
 {
     BoundingBox bb(Point(-30000000, -30000000), Point(30000000, 30000000));
     Slic3r::SVG svg_cummulative(filename, bb);
@@ -176,9 +188,9 @@ void output_expolygons(const ExPolygons& expolys, std::string filename)
     }
 }
 
-std::vector<Vec3d> find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights)
+std::vector<std::pair<ExPolygon, coord_t>> SLAAutoSupports::find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights) const
 {
-    std::vector<Vec3d> support_points_out;
+    std::vector<std::pair<ExPolygon, coord_t>> islands;
 
     struct PointAccessor {
         const Point* operator()(const Point &pt) const { return &pt; }
@@ -192,7 +204,6 @@ std::vector<Vec3d> find_islands(const std::vector<ExPolygons>& slices, const std
         std::string layer_num_str = std::string((i<10 ? "0" : "")) + std::string((i<100 ? "0" : "")) + std::to_string(i);
         output_expolygons(expolys_top, "top" + layer_num_str + ".svg");
         ExPolygons diff = diff_ex(expolys_top, expolys_bottom);
-        ExPolygons islands;
 
         output_expolygons(diff, "diff" + layer_num_str + ".svg");
 
@@ -223,22 +234,86 @@ std::vector<Vec3d> find_islands(const std::vector<ExPolygons>& slices, const std
                 }
 
             if (island) { // all points of the diff polygon are from the top slice
-                islands.push_back(polygon);
+                islands.push_back(std::make_pair(polygon, scale_(i!=0 ? heights[i-1] : heights[0]-(heights[1]-heights[0]))));
             }
             NO_ISLAND: ;// continue with next ExPolygon
         }
 
-        if (!islands.empty())
+        //if (!islands.empty())
-            output_expolygons(islands, "islands" + layer_num_str + ".svg");
+          //  output_expolygons(islands, "islands" + layer_num_str + ".svg");
-        for (const ExPolygon& island : islands) {
-            Point centroid = island.contour.centroid();
-            Vec3d centroid_d(centroid(0), centroid(1), scale_(i!=0 ? heights[i-1] : heights[0]-(heights[1]-heights[0]))) ;
-            support_points_out.push_back(unscale(centroid_d));
-        }
     }
 
-    return support_points_out;
+    return islands;
 }
 
-} // namespace SLAAutoSupports
+std::vector<Vec3d> SLAAutoSupports::uniformly_cover(const std::pair<ExPolygon, coord_t>& island)
+{
+    int num_of_points = std::max(1, (int)(island.first.area()*pow(SCALING_FACTOR, 2) * get_required_density(0)));
+
+    // In case there is just one point to place, we'll place it into the polygon's centroid (unless it lies in a hole).
+    if (num_of_points == 1) {
+        Point out(island.first.contour.centroid());
+        out(2) = island.second;
+
+        for (const auto& hole : island.first.holes)
+            if (hole.contains(out))
+                goto HOLE_HIT;
+        return std::vector<Vec3d>{unscale(out(0), out(1), out(2))};
+    }
+
+HOLE_HIT:
+    // In this case either the centroid lies in a hole, or there are multiple points
+    // to place. We will cover the island another way.
+    // For now we'll just place the points randomly not too close to the others.
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<> dis(0., 1.);
+
+    std::vector<Vec3d> island_new_points;
+    const BoundingBox& bb = get_extents(island.first);
+    const int refused_limit = 30;
+    int refused_points = 0;
+    while (refused_points < refused_limit) {
+        Point out(bb.min(0) + bb.size()(0)  * dis(gen),
+                  bb.min(1) + bb.size()(1)  * dis(gen)) ;
+        Vec3d unscaled_out = unscale(out(0), out(1), island.second);
+        bool add_it = true;
+
+        if (!island.first.contour.contains(out))
+            add_it = false;
+        else
+            for (const Polygon& hole : island.first.holes)
+                if (hole.contains(out))
+                    add_it = false;
+
+        if (add_it) {
+            for (const Vec3d& p : island_new_points) {
+                if ((p - unscaled_out).squaredNorm() < distance_limit(0)) {
+                    add_it = false;
+                    ++refused_points;
+                    break;
+                }
+            }
+        }
+        if (add_it) {
+            out(2) = island.second;
+            island_new_points.emplace_back(unscaled_out);
+        }
+    }
+    return island_new_points;
+}
+
+void SLAAutoSupports::project_upward_onto_mesh(std::vector<Vec3d>& points) const
+{
+    Vec3f dir(0., 0., 1.);
+    igl::Hit hit{0, 0, 0.f, 0.f, 0.f};
+    for (Vec3d& p : points) {
+        igl::ray_mesh_intersect(p.cast<float>(), dir, m_V, m_F, hit);
+        int fid = hit.id;
+        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<double>();
+    }
+}
+
+
 } // namespace Slic3r

src/libslic3r/SLA/SLAAutoSupports.hpp

@@ -3,16 +3,12 @@
 
 #include <libslic3r/Point.hpp>
 #include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/SLA/SLASupportTree.hpp>
 
 
 
 namespace Slic3r {
 
-class ModelObject;
-
-namespace SLAAutoSupports {
-
-
 class SLAAutoSupports {
 public:
     struct Config {
@@ -21,24 +17,30 @@ public:
             float minimal_z;
         };
 
-    SLAAutoSupports(ModelObject& mo, const SLAAutoSupports::Config& c);
+    SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices, const std::vector<float>& heights, const Config& config);
-    void generate();
+    const std::vector<Vec3d>& output() { return m_output; }
 
 private:
+    std::vector<Vec3d> m_output;
+    std::vector<Vec3d> m_normals;
     TriangleMesh mesh;
     static float angle_from_normal(const stl_normal& normal) { return acos((-normal.normalized())(2)); }
     float get_required_density(float angle) const;
-    static float approximate_geodesic_distance(const Vec3f& p1, const Vec3f& p2, Vec3f& n1, Vec3f& n2);
+    float distance_limit(float angle) const;
+    static float approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2);
+    std::vector<std::pair<ExPolygon, coord_t>> find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights) const;
+    void sprinkle_mesh(const TriangleMesh& mesh);
+    std::vector<Vec3d> uniformly_cover(const std::pair<ExPolygon, coord_t>& island);
+    void project_upward_onto_mesh(std::vector<Vec3d>& points) const;
+
+    void output_expolygons(const ExPolygons& expolys, std::string filename) const;
 
-    ModelObject& m_model_object;
     SLAAutoSupports::Config m_config;
+    const Eigen::MatrixXd& m_V;
+    const Eigen::MatrixXi& m_F;
 };
 
 
-std::vector<Vec3d> find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights);
-
-} // namespace SLAAutoSupports
-
 } // namespace Slic3r
 
 

src/libslic3r/SLAPrint.cpp

@@ -511,31 +511,33 @@ void SLAPrint::process()
     // In this step we check the slices, identify island and cover them with
     // support points. Then we sprinkle the rest of the mesh.
     auto support_points = [this, ilh](SLAPrintObject& po) {
-        // find islands to support
+        const ModelObject& mo = *po.m_model_object;
-        double lh = po.m_config.layer_height.getFloat();
-        std::vector<float> heights = calculate_heights(po.transformed_mesh().bounding_box(), po.get_elevation(), ilh, lh);
-        //SLAAutoSupports auto_supports(po.get_model_slices(), heights, *po.m_model_object);
-        std::vector<Vec3d> points = SLAAutoSupports::find_islands(po.get_model_slices(), heights);
-
-        // TODO:
-        
-        // create mesh in igl format
-
-        // cover the islands with points, use igl to get precise z coordinate
-
-        // sprinkle the mesh with points (SLAAutoSupports::generate())
-
-
-        /*for (const auto& p: points)
-            std::cout << p(0) << " " << p(1) << " " << p(2) << std::endl;
-        std::cout << std::endl;
-        */
-        //for (auto& p: points)
-          //  p = po.trafo().inverse() * p;
-
         po.m_supportdata.reset(new SLAPrintObject::SupportData());
         po.m_supportdata->emesh = sla::to_eigenmesh(po.transformed_mesh());
-        po.m_supportdata->support_points = sla::to_point_set(points);
+
+        // If there are no points on the front-end, we will do the autoplacement.
+        // Otherwise we will just blindly copy the frontend data into the backend cache.
+        if(mo.sla_support_points.empty()) {
+            // calculate heights of slices (slices are calculated already)
+            double lh = po.m_config.layer_height.getFloat();
+            std::vector<float> heights = calculate_heights(po.transformed_mesh().bounding_box(), po.get_elevation(), ilh, lh);
+
+            SLAAutoSupports::Config config;
+            const SLAPrintObjectConfig& cfg = po.config();
+            config.density_at_horizontal = cfg.support_density_at_horizontal / 10000.f;
+            config.density_at_45 = cfg.support_density_at_45 / 10000.f;
+            config.minimal_z = cfg.support_minimal_z;
+
+            // Construction of this object does the calculation.
+            SLAAutoSupports auto_supports(po.transformed_mesh(), po.m_supportdata->emesh, po.get_model_slices(), heights, config);
+            // Now let's extract the result.
+            const std::vector<Vec3d>& points = auto_supports.output();
+            po.m_supportdata->support_points = sla::to_point_set(points);
+        }
+        else {
+            // There are some points on the front-end, no calculation will be done.
+            po.m_supportdata->support_points = sla::to_point_set(po.transformed_support_points());
+        }
     };
 
     // In this step we create the supports
@@ -1131,6 +1133,11 @@ const std::vector<ExPolygons> EMPTY_SLICES;
 const TriangleMesh EMPTY_MESH;
 }
 
+const Eigen::MatrixXd& SLAPrintObject::get_support_points() const
+{
+    return m_supportdata->support_points;
+}
+
 const std::vector<ExPolygons> &SLAPrintObject::get_support_slices() const
 {
     // assert(is_step_done(slaposSliceSupports));

src/libslic3r/SLAPrint.hpp

@@ -90,6 +90,9 @@ public:
     const std::vector<ExPolygons>& get_model_slices() const;
     const std::vector<ExPolygons>& get_support_slices() const;
 
+    // This method returns the support points of this SLAPrintObject.
+    const Eigen::MatrixXd& get_support_points() const;
+
     // An index record referencing the slices
     // (get_model_slices(), get_support_slices()) where the keys are the height
     // levels of the model in scaled-clipper coordinates. The levels correspond

src/slic3r/GUI/GLCanvas3D.hpp

@@ -1145,6 +1145,7 @@ private:
 
     static std::vector<float> _parse_colors(const std::vector<std::string>& colors);
 
+public:
     const Print* fff_print() const;
     const SLAPrint* sla_print() const;
 };

src/slic3r/GUI/GLGizmo.cpp

@@ -8,6 +8,7 @@
 #include "libslic3r/Geometry.hpp"
 #include "libslic3r/Utils.hpp"
 #include "libslic3r/SLA/SLASupportTree.hpp"
+#include "libslic3r/SLAPrint.hpp"
 
 #include <cstdio>
 #include <numeric>
@@ -1786,6 +1787,18 @@ void GLGizmoSlaSupports::set_sla_support_data(ModelObject* model_object, const G
     {
         if (is_mesh_update_necessary())
             update_mesh();
+
+        // If there are no points, let's ask the backend if it calculated some.
+        if (model_object->sla_support_points.empty() && m_parent.sla_print()->is_step_done(slaposSupportPoints)) {
+            for (const SLAPrintObject* po : m_parent.sla_print()->objects()) {
+                if (po->model_object()->id() == model_object->id()) {
+                    const Eigen::MatrixXd& points = po->get_support_points();
+                    for (unsigned int i=0; i<points.rows();++i)
+                        model_object->sla_support_points.push_back(Vec3f(po->trafo().inverse().cast<float>() * Vec3f(points(i,0), points(i,1), points(i,2))));
+                        break;
+                }
+            }
+        }
     }
 }
 #else
@@ -2170,6 +2183,9 @@ void GLGizmoSlaSupports::render_tooltip_texture() const {
 #if ENABLE_IMGUI
 void GLGizmoSlaSupports::on_render_input_window(float x, float y, const GLCanvas3D::Selection& selection)
 {
+    bool first_run = true; // This is a hack to redraw the button when all points are removed,
+                           // so it is not delayed until the background process finishes.
+RENDER_AGAIN:
     m_imgui->set_next_window_pos(x, y, ImGuiCond_Always);
     m_imgui->set_next_window_bg_alpha(0.5f);
     m_imgui->begin(on_get_name(), ImGuiWindowFlags_NoMove | ImGuiWindowFlags_NoResize | ImGuiWindowFlags_NoCollapse);
@@ -2184,22 +2200,11 @@ void GLGizmoSlaSupports::on_render_input_window(float x, float y, const GLCanvas
 
     m_imgui->end();
 
-    if (remove_all_clicked)
+    if (remove_all_clicked) {
         delete_current_grabber(true);
-
+        if (first_run) {
-    if (generate) {
+            first_run = false;
-        const DynamicPrintConfig& cfg = *wxGetApp().get_tab(Preset::TYPE_SLA_PRINT)->get_config();
+            goto RENDER_AGAIN;
-        SLAAutoSupports::SLAAutoSupports::Config config;
-        config.density_at_horizontal = cfg.opt_int("support_density_at_horizontal") / 10000.f;
-        config.density_at_45 = cfg.opt_int("support_density_at_45") / 10000.f;
-        config.minimal_z = cfg.opt_float("support_minimal_z");
-
-        SLAAutoSupports::SLAAutoSupports sas(*m_model_object, config);
-        sas.generate();
-        m_grabbers.clear();
-        for (const Vec3f& point : m_model_object->sla_support_points) {
-            m_grabbers.push_back(Grabber());
-            m_grabbers.back().center = point.cast<double>();
         }
     }