Add dedicated tests for support point generation

src/libslic3r/BoundingBox.hpp

@@ -192,6 +192,20 @@ inline BoundingBox3 scaled(const BoundingBoxf3 &bb) { return {scaled(bb.min), sc
 inline BoundingBoxf unscaled(const BoundingBox &bb) { return {unscaled(bb.min), unscaled(bb.max)}; }
 inline BoundingBoxf3 unscaled(const BoundingBox3 &bb) { return {unscaled(bb.min), unscaled(bb.max)}; }
 
+template<class Tout, class Tin>
+auto cast(const BoundingBoxBase<Tin> &b)
+{
+    return BoundingBoxBase<Vec<3, Tout>>{b.min.template cast<Tout>(),
+                                         b.max.template cast<Tout>()};
+}
+
+template<class Tout, class Tin>
+auto cast(const BoundingBox3Base<Tin> &b)
+{
+    return BoundingBox3Base<Vec<3, Tout>>{b.min.template cast<Tout>(),
+                                          b.max.template cast<Tout>()};
+}
+
 } // namespace Slic3r
 
 // Serialization through the Cereal library

src/libslic3r/Line.cpp

@@ -33,24 +33,6 @@ bool Line::intersection_infinite(const Line &other, Point* point) const
     return true;
 }
 
