diff --git a/src/libslic3r/SLA/SLAAutoSupports.cpp b/src/libslic3r/SLA/SLAAutoSupports.cpp index c87f985dac..4389fd71f2 100644 --- a/src/libslic3r/SLA/SLAAutoSupports.cpp +++ b/src/libslic3r/SLA/SLAAutoSupports.cpp @@ -9,6 +9,8 @@ #include "SVG.hpp" #include "Point.hpp" #include "ClipperUtils.hpp" +#include "Tesselate.hpp" +#include "libslic3r.h" #include #include @@ -121,15 +123,15 @@ void SLAAutoSupports::process(const std::vector& slices, const std:: Structure& top = structures_new.back(); //FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands. // At least it is now using a bounding box check for pre-filtering. - for (Structure& bottom : structures_old) - if (top.overlaps(bottom)) { - top.structures_below.push_back(&bottom); + for (Structure& bottom : structures_old) + if (top.overlaps(bottom)) { + top.structures_below.push_back(&bottom); float centroids_dist = (bottom.centroid - top.centroid).norm(); // Penalization resulting from centroid offset: -// bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, (1600.f * layer_height) * centroids_dist * centroids_dist / bottom.area)); +// bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, (1600.f * layer_height) * centroids_dist * centroids_dist / bottom.area)); bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, 80.f * centroids_dist * centroids_dist / bottom.area)); // Penalization resulting from increasing polygon area: - bottom.supports_force *= std::min(1.f, 20.f * bottom.area / top.area); + bottom.supports_force *= std::min(1.f, 20.f * bottom.area / top.area); } } @@ -154,7 +156,7 @@ void SLAAutoSupports::process(const std::vector& slices, const std:: else // Let's see if there's anything that overlaps enough to need supports: // What we now have in polygons needs support, regardless of what the forces are, so we can add them. - for (const ExPolygon& p : diff_ex(to_polygons(*s.polygon), offset(s.expolygons_below(), between_layers_offset))) + for (const ExPolygon& p : diff_ex(to_polygons(*s.polygon), offset(s.expolygons_below(), between_layers_offset))) //FIXME is it an island point or not? Vojtech thinks it is. uniformly_cover(p, s); } @@ -163,12 +165,12 @@ void SLAAutoSupports::process(const std::vector& slices, const std:: for (Structure& s : structures_new) { // Areas not supported by the areas below. ExPolygons e = diff_ex(to_polygons(*s.polygon), s.polygons_below()); - float e_area = 0.f; - for (const ExPolygon &ex : e) - e_area += float(ex.area()); + float e_area = 0.f; + for (const ExPolygon &ex : e) + e_area += float(ex.area()); // Penalization resulting from large diff from the last layer: // s.supports_force /= std::max(1.f, (layer_height / 0.3f) * e_area / s.area); - s.supports_force /= std::max(1.f, 0.17f * (e_area * float(SCALING_FACTOR * SCALING_FACTOR)) / s.area); + s.supports_force /= std::max(1.f, 0.17f * (e_area * float(SCALING_FACTOR * SCALING_FACTOR)) / s.area); if (s.area * m_config.tear_pressure > s.supports_force) { //FIXME Don't calculate area inside the compare function! @@ -196,64 +198,204 @@ void SLAAutoSupports::process(const std::vector& slices, const std:: } } +std::vector sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng) +{ + // Triangulate the polygon with holes into triplets of 3D points. + std::vector triangles = Slic3r::triangulate_expolygon_2f(expoly); + + std::vector out; + if (! triangles.empty()) + { + // Calculate area of each triangle. + std::vector areas; + areas.reserve(triangles.size() / 3); + for (size_t i = 0; i < triangles.size(); ) { + const Vec2f &a = triangles[i ++]; + const Vec2f v1 = triangles[i ++] - a; + const Vec2f v2 = triangles[i ++] - a; + areas.emplace_back(0.5f * std::abs(cross2(v1, v2))); + if (i != 3) + // Prefix sum of the areas. + areas.back() += areas[areas.size() - 2]; + } + + size_t num_samples = size_t(ceil(areas.back() * samples_per_mm2)); + std::uniform_real_distribution<> random_triangle(0., double(areas.back())); + std::uniform_real_distribution<> random_float(0., 1.); + for (size_t i = 0; i < num_samples; ++ i) { + double r = random_triangle(rng); + size_t idx_triangle = std::min(std::upper_bound(areas.begin(), areas.end(), (float)r) - areas.begin(), areas.size() - 1) * 3; + // Select a random point on the triangle. + double u = float(sqrt(random_float(rng))); + double v = float(random_float(rng)); + const Vec2f &a = triangles[idx_triangle ++]; + const Vec2f &b = triangles[idx_triangle++]; + const Vec2f &c = triangles[idx_triangle]; + const Vec2f x = a * (1.f - u) + b * (u * (1.f - v)) + c * (v * u); + out.emplace_back(x); + } + } + return out; +} + +std::vector sample_expolygon_with_boundary(const ExPolygon &expoly, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng) +{ + std::vector out = sample_expolygon(expoly, samples_per_mm2, rng); + double point_stepping_scaled = scale_(1.f) / samples_per_mm_boundary; + for (size_t i_contour = 0; i_contour <= expoly.holes.size(); ++ i_contour) { + const Polygon &contour = (i_contour == 0) ? expoly.contour : expoly.holes[i_contour - 1]; + const Points pts = contour.equally_spaced_points(point_stepping_scaled); + for (size_t i = 0; i < pts.size(); ++ i) + out.emplace_back(unscale(pts[i].x()), unscale(pts[i].y())); + } + return out; +} + +std::vector poisson_disk_from_samples(const std::vector &raw_samples, float radius) +{ + Vec2f corner_min(FLT_MAX, FLT_MAX); + for (const Vec2f &pt : raw_samples) { + corner_min.x() = std::min(corner_min.x(), pt.x()); + corner_min.y() = std::min(corner_min.y(), pt.y()); + } + + // Assign the raw samples to grid cells, sort the grid cells lexicographically. + struct RawSample { + Vec2f coord; + Vec2i cell_id; + }; + std::vector raw_samples_sorted; + RawSample sample; + for (const Vec2f &pt : raw_samples) { + sample.coord = pt; + sample.cell_id = ((pt - corner_min) / radius).cast(); + raw_samples_sorted.emplace_back(sample); + } + std::sort(raw_samples_sorted.begin(), raw_samples_sorted.end(), [](const RawSample &lhs, const RawSample &rhs) + { return lhs.cell_id.x() < rhs.cell_id.x() || (lhs.cell_id.x() == rhs.cell_id.x() && lhs.cell_id.y() < rhs.cell_id.y()); }); + + struct PoissonDiskGridEntry { + // Resulting output sample points for this cell: + enum { + max_positions = 4 + }; + Vec2f poisson_samples[max_positions]; + int num_poisson_samples = 0; + + // Index into raw_samples: + int first_sample_idx; + int sample_cnt; + }; + + struct CellIDHash { + std::size_t operator()(const Vec2i &cell_id) { + return std::hash()(cell_id.x()) ^ std::hash()(cell_id.y() * 593); + } + }; + + // Map from cell IDs to hash_data. Each hash_data points to the range in raw_samples corresponding to that cell. + // (We could just store the samples in hash_data. This implementation is an artifact of the reference paper, which + // is optimizing for GPU acceleration that we haven't implemented currently.) + typedef std::unordered_map Cells; + std::unordered_map cells; + { + Cells::iterator last_cell_id_it; + Vec2i last_cell_id(-1, -1); + for (int i = 0; i < raw_samples_sorted.size(); ++ i) { + const RawSample &sample = raw_samples_sorted[i]; + if (sample.cell_id == last_cell_id) { + // This sample is in the same cell as the previous, so just increase the count. Cells are + // always contiguous, since we've sorted raw_samples_sorted by cell ID. + ++ last_cell_id_it->second.sample_cnt; + } else { + // This is a new cell. + PoissonDiskGridEntry data; + data.first_sample_idx = i; + data.sample_cnt = 1; + auto result = cells.insert({sample.cell_id, data}); + last_cell_id = sample.cell_id; + last_cell_id_it = result.first; + } + } + } + + const int max_trials = 5; + const float radius_squared = radius * radius; + for (int trial = 0; trial < max_trials; ++ trial) { + // Create sample points for each entry in cells. + for (auto &it : cells) { + const Vec2i &cell_id = it.first; + PoissonDiskGridEntry &cell_data = it.second; + // This cell's raw sample points start at first_sample_idx. On trial 0, try the first one. On trial 1, try first_sample_idx + 1. + int next_sample_idx = cell_data.first_sample_idx + trial; + if (trial >= cell_data.sample_cnt) + // There are no more points to try for this cell. + continue; + const RawSample &candidate = raw_samples_sorted[next_sample_idx]; + // See if this point conflicts with any other points in this cell, or with any points in + // neighboring cells. Note that it's possible to have more than one point in the same cell. + bool conflict = false; + for (int i = -1; i < 2 && ! conflict; ++ i) { + for (int j = -1; j < 2; ++ j) { + const auto &it_neighbor = cells.find(cell_id + Vec2i(i, j)); + if (it_neighbor != cells.end()) { + const PoissonDiskGridEntry &neighbor = it_neighbor->second; + for (int i_sample = 0; i_sample < neighbor.num_poisson_samples; ++ i_sample) + if ((neighbor.poisson_samples[i_sample] - candidate.coord).squaredNorm() < radius_squared) { + conflict = true; + break; + } + } + } + } + if (! conflict) { + // Store the new sample. + assert(cell_data.num_poisson_samples < cell_data.max_positions); + if (cell_data.num_poisson_samples < cell_data.max_positions) + cell_data.poisson_samples[cell_data.num_poisson_samples ++] = candidate.coord; + } + } + } + + // Copy the results to the output. + std::vector out; + for (const auto it : cells) + for (int i = 0; i < it.second.num_poisson_samples; ++ i) + out.emplace_back(it.second.poisson_samples[i]); + return out; +} + void SLAAutoSupports::uniformly_cover(const ExPolygon& island, Structure& structure, bool is_new_island, bool just_one) { //int num_of_points = std::max(1, (int)((island.area()*pow(SCALING_FACTOR, 2) * m_config.tear_pressure)/m_config.support_force)); - + const float density_horizontal = m_config.tear_pressure / m_config.support_force; + //FIXME why? + const float poisson_radius = 1.f / (5.f * density_horizontal); +// const float poisson_radius = 1.f / (15.f * density_horizontal); + const float samples_per_mm2 = 30.f / (float(M_PI) * poisson_radius * poisson_radius); - // We will cover the island another way. - // For now we'll just place the points randomly not too close to the others. //FIXME share the random generator. The random generator may be not so cheap to initialize, also we don't want the random generator to be restarted for each polygon. - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_real_distribution<> dis(0., 1.); + std::random_device rd; + std::mt19937 rng(rd()); + std::vector raw_samples = sample_expolygon_with_boundary(island, samples_per_mm2, 5.f / poisson_radius, rng); + std::vector poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius); - std::vector island_new_points; - const BoundingBox bb = get_extents(island); - const int refused_limit = (int)floor(30.f * ((float)bb.size()(0)*bb.size()(1) / (float)island.area()) + 0.5f); - int refused_points = 0; - //FIXME this is very inefficient (may be blind) for long narrow polygons (thin crescent, thin ring). Use some search structure: Triangulate the polygon first? - //FIXME use a low discrepancy sequence, Poisson sampling. - // Poisson sampling (adapt from 3D to 2D by triangulation): https://github.com/zewt/maya-implicit-skinning/blob/master/src/meshes/vcg_lib/utils_sampling.cpp - while (refused_points < refused_limit) { - Point out; - if (refused_points == 0 && island_new_points.empty()) // first iteration - out = island.contour.centroid(); - else - out = Point(bb.min(0) + bb.size()(0) * dis(gen), bb.min(1) + bb.size()(1) * dis(gen)); +#ifdef SLA_AUTOSUPPORTS_DEBUG + { + static int irun = 0; + Slic3r::SVG svg(debug_out_path("SLA_supports-uniformly_cover-%d.svg", irun ++), get_extents(island)); + svg.draw(island); + for (const Vec2f &pt : raw_samples) + svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "red"); + for (const Vec2f &pt : poisson_samples) + svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "blue"); + } +#endif /* NDEBUG */ - Vec3d unscaled_out = unscale(out(0), out(1), 0); - bool add_it = true; - - if (!island.contour.contains(out)) - add_it = false; - else - for (const Polygon& hole : island.holes) - if (hole.contains(out)) { - add_it = false; - break; - } - - if (add_it) { - for (const Vec3d& p : island_new_points) { - if ((p - unscaled_out).squaredNorm() < 1./(2.4*density_horizontal)) { - add_it = false; - break; - } - } - } - if (add_it) { - island_new_points.emplace_back(unscaled_out); - if (just_one) - break; - } - else - ++refused_points; - } - - for (const Vec3d& p : island_new_points) { - m_output.emplace_back(float(p(0)), float(p(1)), structure.height, 0.4f, is_new_island); + assert(! poisson_samples.empty()); + for (const Vec2f &pt : poisson_samples) { + m_output.emplace_back(float(pt(0)), float(pt(1)), structure.height, 0.4f, is_new_island); structure.supports_force += m_config.support_force; } } diff --git a/src/libslic3r/SLA/SLACommon.hpp b/src/libslic3r/SLA/SLACommon.hpp index dcf991aec6..ab535fb9ed 100644 --- a/src/libslic3r/SLA/SLACommon.hpp +++ b/src/libslic3r/SLA/SLACommon.hpp @@ -3,6 +3,7 @@ #include +// #define SLIC3R_SLA_NEEDS_WINDTREE namespace Slic3r { @@ -11,6 +12,7 @@ typedef Eigen::Matrix Vec3f; typedef Eigen::Matrix Vec3d; class TriangleMesh; + namespace sla { struct SupportPoint { @@ -115,11 +117,13 @@ public: int F_idx() const { return m_fidx; } }; +#ifdef SLIC3R_SLA_NEEDS_WINDTREE // The signed distance from a point to the mesh. Outputs the distance, // the index of the triangle and the closest point in mesh coordinate space. si_result signed_distance(const Vec3d& p) const; bool inside(const Vec3d& p) const; +#endif /* SLIC3R_SLA_NEEDS_WINDTREE */ }; diff --git a/src/libslic3r/SLA/SLASupportTree.cpp b/src/libslic3r/SLA/SLASupportTree.cpp index 87c967f8bb..3e36dfd851 100644 --- a/src/libslic3r/SLA/SLASupportTree.cpp +++ b/src/libslic3r/SLA/SLASupportTree.cpp @@ -552,14 +552,12 @@ enum { // For indexing Eigen vectors as v(X), v(Y), v(Z) instead of numbers PointSet to_point_set(const std::vector &v) { - PointSet ret(v.size(), 5); + PointSet ret(v.size(), 3); long i = 0; for(const SupportPoint& support_point : v) { ret.row(i)(0) = support_point.pos(0); ret.row(i)(1) = support_point.pos(1); ret.row(i)(2) = support_point.pos(2); - ret.row(i)(3) = support_point.head_front_radius; - ret.row(i)(4) = (float)support_point.is_new_island; ++i; } return ret; diff --git a/src/libslic3r/SLA/SLASupportTreeIGL.cpp b/src/libslic3r/SLA/SLASupportTreeIGL.cpp index d3af1eac89..6a3b71e7d6 100644 --- a/src/libslic3r/SLA/SLASupportTreeIGL.cpp +++ b/src/libslic3r/SLA/SLASupportTreeIGL.cpp @@ -95,7 +95,9 @@ size_t SpatIndex::size() const class EigenMesh3D::AABBImpl: public igl::AABB { public: +#ifdef SLIC3R_SLA_NEEDS_WINDTREE igl::WindingNumberAABB windtree; +#endif /* SLIC3R_SLA_NEEDS_WINDTREE */ }; EigenMesh3D::EigenMesh3D(const TriangleMesh& tmesh): m_aabb(new AABBImpl()) { @@ -136,7 +138,9 @@ EigenMesh3D::EigenMesh3D(const TriangleMesh& tmesh): m_aabb(new AABBImpl()) { // Build the AABB accelaration tree m_aabb->init(m_V, m_F); +#ifdef SLIC3R_SLA_NEEDS_WINDTREE m_aabb->windtree.set_mesh(m_V, m_F); +#endif /* SLIC3R_SLA_NEEDS_WINDTREE */ } EigenMesh3D::~EigenMesh3D() {} @@ -168,6 +172,7 @@ EigenMesh3D::query_ray_hit(const Vec3d &s, const Vec3d &dir) const return ret; } +#ifdef SLIC3R_SLA_NEEDS_WINDTREE EigenMesh3D::si_result EigenMesh3D::signed_distance(const Vec3d &p) const { double sign = 0; double sqdst = 0; int i = 0; Vec3d c; igl::signed_distance_winding_number(*m_aabb, m_V, m_F, m_aabb->windtree, @@ -179,6 +184,7 @@ EigenMesh3D::si_result EigenMesh3D::signed_distance(const Vec3d &p) const { bool EigenMesh3D::inside(const Vec3d &p) const { return m_aabb->windtree.inside(p); } +#endif /* SLIC3R_SLA_NEEDS_WINDTREE */ /* **************************************************************************** * Misc functions diff --git a/src/libslic3r/SLAPrint.cpp b/src/libslic3r/SLAPrint.cpp index e00437862f..4865150c21 100644 --- a/src/libslic3r/SLAPrint.cpp +++ b/src/libslic3r/SLAPrint.cpp @@ -25,7 +25,7 @@ using SupportTreePtr = std::unique_ptr; class SLAPrintObject::SupportData { public: sla::EigenMesh3D emesh; // index-triangle representation - sla::PointSet support_points; // all the support points (manual/auto) + std::vector support_points; // all the support points (manual/auto) SupportTreePtr support_tree_ptr; // the supports SlicedSupports support_slices; // sliced supports std::vector level_ids; @@ -548,18 +548,17 @@ void SLAPrint::process() // Now let's extract the result. const std::vector& points = auto_supports.output(); this->throw_if_canceled(); - po.m_supportdata->support_points = sla::to_point_set(points); + po.m_supportdata->support_points = points; BOOST_LOG_TRIVIAL(debug) << "Automatic support points: " - << po.m_supportdata->support_points.rows(); + << po.m_supportdata->support_points.size(); // Using RELOAD_SLA_SUPPORT_POINTS to tell the Plater to pass the update status to GLGizmoSlaSupports report_status(*this, -1, L("Generating support points"), SlicingStatus::RELOAD_SLA_SUPPORT_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()); + po.m_supportdata->support_points = po.transformed_support_points(); } }; @@ -597,7 +596,7 @@ void SLAPrint::process() ctl.cancelfn = [this]() { throw_if_canceled(); }; po.m_supportdata->support_tree_ptr.reset( - new SLASupportTree(po.m_supportdata->support_points, + new SLASupportTree(sla::to_point_set(po.m_supportdata->support_points), po.m_supportdata->emesh, scfg, ctl)); // Create the unified mesh @@ -608,7 +607,7 @@ void SLAPrint::process() po.m_supportdata->support_tree_ptr->merged_mesh(); BOOST_LOG_TRIVIAL(debug) << "Processed support point count " - << po.m_supportdata->support_points.rows(); + << po.m_supportdata->support_points.size(); // Check the mesh for later troubleshooting. if(po.support_mesh().empty()) @@ -1168,7 +1167,7 @@ const std::vector EMPTY_SLICES; const TriangleMesh EMPTY_MESH; } -const Eigen::MatrixXd& SLAPrintObject::get_support_points() const +const std::vector& SLAPrintObject::get_support_points() const { return m_supportdata->support_points; } diff --git a/src/libslic3r/SLAPrint.hpp b/src/libslic3r/SLAPrint.hpp index fc4e9a4d4d..9aaba6a479 100644 --- a/src/libslic3r/SLAPrint.hpp +++ b/src/libslic3r/SLAPrint.hpp @@ -91,7 +91,7 @@ public: const std::vector& get_support_slices() const; // This method returns the support points of this SLAPrintObject. - const Eigen::MatrixXd& get_support_points() const; + const std::vector& get_support_points() const; // An index record referencing the slices // (get_model_slices(), get_support_slices()) where the keys are the height diff --git a/src/slic3r/GUI/GLGizmo.cpp b/src/slic3r/GUI/GLGizmo.cpp index f9191b69b6..eecf1de2c2 100644 --- a/src/slic3r/GUI/GLGizmo.cpp +++ b/src/slic3r/GUI/GLGizmo.cpp @@ -1793,10 +1793,10 @@ void GLGizmoSlaSupports::set_sla_support_data(ModelObject* model_object, const G if (m_editing_mode_cache.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; itrafo().inverse().cast() * Vec3f(points(i,0), points(i,1), points(i,2)), - points(i, 3), points(i, 4)); + const std::vector& points = po->get_support_points(); + auto mat = po->trafo().inverse().cast(); + for (unsigned int i=0; i