From f1f0fbc9ad19567347b2883a1c33e5bfe6cec0f6 Mon Sep 17 00:00:00 2001 From: Filip Sykala - NTB T15p Date: Wed, 21 Aug 2024 15:31:12 +0200 Subject: [PATCH] Add new implementation of samling --- src/libslic3r/CMakeLists.txt | 2 + .../SLA/SupportPointGeneratorNew.cpp | 677 ++++++++++++++++++ .../SLA/SupportPointGeneratorNew.hpp | 281 ++++++++ 3 files changed, 960 insertions(+) create mode 100644 src/libslic3r/SLA/SupportPointGeneratorNew.cpp create mode 100644 src/libslic3r/SLA/SupportPointGeneratorNew.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index d8c8d1e069..46d81c2367 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -444,6 +444,8 @@ set(SLIC3R_SOURCES SLA/SupportPoint.hpp SLA/SupportPointGenerator.hpp SLA/SupportPointGenerator.cpp + SLA/SupportPointGeneratorNew.hpp + SLA/SupportPointGeneratorNew.cpp SLA/Clustering.hpp SLA/Clustering.cpp SLA/ReprojectPointsOnMesh.hpp diff --git a/src/libslic3r/SLA/SupportPointGeneratorNew.cpp b/src/libslic3r/SLA/SupportPointGeneratorNew.cpp new file mode 100644 index 0000000000..473f4819c3 --- /dev/null +++ b/src/libslic3r/SLA/SupportPointGeneratorNew.cpp @@ -0,0 +1,677 @@ +///|/ Copyright (c) Prusa Research 2019 - 2023 Tomáš Mészáros @tamasmeszaros, Vojtěch Bubník @bubnikv, Lukáš Matěna @lukasmatena +///|/ +///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher +///|/ +#include +#include +#include +#include +#include +#include + +#include "SupportPointGenerator.hpp" +#include "libslic3r/Execution/ExecutionTBB.hpp" +#include "libslic3r/Geometry/ConvexHull.hpp" +#include "libslic3r/Model.hpp" +#include "libslic3r/ExPolygon.hpp" +#include "libslic3r/Point.hpp" +#include "libslic3r/ClipperUtils.hpp" +#include "libslic3r/Tesselate.hpp" +#include "libslic3r/MinAreaBoundingBox.hpp" +#include "libslic3r/libslic3r.h" +#include "libslic3r/AABBMesh.hpp" +#include "libslic3r/Execution/Execution.hpp" +#include "libslic3r/SLA/SupportPoint.hpp" +#include "libslic3r/BoundingBox.hpp" + +namespace Slic3r { +namespace sla { + +/*float SupportPointGenerator::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2) +{ + n1.normalize(); + n2.normalize(); + + Vec3d v = (p2-p1); + v.normalize(); + + float c1 = n1.dot(v); + float c2 = n2.dot(v); + float result = pow(p1(0)-p2(0), 2) + pow(p1(1)-p2(1), 2) + pow(p1(2)-p2(2), 2); + // Check for division by zero: + if(fabs(c1 - c2) > 0.0001) + result *= (asin(c1) - asin(c2)) / (c1 - c2); + return result; +} + + +float SupportPointGenerator::get_required_density(float angle) const +{ + // calculation would be density_0 * cos(angle). To provide one more degree of freedom, we will scale the angle + // to get the user-set density for 45 deg. So it ends up as density_0 * cos(K * angle). + float K = 4.f * float(acos(m_config.density_at_45/m_config.density_at_horizontal) / M_PI); + return std::max(0.f, float(m_config.density_at_horizontal * cos(K*angle))); +} + +float SupportPointGenerator::distance_limit(float angle) const +{ + return 1./(2.4*get_required_density(angle)); +}*/ + +SupportPointGenerator::SupportPointGenerator( + const AABBMesh &emesh, + const std::vector &slices, + const std::vector & heights, + const Config & config, + std::function throw_on_cancel, + std::function statusfn) + : m_config(config) + , m_emesh(emesh) + , m_throw_on_cancel(throw_on_cancel) + , m_statusfn(statusfn) +{ + std::random_device rd; + m_rng.seed(rd()); + execute(slices, heights); +} + +void SupportPointGenerator::execute(const std::vector &slices, + const std::vector & heights) +{ + process(slices, heights); + project_onto_mesh(m_output); +} + +void SupportPointGenerator::project_onto_mesh(std::vector& points) const +{ + // The function makes sure that all the points are really exactly placed on the mesh. + + // Use a reasonable granularity to account for the worker thread synchronization cost. + static constexpr size_t gransize = 64; + + execution::for_each(ex_tbb, size_t(0), points.size(), [this, &points](size_t idx) + { + if ((idx % 16) == 0) + // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves. + m_throw_on_cancel(); + + Vec3f& p = points[idx].pos; + // Project the point upward and downward and choose the closer intersection with the mesh. + AABBMesh::hit_result hit_up = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., 1.)); + AABBMesh::hit_result hit_down = m_emesh.query_ray_hit(p.cast(), Vec3d(0., 0., -1.)); + + bool up = hit_up.is_hit(); + bool down = hit_down.is_hit(); + + if (!up && !down) + return; + + AABBMesh::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down; + p = p + (hit.distance() * hit.direction()).cast(); + }, gransize); +} + +static std::vector make_layers( + const std::vector& slices, const std::vector& heights, + std::function throw_on_cancel, const sla::SupportPointGenerator::Config& config) +{ + assert(slices.size() == heights.size()); + + // Allocate empty layers. + std::vector layers; + layers.reserve(slices.size()); + for (size_t i = 0; i < slices.size(); ++ i) + layers.emplace_back(i, heights[i]); + + execution::for_each(ex_tbb, size_t(0), layers.size(), + [&layers, &slices, min_area = config.minimal_island_area, throw_on_cancel](size_t layer_id) + { + if ((layer_id % 8) == 0) + // Don't call the following function too often as it flushes + // CPU write caches due to synchronization primitves. + throw_on_cancel(); + + SupportPointGenerator::MyLayer &layer = layers[layer_id]; + const ExPolygons & islands = slices[layer_id]; + layer.islands.reserve(islands.size()); + for (const ExPolygon &island : islands) { + float area = float(island.area() * SCALING_FACTOR * SCALING_FACTOR); + + // Skip too small islands (non-printable) + if (area < min_area) + continue; + + layer.islands.emplace_back(layer, island, get_extents(island.contour), area); + } + }, 32 /*gransize*/); + + // Calculate overlap of successive layers. Link overlapping islands. + execution::for_each(ex_tbb, size_t(1), layers.size(), + [&layers, &heights, throw_on_cancel] (size_t layer_id) + { + if ((layer_id % 2) == 0) + // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves. + throw_on_cancel(); + SupportPointGenerator::MyLayer &layer_above = layers[layer_id]; + SupportPointGenerator::MyLayer &layer_below = layers[layer_id - 1]; + //FIXME WTF? + const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]); + const float safe_angle = 35.f * (float(M_PI)/180.f); // smaller number - less supports + const float between_layers_offset = scaled(layer_height * std::tan(safe_angle)); + const float slope_angle = 75.f * (float(M_PI)/180.f); // smaller number - less supports + const float slope_offset = scaled(layer_height * std::tan(slope_angle)); + //FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands. + for (SupportPointGenerator::Structure &top : layer_above.islands) { + for (SupportPointGenerator::Structure &bottom : layer_below.islands) { + float overlap_area = top.overlap_area(bottom); + if (overlap_area > 0) { + top.islands_below.emplace_back(&bottom, overlap_area); + bottom.islands_above.emplace_back(&top, overlap_area); + } + } + if (! top.islands_below.empty()) { + // Why only polygons?? (some sort of speed up?) + Polygons bottom_polygons = top.polygons_below(); + top.overhangs = diff_ex(*top.polygon, bottom_polygons); + if (! top.overhangs.empty()) { + + // Produce 2 bands around the island, a safe band for dangling overhangs + // and an unsafe band for sloped overhangs. + // These masks include the original island + auto dangl_mask = expand(bottom_polygons, between_layers_offset, ClipperLib::jtSquare); + auto overh_mask = expand(bottom_polygons, slope_offset, ClipperLib::jtSquare); + + // Absolutely hopeless overhangs are those outside the unsafe band + top.overhangs = diff_ex(*top.polygon, overh_mask); + + // Now cut out the supported core from the safe band + // and cut the safe band from the unsafe band to get distinct + // zones. + overh_mask = diff(overh_mask, dangl_mask); + dangl_mask = diff(dangl_mask, bottom_polygons); + + top.dangling_areas = intersection_ex(*top.polygon, dangl_mask); + top.overhangs_slopes = intersection_ex(*top.polygon, overh_mask); + + top.overhangs_area = 0.f; + + // Sort overhangs by area + std::vector> expolys_with_areas; + for (ExPolygon &ex : top.overhangs) { + float area = float(ex.area()); + expolys_with_areas.emplace_back(&ex, area); + top.overhangs_area += area; + } + std::sort(expolys_with_areas.begin(), expolys_with_areas.end(), + [](const std::pair &p1, const std::pair &p2) + { return p1.second > p2.second; }); + ExPolygons overhangs_sorted; + for (auto &p : expolys_with_areas) + overhangs_sorted.emplace_back(std::move(*p.first)); + top.overhangs = std::move(overhangs_sorted); + top.overhangs_area *= float(SCALING_FACTOR * SCALING_FACTOR); + } + } + } + }, 8 /* gransize */); + + return layers; +} + +void SupportPointGenerator::process(const std::vector& slices, const std::vector& heights) +{ +#ifdef SLA_SUPPORTPOINTGEN_DEBUG + std::vector> islands; +#endif /* SLA_SUPPORTPOINTGEN_DEBUG */ + + std::vector layers = make_layers(slices, heights, m_throw_on_cancel, m_config); + + PointGrid3D point_grid; + point_grid.cell_size = Vec3f(10.f, 10.f, 10.f); + + double increment = 100.0 / layers.size(); + double status = 0; + + std::vector support_force_bottom; + for (unsigned int layer_id = 0; layer_id < layers.size(); ++ layer_id) { + bool is_first_layer = layer_id == 0; + SupportPointGenerator::MyLayer *layer_top = &layers[layer_id]; + SupportPointGenerator::MyLayer *layer_bottom = (!is_first_layer) ? &layers[layer_id - 1] : nullptr; + if (!is_first_layer) { + support_force_bottom.resize(layer_bottom->islands.size()); + for (size_t i = 0; i < layer_bottom->islands.size(); ++i) + support_force_bottom[i] = layer_bottom->islands[i].supports_force_total(); + } + for (const Structure &top : layer_top->islands) + for (const Structure::Link &bottom_link : top.islands_below) { + const Structure &bottom = *bottom_link.island; + //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)); + size_t bottom_island_index = &bottom - layer_bottom->islands.data(); + float &support_force = support_force_bottom[bottom_island_index]; + //FIXME this condition does not reflect a bifurcation into a one large island and one tiny island well, it incorrectly resets the support force to zero. +// One should rather work with the overlap area vs overhang area. +// support_force *= std::min(1.f, 1.f - std::min(1.f, 0.1f * centroids_dist * centroids_dist / bottom.area)); + // Penalization resulting from increasing polygon area: + support_force *= std::min(1.f, 20.f * bottom.area / top.area); + } + // Let's assign proper support force to each of them: + if (!is_first_layer) { + for (Structure &below : layer_bottom->islands) { + float below_support_force = support_force_bottom[&below - layer_bottom->islands.data()]; + float above_overlap_area = 0.f; + for (Structure::Link &above_link : below.islands_above) + above_overlap_area += above_link.overlap_area; + for (Structure::Link &above_link : below.islands_above) + above_link.island->supports_force_inherited += below_support_force * above_link.overlap_area / above_overlap_area; + } + } + // Now iterate over all polygons and append new points if needed. + for (Structure &s : layer_top->islands) { + // Penalization resulting from large diff from the last layer: + s.supports_force_inherited /= std::max(1.f, 0.17f * s.overhangs_area / s.area); + + add_support_points(s, point_grid); + } + + m_throw_on_cancel(); + + status += increment; + m_statusfn(int(std::round(status))); + } +} + +void SupportPointGenerator::add_support_points(SupportPointGenerator::Structure &s, SupportPointGenerator::PointGrid3D &grid3d) +{ + // Select each type of surface (overrhang, dangling, slope), derive the support + // force deficit for it and call uniformly conver with the right params + if (s.islands_below.empty()) { + // completely new island - needs support no doubt + // deficit is full, there is nothing below that would hold this island + uniformly_cover({ *s.polygon }, s, s.area, grid3d, IslandCoverageFlags(icfIsNew | icfWithBoundary) ); + return; + } + + if (! s.overhangs.empty()) { + uniformly_cover(s.overhangs, s, s.overhangs_area, grid3d); + } + + auto areafn = [](double sum, auto &p) { return sum + p.area() * SCALING_FACTOR * SCALING_FACTOR; }; + + float current = s.supports_force_total(); + if (! s.dangling_areas.empty()) { + // 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. + + // Before Tamas changes + // a * tp - current * .09f *std::sqrt(1. - a / s.area) + + // just befor + // a * ( 1 - current * s.area); + + double sum_of_dangling_area = std::accumulate(s.dangling_areas.begin(), s.dangling_areas.end(), 0., areafn); + double dangling_ratio = sum_of_dangling_area / s.area; + float deficit = current * .09f * dangling_ratio; + //uniformly_cover(s.dangling_areas, s, deficit, grid3d, icfWithBoundary); + } + + current = s.supports_force_total(); + if (! s.overhangs_slopes.empty()) { + double sum_of_overhang_area = std::accumulate(s.overhangs_slopes.begin(), s.overhangs_slopes.end(), 0., areafn); + double overhang_ratio = sum_of_overhang_area / s.area; + float deficit = current * .0015f * overhang_ratio; + //uniformly_cover(s.overhangs_slopes, s, deficit, grid3d, icfWithBoundary); + } +} + +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. + auto areas = reserve_vector(triangles.size() / 3); + double aback = 0.; + 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; + + // Prefix sum of the areas. + areas.emplace_back(aback + 0.5f * std::abs(cross2(v1, v2))); + aback = areas.back(); + } + + 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. + const Vec2f &a = triangles[idx_triangle ++]; + const Vec2f &b = triangles[idx_triangle++]; + const Vec2f &c = triangles[idx_triangle]; +#if 1 + // https://www.cs.princeton.edu/~funk/tog02.pdf + // page 814, formula 1. + double u = float(std::sqrt(random_float(rng))); + double v = float(random_float(rng)); + out.emplace_back(a * (1.f - u) + b * (u * (1.f - v)) + c * (v * u)); +#else + // Greg Turk, Graphics Gems + // https://devsplorer.wordpress.com/2019/08/07/find-a-random-point-on-a-plane-using-barycentric-coordinates-in-unity/ + double u = float(random_float(rng)); + double v = float(random_float(rng)); + if (u + v >= 1.f) { + u = 1.f - u; + v = 1.f - v; + } + out.emplace_back(a + u * (b - a) + v * (c - a)); +#endif + } + } + return out; +} + + +std::vector sample_expolygon(const ExPolygons &expolys, float samples_per_mm2, std::mt19937 &rng) +{ + std::vector out; + for (const ExPolygon &expoly : expolys) + append(out, sample_expolygon(expoly, samples_per_mm2, rng)); + + return out; +} + +void sample_expolygon_boundary(const ExPolygon & expoly, + float samples_per_mm, + std::vector &out, + std::mt19937 & /*rng*/) +{ + double point_stepping_scaled = scale_(1.f) / samples_per_mm; + 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())); + } +} + +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); + sample_expolygon_boundary(expoly, samples_per_mm_boundary, out, rng); + return out; +} + +std::vector sample_expolygon_with_boundary(const ExPolygons &expolys, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng) +{ + std::vector out; + for (const ExPolygon &expoly : expolys) + append(out, sample_expolygon_with_boundary(expoly, samples_per_mm2, samples_per_mm_boundary, rng)); + return out; +} + +template +static inline std::vector poisson_disk_from_samples(const std::vector &raw_samples, float radius, REFUSE_FUNCTION refuse_function) +{ + Vec2f corner_min(std::numeric_limits::max(), std::numeric_limits::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; + RawSample(const Vec2f &crd = {}, const Vec2i &id = {}): coord{crd}, cell_id{id} {} + }; + + auto raw_samples_sorted = reserve_vector(raw_samples.size()); + for (const Vec2f &pt : raw_samples) + raw_samples_sorted.emplace_back(pt, ((pt - corner_min) / radius).cast()); + + 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) const { + 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; + Cells cells; + { + typename Cells::iterator last_cell_id_it; + Vec2i last_cell_id(-1, -1); + for (size_t 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 = int(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 = refuse_function(candidate.coord); + 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 SupportPointGenerator::uniformly_cover(const ExPolygons& islands, Structure& structure, float deficit, PointGrid3D &grid3d, IslandCoverageFlags flags) +{ + //int num_of_points = std::max(1, (int)((island.area()*pow(SCALING_FACTOR, 2) * m_config.tear_pressure)/m_config.support_force)); + + float support_force_deficit = deficit; +// auto bb = get_extents(islands); + + if (flags & icfIsNew) { + auto chull = Geometry::convex_hull(islands); + auto rotbox = MinAreaBoundigBox{chull, MinAreaBoundigBox::pcConvex}; + Vec2d bbdim = {unscaled(rotbox.width()), unscaled(rotbox.height())}; + + if (bbdim.x() > bbdim.y()) std::swap(bbdim.x(), bbdim.y()); + double aspectr = bbdim.y() / bbdim.x(); + + support_force_deficit *= (1 + aspectr / 2.); + } + + if (support_force_deficit < 0) + return; + + // Number of newly added points. + const size_t poisson_samples_target = size_t(ceil(support_force_deficit / m_config.support_force())); + + const float density_horizontal = 1. / m_config.support_force(); + //FIXME why? + float poisson_radius = std::max(m_config.minimal_distance, 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); + // Minimum distance between samples, in 3D space. +// float min_spacing = poisson_radius / 3.f; + float min_spacing = poisson_radius; + + //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::vector raw_samples = + flags & icfWithBoundary ? + sample_expolygon_with_boundary(islands, samples_per_mm2, + 5.f / poisson_radius, m_rng) : + sample_expolygon(islands, samples_per_mm2, m_rng); + + std::vector poisson_samples; + for (size_t iter = 0; iter < 4; ++ iter) { + poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius, + [&structure, &grid3d, min_spacing](const Vec2f &pos) { + return grid3d.collides_with(pos, structure.layer->print_z, min_spacing); + }); + if (poisson_samples.size() >= poisson_samples_target || m_config.minimal_distance > poisson_radius-EPSILON) + break; + float coeff = 0.5f; + if (poisson_samples.size() * 2 > poisson_samples_target) + coeff = float(poisson_samples.size()) / float(poisson_samples_target); + poisson_radius = std::max(m_config.minimal_distance, poisson_radius * coeff); + min_spacing = std::max(m_config.minimal_distance, min_spacing * coeff); + } + +#ifdef SLA_SUPPORTPOINTGEN_DEBUG + { + static int irun = 0; + Slic3r::SVG svg(debug_out_path("SLA_supports-uniformly_cover-%d.svg", irun ++), get_extents(islands)); + for (const ExPolygon &island : islands) + 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 */ + +// assert(! poisson_samples.empty()); + if (poisson_samples_target < poisson_samples.size()) { + std::shuffle(poisson_samples.begin(), poisson_samples.end(), m_rng); + poisson_samples.erase(poisson_samples.begin() + poisson_samples_target, poisson_samples.end()); + } + for (const Vec2f &pt : poisson_samples) { + m_output.emplace_back( + float(pt(0)), float(pt(1)), structure.layer->print_z/*structure.zlevel*/, m_config.head_diameter / 2.f, + flags & icfIsNew + ); + structure.supports_force_this_layer += m_config.support_force(); + grid3d.insert(pt, &structure); + } +} + + +void remove_bottom_points(std::vector &pts, float lvl) +{ + // get iterator to the reorganized vector end + auto endit = std::remove_if(pts.begin(), pts.end(), [lvl] + (const sla::SupportPoint &sp) { + return sp.pos.z() <= lvl; + }); + + // erase all elements after the new end + pts.erase(endit, pts.end()); +} + +#ifdef SLA_SUPPORTPOINTGEN_DEBUG +void SupportPointGenerator::output_structures(const std::vector& structures) +{ + for (unsigned int i=0 ; i{*structures[i].polygon}, ss.str()); + } +} + +void SupportPointGenerator::output_expolygons(const ExPolygons& expolys, const std::string &filename) +{ + BoundingBox bb(Point(-30000000, -30000000), Point(30000000, 30000000)); + Slic3r::SVG svg_cummulative(filename, bb); + for (size_t i = 0; i < expolys.size(); ++ i) { + /*Slic3r::SVG svg("single"+std::to_string(i)+".svg", bb); + svg.draw(expolys[i]); + svg.draw_outline(expolys[i].contour, "black", scale_(0.05)); + svg.draw_outline(expolys[i].holes, "blue", scale_(0.05)); + svg.Close();*/ + + svg_cummulative.draw(expolys[i]); + svg_cummulative.draw_outline(expolys[i].contour, "black", scale_(0.05)); + svg_cummulative.draw_outline(expolys[i].holes, "blue", scale_(0.05)); + } +} +#endif + +SupportPoints transformed_support_points(const ModelObject &mo, + const Transform3d &trafo) +{ + auto spts = mo.sla_support_points; + Transform3f tr = trafo.cast(); + for (sla::SupportPoint& suppt : spts) { + suppt.pos = tr * suppt.pos; + } + + return spts; +} + +} // namespace sla +} // namespace Slic3r diff --git a/src/libslic3r/SLA/SupportPointGeneratorNew.hpp b/src/libslic3r/SLA/SupportPointGeneratorNew.hpp new file mode 100644 index 0000000000..466dab2503 --- /dev/null +++ b/src/libslic3r/SLA/SupportPointGeneratorNew.hpp @@ -0,0 +1,281 @@ +///|/ Copyright (c) Prusa Research 2024 Filip Sykala @Jony01 +///|/ +///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher +///|/ +#ifndef SLA_SUPPORTPOINTGENERATOR_HPP +#define SLA_SUPPORTPOINTGENERATOR_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "libslic3r/ExPolygon.hpp" +#include "libslic3r/Polygon.hpp" +#include "libslic3r/libslic3r.h" + +namespace Slic3r { +class AABBMesh; +} // namespace Slic3r + +// #define SLA_SUPPORTPOINTGEN_DEBUG + +namespace Slic3r { namespace sla { + +/// +/// Configuration for automatic support placement +/// +struct SupportPointGeneratorConfig{ + + /// + /// 0 mean only one support point for each island + /// lower than one mean less amount of support points + /// 1 mean fine tuned sampling + /// more than one mean bigger amout of support points + /// + float density_relative{1.f}; + + /// + /// Minimal distance between support points + /// Should correspond with density relative + /// + float minimal_distance{1.f}; + + /// + /// size of supported point - range in future + /// + float head_diameter{0.4f}; + + // FIXME: calculate actual pixel area from printer config: + // const float pixel_area = + // pow(wxGetApp().preset_bundle->project_config.option("display_width") / + // wxGetApp().preset_bundle->project_config.option("display_pixels_x"), 2.f); // + // Minimal island Area to print - TODO: Should be modifiable from UI + float minimal_island_area = pow(0.047f, 2.f); // [in mm^2] pixel_area +}; + + +class SupportPointGenerator { +public: + + + SupportPointGenerator(const AABBMesh& emesh, const std::vector& slices, const std::vector& heights, + std::function throw_on_cancel, std::function statusfn); + + const std::vector& output() const { return m_output; } + std::vector& output() { return m_output; } + + struct MyLayer; + + // Keep data for one area(ExPlygon) on the layer(on slice Expolygons) + struct Structure { + Structure(MyLayer &layer, const ExPolygon& poly, const BoundingBox &bbox, float area) : + layer(&layer), polygon(&poly), bbox(bbox), area(area) +#ifdef SLA_SUPPORTPOINTGEN_DEBUG + , unique_id(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch())) +#endif /* SLA_SUPPORTPOINTGEN_DEBUG */ + {} + // Parent layer - with all ExPolygons in layer + layer_height + MyLayer *layer; + // Source ExPolygon + const ExPolygon* polygon = nullptr; + // Cache bounding box of polygon + const BoundingBox bbox; + // area of polygon [in mm^2] without holes + const float area = 0.f; + + // How well is this ExPolygon held to the print base? + // Positive number, the higher the better. + float supports_force_this_layer = 0.f; + float supports_force_inherited = 0.f; + float supports_force_total() const { return this->supports_force_this_layer + this->supports_force_inherited; } +#ifdef SLA_SUPPORTPOINTGEN_DEBUG + std::chrono::milliseconds unique_id; +#endif /* SLA_SUPPORTPOINTGEN_DEBUG */ + + struct Link { + Link(Structure *island, float overlap_area) : island(island), overlap_area(overlap_area) {} + Structure *island; + float overlap_area; + }; + +#ifdef NDEBUG + // In release mode, use the optimized container. + boost::container::small_vector islands_above; + boost::container::small_vector islands_below; +#else + // In debug mode, use the standard vector, which is well handled by debugger visualizer. + std::vector islands_above; + std::vector islands_below; +#endif + // Overhangs, that are dangling considerably. + ExPolygons dangling_areas; + // Complete overhangs. + ExPolygons overhangs; + // Overhangs, where the surface must slope. + ExPolygons overhangs_slopes; + // Sum of all overhang areas from structure + float overhangs_area = 0.f; // [in mm^2] + + bool overlaps(const Structure &rhs) const { + return this->bbox.overlap(rhs.bbox) && this->polygon->overlaps(*rhs.polygon); + } + float overlap_area(const Structure &rhs) const { + double out = 0.; + if (this->bbox.overlap(rhs.bbox)) { + Polygons polys = intersection(*this->polygon, *rhs.polygon); + for (const Polygon &poly : polys) + out += poly.area(); + } + return float(out); + } + float area_below() const { + float area = 0.f; + for (const Link &below : this->islands_below) + area += below.island->area; + return area; + } + Polygons polygons_below() const { + size_t cnt = 0; + for (const Link &below : this->islands_below) + cnt += 1 + below.island->polygon->holes.size(); + Polygons out; + out.reserve(cnt); + for (const Link &below : this->islands_below) { + out.emplace_back(below.island->polygon->contour); + append(out, below.island->polygon->holes); + } + return out; + } + ExPolygons expolygons_below() const { + ExPolygons out; + out.reserve(this->islands_below.size()); + for (const Link &below : this->islands_below) + out.emplace_back(*below.island->polygon); + return out; + } + // Positive deficit of the supports. If negative, this area is well supported. If positive, more supports need to be added. + float support_force_deficit(const float tear_pressure) const { return this->area * tear_pressure - this->supports_force_total(); } + }; + + struct MyLayer { + MyLayer(const size_t layer_id, coordf_t print_z) : layer_id(layer_id), print_z(print_z) {} + // index into heights + slices + size_t layer_id; + // Absolute distance from Zero - copy value from heights + coordf_t print_z; // [in mm] + std::vector islands; + }; + + struct RichSupportPoint { + Vec3f position; + Structure *island; + }; + + struct PointGrid3D { + struct GridHash { + std::size_t operator()(const Vec3i &cell_id) const { + return std::hash()(cell_id.x()) ^ std::hash()(cell_id.y() * 593) ^ std::hash()(cell_id.z() * 7919); + } + }; + typedef std::unordered_multimap Grid; + + Vec3f cell_size; + Grid grid; + + Vec3i cell_id(const Vec3f &pos) { + return Vec3i(int(floor(pos.x() / cell_size.x())), + int(floor(pos.y() / cell_size.y())), + int(floor(pos.z() / cell_size.z()))); + } + + void insert(const Vec2f &pos, Structure *island) { + RichSupportPoint pt; + pt.position = Vec3f(pos.x(), pos.y(), float(island->layer->print_z)); + pt.island = island; + grid.emplace(cell_id(pt.position), pt); + } + + bool collides_with(const Vec2f &pos, float print_z, float radius) { + Vec3f pos3d(pos.x(), pos.y(), print_z); + Vec3i cell = cell_id(pos3d); + std::pair it_pair = grid.equal_range(cell); + if (collides_with(pos3d, radius, it_pair.first, it_pair.second)) + return true; + for (int i = -1; i < 2; ++ i) + for (int j = -1; j < 2; ++ j) + for (int k = -1; k < 1; ++ k) { + if (i == 0 && j == 0 && k == 0) + continue; + it_pair = grid.equal_range(cell + Vec3i(i, j, k)); + if (collides_with(pos3d, radius, it_pair.first, it_pair.second)) + return true; + } + return false; + } + + private: + bool collides_with(const Vec3f &pos, float radius, Grid::const_iterator it_begin, Grid::const_iterator it_end) { + for (Grid::const_iterator it = it_begin; it != it_end; ++ it) { + float dist2 = (it->second.position - pos).squaredNorm(); + if (dist2 < radius * radius) + return true; + } + return false; + } + }; + + void execute(const std::vector &slices, + const std::vector & heights); + + void seed(std::mt19937::result_type s) { m_rng.seed(s); } +private: + std::vector m_output; + + // Configuration + SupportPointGenerator::Config m_config; + + void process(const std::vector& slices, const std::vector& heights); + +public: + enum IslandCoverageFlags : uint8_t { icfNone = 0x0, icfIsNew = 0x1, icfWithBoundary = 0x2 }; + +private: + + void uniformly_cover(const ExPolygons& islands, Structure& structure, float deficit, PointGrid3D &grid3d, IslandCoverageFlags flags = icfNone); + + void add_support_points(Structure& structure, PointGrid3D &grid3d); + + void project_onto_mesh(std::vector& points) const; + +#ifdef SLA_SUPPORTPOINTGEN_DEBUG + static void output_expolygons(const ExPolygons& expolys, const std::string &filename); + static void output_structures(const std::vector &structures); +#endif // SLA_SUPPORTPOINTGEN_DEBUG + + const AABBMesh& m_emesh; + std::function m_throw_on_cancel; + std::function m_statusfn; + + std::mt19937 m_rng; +}; + +void remove_bottom_points(std::vector &pts, float lvl); + +std::vector sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng); +void sample_expolygon_boundary(const ExPolygon &expoly, float samples_per_mm, std::vector &out, std::mt19937 &rng); + +}} // namespace Slic3r::sla + +#endif // SUPPORTPOINTGENERATOR_HPP