-// Distance to the closest point of line.
-double Line::distance_to_squared(const Point &point, const Point &a, const Point &b)
-{
-    const Vec2d   v  = (b - a).cast<double>();
-    const Vec2d   va = (point  - a).cast<double>();
-    const double  l2 = v.squaredNorm();  // avoid a sqrt
-    if (l2 == 0.0) 
-        // a == b case
-        return va.squaredNorm();
-    // Consider the line extending the segment, parameterized as a + t (b - a).
-    // We find projection of this point onto the line. 
-    // It falls where t = [(this-a) . (b-a)] / |b-a|^2
-    const double t = va.dot(v) / l2;
-    if (t < 0.0)      return va.squaredNorm();  // beyond the 'a' end of the segment
-    else if (t > 1.0) return (point - b).cast<double>().squaredNorm();  // beyond the 'b' end of the segment
-    return (t * v - va).squaredNorm();
-}
-
 double Line::perp_distance_to(const Point &point) const
 {
     const Line  &line = *this;

src/libslic3r/Line.hpp

@@ -18,6 +18,35 @@ typedef std::vector<ThickLine> ThickLines;
 
 Linef3 transform(const Linef3& line, const Transform3d& t);
 
+namespace line_alg {
+
+// Distance to the closest point of line.
+template<class L, class T, int N>
+double distance_to_squared(const L &line, const Vec<N, T> &point)
+{
+    const Vec<N, double>  v  = line.vector().template cast<double>();
+    const Vec<N, double>  va = (point  - line.a).template cast<double>();
+    const double  l2 = v.squaredNorm();  // avoid a sqrt
+    if (l2 == 0.0)
+        // a == b case
+        return va.squaredNorm();
+    // Consider the line extending the segment, parameterized as a + t (b - a).
+    // We find projection of this point onto the line.
+    // It falls where t = [(this-a) . (b-a)] / |b-a|^2
+    const double t = va.dot(v) / l2;
+    if (t < 0.0)      return va.squaredNorm();  // beyond the 'a' end of the segment
+    else if (t > 1.0) return (point - line.b).template cast<double>().squaredNorm();  // beyond the 'b' end of the segment
+    return (t * v - va).squaredNorm();
+}
+
+template<class L, class T, int N>
+double distance_to(const L &line, const Vec<N, T> &point)
+{
+    return std::sqrt(distance_to_squared(line, point));
+}
+
+} // namespace line_alg
+
 class Line
 {
 public:
@@ -47,7 +76,7 @@ public:
     // Clip a line with a bounding box. Returns false if the line is completely outside of the bounding box.
 	bool   clip_with_bbox(const BoundingBox &bbox);
 
-    static double distance_to_squared(const Point &point, const Point &a, const Point &b);
+    static inline double distance_to_squared(const Point &point, const Point &a, const Point &b) { return line_alg::distance_to_squared(Line{a, b}, Vec<2, coord_t>{point}); }
     static double distance_to(const Point &point, const Point &a, const Point &b) { return sqrt(distance_to_squared(point, a, b)); }
 
     Point a;

src/libslic3r/Point.hpp

@@ -88,6 +88,8 @@ inline std::string to_string(const Vec3d   &pt) { return std::string("[") + std:
 std::vector<Vec3f> transform(const std::vector<Vec3f>& points, const Transform3f& t);
 Pointf3s transform(const Pointf3s& points, const Transform3d& t);
 
+template<int N, class T> using Vec = Eigen::Matrix<T,  N, 1, Eigen::DontAlign>;
+
 class Point : public Vec2crd
 {
 public:

src/libslic3r/SLA/Hollowing.cpp

@@ -273,4 +273,13 @@ void cut_drainholes(std::vector<ExPolygons> & obj_slices,
         obj_slices[i] = diff_ex(obj_slices[i], hole_slices[i]);
 }
 
+void hollow_mesh(TriangleMesh &mesh, const HollowingConfig &cfg)
+{
+    std::unique_ptr<Slic3r::TriangleMesh> inter_ptr =
+            Slic3r::sla::generate_interior(mesh);
+
+    if (inter_ptr) mesh.merge(*inter_ptr);
+    mesh.require_shared_vertices();
+}
+
 }} // namespace Slic3r::sla

src/libslic3r/SLA/Hollowing.hpp

@@ -62,6 +62,8 @@ std::unique_ptr<TriangleMesh> generate_interior(const TriangleMesh &mesh,
                                                 const HollowingConfig &  = {},
                                                 const JobController &ctl = {});
 
+void hollow_mesh(TriangleMesh &mesh, const HollowingConfig &cfg);
+
 void cut_drainholes(std::vector<ExPolygons> & obj_slices,
                     const std::vector<float> &slicegrid,
                     float                     closing_radius,

tests/sla_print/CMakeLists.txt

@@ -1,7 +1,8 @@
 get_filename_component(_TEST_NAME ${CMAKE_CURRENT_LIST_DIR} NAME)
-add_executable(${_TEST_NAME}_tests ${_TEST_NAME}_tests_main.cpp 
+add_executable(${_TEST_NAME}_tests ${_TEST_NAME}_tests_main.cpp
-    sla_print_tests.cpp 
+    sla_print_tests.cpp
     sla_test_utils.hpp sla_test_utils.cpp sla_treebuilder_tests.cpp
+    sla_supptgen_tests.cpp
     sla_raycast_tests.cpp)
 target_link_libraries(${_TEST_NAME}_tests test_common libslic3r)
 set_property(TARGET ${_TEST_NAME}_tests PROPERTY FOLDER "tests")

tests/sla_print/sla_supptgen_tests.cpp

@@ -0,0 +1,148 @@
+#include <catch2/catch.hpp>
+#include <test_utils.hpp>
+
+#include <libslic3r/ExPolygon.hpp>
+#include <libslic3r/BoundingBox.hpp>
+
+#include "sla_test_utils.hpp"
+
+namespace Slic3r { namespace sla {
+
+TEST_CASE("Overhanging point should be supported", "[SupGen]") {
+
+    // Pyramid with 45 deg slope
+    TriangleMesh mesh = make_pyramid(10.f, 10.f);
+    mesh.rotate_y(PI);
+    mesh.require_shared_vertices();
+    mesh.WriteOBJFile("Pyramid.obj");
+
+    sla::SupportPoints pts = calc_support_pts(mesh);
+
+    // The overhang, which is the upside-down pyramid's edge
+    Vec3f overh{0., 0., -10.};
+
+    REQUIRE(!pts.empty());
+
+    float dist = (overh - pts.front().pos).norm();
+
+    for (const auto &pt : pts)
+        dist = std::min(dist, (overh - pt.pos).norm());
+
+    // Should require exactly one support point at the overhang
+    REQUIRE(pts.size() > 0);
+    REQUIRE(dist < 1.f);
+}
+
+double min_point_distance(const sla::SupportPoints &pts)
+{
+    sla::PointIndex index;
+
+    for (size_t i = 0; i < pts.size(); ++i)
+        index.insert(pts[i].pos.cast<double>(), i);
+
+    auto d = std::numeric_limits<double>::max();
+    index.foreach([&d, &index](const sla::PointIndexEl &el) {
+        auto res = index.nearest(el.first, 2);
+        for (const sla::PointIndexEl &r : res)
+            if (r.second != el.second)
+                d = std::min(d, (el.first - r.first).norm());
+    });
+
+    return d;
+}
+
+TEST_CASE("Overhanging horizontal surface should be supported", "[SupGen]") {
+    double width = 10., depth = 10., height = 1.;
+
+    TriangleMesh mesh = make_cube(width, depth, height);
+    mesh.translate(0., 0., 5.); // lift up
+    mesh.require_shared_vertices();
+    mesh.WriteOBJFile("Cuboid.obj");
+
+    sla::SupportPointGenerator::Config cfg;
+    sla::SupportPoints pts = calc_support_pts(mesh, cfg);
+
+    double mm2 = width * depth;
+
+    REQUIRE(!pts.empty());
+    REQUIRE(pts.size() * cfg.support_force() > mm2 * cfg.tear_pressure());
+    REQUIRE(min_point_distance(pts) >= cfg.minimal_distance);
+}
+
+template<class M> auto&& center_around_bb(M &&mesh)
+{
+    auto bb = mesh.bounding_box();
+    mesh.translate(-bb.center().template cast<float>());
+
+    return std::forward<M>(mesh);
+}
+
+TEST_CASE("Overhanging edge should be supported", "[SupGen]") {
+    float width = 10.f, depth = 10.f, height = 5.f;
+
+    TriangleMesh mesh = make_prism(width, depth, height);
+    mesh.rotate_y(PI); // rotate on its back
+    mesh.translate(0., 0., height);
+    mesh.require_shared_vertices();
+    mesh.WriteOBJFile("Prism.obj");
+
+    sla::SupportPointGenerator::Config cfg;
+    sla::SupportPoints pts = calc_support_pts(mesh, cfg);
+
+    REQUIRE(min_point_distance(pts) >= cfg.minimal_distance);
+
+    Linef3 overh{ {0.f, -depth / 2.f, 0.f}, {0.f, depth / 2.f, 0.f}};
+
+    // Get all the points closer that 1 mm to the overhanging edge:
+    sla::SupportPoints overh_pts; overh_pts.reserve(pts.size());
+
+    std::copy_if(pts.begin(), pts.end(), std::back_inserter(overh_pts),
+                 [&overh](const sla::SupportPoint &pt){
+                     return line_alg::distance_to(overh, Vec3d{pt.pos.cast<double>()}) < 1.;
+                 });
+
+    REQUIRE(overh_pts.size() * cfg.support_force() > overh.length() * cfg.tear_pressure());
+    REQUIRE(min_point_distance(pts) >= cfg.minimal_distance);
+}
+
+// FIXME: Not working yet
+//TEST_CASE("Hollowed cube should be supported from the inside", "[SupGen][Hollowed]") {
+//    TriangleMesh mesh = make_cube(20., 20., 20.);
+
+//    hollow_mesh(mesh, HollowingConfig{});
+
+//    mesh.WriteOBJFile("cube_hollowed.obj");
+
+//    auto bb = mesh.bounding_box();
+//    auto h  = float(bb.max.z() - bb.min.z());
+//    Vec3f mv = bb.center().cast<float>() - Vec3f{0.f, 0.f, 0.5f * h};
+//    mesh.translate(-mv);
+//    mesh.require_shared_vertices();
+
+//    sla::SupportPointGenerator::Config cfg;
+//    sla::SupportPoints pts = calc_support_pts(mesh, cfg);
+//    sla::remove_bottom_points(pts, mesh.bounding_box().min.z() + EPSILON);
+
+//    REQUIRE(!pts.empty());
+//}
+
+TEST_CASE("Two parallel plates should be supported", "[SupGen][Hollowed]")
+{
+    double width = 20., depth = 20., height = 1.;
+
+    TriangleMesh mesh = center_around_bb(make_cube(width + 5., depth + 5., height));
+    TriangleMesh mesh_high = center_around_bb(make_cube(width, depth, height));
+    mesh_high.translate(0., 0., 10.); // lift up
+    mesh.merge(mesh_high);
+    mesh.require_shared_vertices();
+
+    mesh.WriteOBJFile("parallel_plates.obj");
+
+    sla::SupportPointGenerator::Config cfg;
+    sla::SupportPoints pts = calc_support_pts(mesh, cfg);
+    sla::remove_bottom_points(pts, mesh.bounding_box().min.z() + EPSILON);
+
+    REQUIRE(!pts.empty());
+}
+
+}} // namespace Slic3r::sla

tests/sla_print/sla_test_utils.cpp

@@ -411,3 +411,71 @@ double predict_error(const ExPolygon &p, const sla::RasterBase::PixelDim &pd)
     
     return error;
 }
+
+
+// Make a 3D pyramid
+TriangleMesh make_pyramid(float base, float height)
+{
+    float a = base / 2.f;
+
+    TriangleMesh mesh(
+        {
+            {-a, -a, 0}, {a, -a, 0}, {a, a, 0},
+            {-a, a, 0}, {0.f, 0.f, height}
+        },
+        {
+            {0, 1, 2},
+            {0, 2, 3},
+            {0, 1, 4},
+            {1, 2, 4},
+            {2, 3, 4},
+            {3, 0, 4}
+        });
+
+    mesh.repair();
+
+    return mesh;
+}
+
+    TriangleMesh make_prism(double width, double length, double height)
+{
+    // We need two upward facing triangles
+
+        double x = width / 2., y = length / 2.;
+
+        TriangleMesh mesh(
+            {
+                {-x, -y, 0.}, {x, -y, 0.}, {0., -y, height},
+                {-x, y, 0.}, {x, y, 0.}, {0., y, height},
+                },
+            {
+                {0, 1, 2}, // side 1
+                {4, 3, 5}, // side 2
+                {1, 4, 2}, {2, 4, 5}, // roof 1
+                {0, 2, 5}, {0, 5, 3}, // roof 2
+                {3, 4, 1}, {3, 1, 0} // bottom
+            });
+
+    return mesh;
+}
+
+sla::SupportPoints calc_support_pts(
+    const TriangleMesh &                      mesh,
+    const sla::SupportPointGenerator::Config &cfg)
+{
+    // Prepare the slice grid and the slices
+    std::vector<ExPolygons> slices;
+    auto                    bb      = cast<float>(mesh.bounding_box());
+    std::vector<float>      heights = grid(bb.min.z(), bb.max.z(), 0.1f);
+    slice_mesh(mesh, heights, slices, CLOSING_RADIUS, [] {});
+
+    // Prepare the support point calculator
+    sla::IndexedMesh emesh{mesh};
+    sla::SupportPointGenerator spgen{emesh, cfg, []{}, [](int){}};
+
+    // Calculate the support points
+    spgen.seed(0);
+    spgen.execute(slices, heights);
+
+    return spgen.output();
+}

tests/sla_print/sla_test_utils.hpp

@@ -185,4 +185,13 @@ long raster_pxsum(const sla::RasterGrayscaleAA &raster);
 
 double predict_error(const ExPolygon &p, const sla::RasterBase::PixelDim &pd);
 
+// Make a 3D pyramid
+TriangleMesh make_pyramid(float base, float height);
+
+TriangleMesh make_prism(double width, double length, double height);
+
+sla::SupportPoints calc_support_pts(
+    const TriangleMesh &                      mesh,
+    const sla::SupportPointGenerator::Config &cfg = {});
+
 #endif // SLA_TEST_UTILS_HPP