From 672cf5d45f913e3eb74722bd9fcfe113f8ba5d96 Mon Sep 17 00:00:00 2001 From: bubnikv Date: Mon, 18 Feb 2019 11:46:06 +0100 Subject: [PATCH] Parallelization of large part of the SLA support point calculation. New 3D grid to check, whether two support points are not too close. --- src/libslic3r/SLA/SLAAutoSupports.cpp | 204 +++++++++++++++++--------- src/libslic3r/SLA/SLAAutoSupports.hpp | 100 +++++++++++-- 2 files changed, 220 insertions(+), 84 deletions(-) diff --git a/src/libslic3r/SLA/SLAAutoSupports.cpp b/src/libslic3r/SLA/SLAAutoSupports.cpp index 38add28eec..7321387af4 100644 --- a/src/libslic3r/SLA/SLAAutoSupports.cpp +++ b/src/libslic3r/SLA/SLAAutoSupports.cpp @@ -93,99 +93,152 @@ void SLAAutoSupports::project_onto_mesh(std::vector& points) }); } +static std::vector make_layers( + const std::vector& slices, const std::vector& heights, + std::function throw_on_cancel) +{ + 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]); + + // 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); // + const float pixel_area = pow(0.047f, 2.f); + + // Use a reasonable granularity to account for the worker thread synchronization cost. + tbb::parallel_for(tbb::blocked_range(0, layers.size(), 32), + [&layers, &slices, &heights, pixel_area, throw_on_cancel](const tbb::blocked_range& range) { + for (size_t layer_id = range.begin(); layer_id < range.end(); ++ 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(); + SLAAutoSupports::MyLayer &layer = layers[layer_id]; + const ExPolygons &islands = slices[layer_id]; + //FIXME WTF? + const float height = (layer_id>2 ? heights[layer_id-3] : heights[0]-(heights[1]-heights[0])); + layer.islands.reserve(islands.size()); + for (const ExPolygon &island : islands) { + float area = float(island.area() * SCALING_FACTOR * SCALING_FACTOR); + if (area >= pixel_area) + //FIXME this is not a correct centroid of a polygon with holes. + layer.islands.emplace_back(layer, island, get_extents(island.contour), Slic3r::unscale(island.contour.centroid()).cast(), area, height); + } + } + }); + + // Calculate overlap of successive layers. Link overlapping islands. + tbb::parallel_for(tbb::blocked_range(1, layers.size(), 8), + [&layers, &heights, throw_on_cancel](const tbb::blocked_range& range) { + for (size_t layer_id = range.begin(); layer_id < range.end(); ++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(); + SLAAutoSupports::MyLayer &layer_above = layers[layer_id]; + SLAAutoSupports::MyLayer &layer_below = layers[layer_id - 1]; + //FIXME WTF? + const float height = (layer_id>2 ? heights[layer_id-3] : heights[0]-(heights[1]-heights[0])); + const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]); + const float safe_angle = 5.f * (float(M_PI)/180.f); // smaller number - less supports + const float between_layers_offset = float(scale_(layer_height / std::tan(safe_angle))); + //FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands. + for (SLAAutoSupports::Structure &top : layer_above.islands) { + for (SLAAutoSupports::Structure &bottom : layer_below.islands) + if (top.overlaps(bottom)) { + top.islands_below.emplace_back(&bottom); + bottom.islands_above.emplace_back(&top); + } + if (! top.islands_below.empty()) { + Polygons top_polygons = to_polygons(*top.polygon); + Polygons bottom_polygons = top.polygons_below(); + top.overhangs = diff_ex(top_polygons, bottom_polygons); + if (! top.overhangs.empty()) { + top.overhangs_area = 0.f; + 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); + top.dangling_areas = diff_ex(top_polygons, offset(bottom_polygons, between_layers_offset)); + } + } + } + } + }); + + return layers; +} + void SLAAutoSupports::process(const std::vector& slices, const std::vector& heights) { +#ifdef SLA_AUTOSUPPORTS_DEBUG std::vector> islands; - std::vector structures_old; - std::vector structures_new; +#endif /* SLA_AUTOSUPPORTS_DEBUG */ - for (unsigned int i = 0; i layers = make_layers(slices, heights, m_throw_on_cancel); - //FIXME WTF? - const float height = (i>2 ? heights[i-3] : heights[0]-(heights[1]-heights[0])); - const float layer_height = (i!=0 ? heights[i]-heights[i-1] : heights[0]); - - const float safe_angle = 5.f * (float(M_PI)/180.f); // smaller number - less supports - const float between_layers_offset = float(scale_(layer_height / std::tan(safe_angle))); + PointGrid3D point_grid; + point_grid.cell_size = Vec3f(10.f, 10.f, 10.f); - // 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); // - const float pixel_area = pow(0.047f, 2.f); - - // Check all ExPolygons on this slice and check whether they are new or belonging to something below. - for (const ExPolygon& polygon : expolys_top) { - float area = float(polygon.area() * SCALING_FACTOR * SCALING_FACTOR); - if (area < pixel_area) - continue; - //FIXME this is not a correct centroid of a polygon with holes. - structures_new.emplace_back(polygon, get_extents(polygon.contour), Slic3r::unscale(polygon.contour.centroid()).cast(), area, height); - 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); - float centroids_dist = (bottom.centroid - top.centroid).norm(); - // Penalization resulting from centroid offset: + for (unsigned int layer_id = 0; layer_id < layers.size(); ++ layer_id) { + SLAAutoSupports::MyLayer &layer_top = layers[layer_id]; + for (Structure &top : layer_top.islands) + for (Structure *bottom : top.islands_below) { + 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, 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, 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); + } // Let's assign proper support force to each of them: - for (const Structure& below : structures_old) { - std::vector above_list; - float above_area = 0.f; - for (Structure& new_str : structures_new) - for (const Structure* below1 : new_str.structures_below) - if (&below == below1) { - above_list.push_back(&new_str); - above_area += above_list.back()->area; - } - for (Structure* above : above_list) - above->supports_force += below.supports_force * above->area / above_area; + if (layer_id > 0) { + for (Structure &below : layers[layer_id - 1].islands) { + float above_area = 0.f; + for (Structure *above : below.islands_above) + above_area += above->area; + for (Structure *above : below.islands_above) + above->supports_force += below.supports_force * above->area / above_area; + } } - // Now iterate over all polygons and append new points if needed. - for (Structure& s : structures_new) { - if (s.structures_below.empty()) // completely new island - needs support no doubt - uniformly_cover(*s.polygon, s, true); + for (Structure &s : layer_top.islands) { + if (s.islands_below.empty()) // completely new island - needs support no doubt + uniformly_cover(*s.polygon, s, point_grid, true); 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 : s.dangling_areas) //FIXME is it an island point or not? Vojtech thinks it is. - uniformly_cover(p, s); + uniformly_cover(p, s, point_grid); } // We should also check if current support is enough given the polygon area. - 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()); + for (Structure& s : layer_top.islands) { // 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 * (s.overhangs_area) / s.area); if (s.area * m_config.tear_pressure > s.supports_force) { //FIXME Don't calculate area inside the compare function! //FIXME Cover until the force deficit is covered. Cover multiple areas, sort by decreasing area. - ExPolygons::iterator largest_it = std::max_element(e.begin(), e.end(), [](const ExPolygon& a, const ExPolygon& b) { return a.area() < b.area(); }); - if (!e.empty()) + if (! s.overhangs.empty()) //FIXME add the support force deficit as a parameter, only cover until the defficiency is covered. - uniformly_cover(*largest_it, s); + uniformly_cover(s.overhangs.front(), s, point_grid); } } - // All is done. Prepare to advance to the next layer. - structures_old = std::move(structures_new); - structures_new.clear(); - m_throw_on_cancel(); #ifdef SLA_AUTOSUPPORTS_DEBUG @@ -251,7 +304,8 @@ std::vector sample_expolygon_with_boundary(const ExPolygon &expoly, float return out; } -std::vector poisson_disk_from_samples(const std::vector &raw_samples, float radius) +template +static inline std::vector poisson_disk_from_samples(const std::vector &raw_samples, float radius, REFUSE_FUNCTION refuse_function) { Vec2f corner_min(FLT_MAX, FLT_MAX); for (const Vec2f &pt : raw_samples) { @@ -334,7 +388,7 @@ std::vector poisson_disk_from_samples(const std::vector &raw_sampl 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; + 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)); @@ -365,7 +419,7 @@ std::vector poisson_disk_from_samples(const std::vector &raw_sampl return out; } -void SLAAutoSupports::uniformly_cover(const ExPolygon& island, Structure& structure, bool is_new_island, bool just_one) +void SLAAutoSupports::uniformly_cover(const ExPolygon& island, Structure& structure, PointGrid3D &grid3d, 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)); @@ -374,12 +428,17 @@ void SLAAutoSupports::uniformly_cover(const ExPolygon& island, Structure& struct 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); + // Minimum distance between samples, in 3D space. + const float min_spacing = poisson_radius / 3.f; //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 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 poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius, + [&structure, &grid3d, min_spacing](const Vec2f &pos) { + return grid3d.collides_with(pos, &structure, min_spacing); + }); #ifdef SLA_AUTOSUPPORTS_DEBUG { @@ -393,10 +452,11 @@ void SLAAutoSupports::uniformly_cover(const ExPolygon& island, Structure& struct } #endif /* NDEBUG */ - assert(! poisson_samples.empty()); +// assert(! poisson_samples.empty()); for (const Vec2f &pt : poisson_samples) { m_output.emplace_back(float(pt(0)), float(pt(1)), structure.height, 0.2f, is_new_island); structure.supports_force += m_config.support_force; + grid3d.insert(pt, &structure); } } diff --git a/src/libslic3r/SLA/SLAAutoSupports.hpp b/src/libslic3r/SLA/SLAAutoSupports.hpp index fb5aceb4a2..17f94bbd6e 100644 --- a/src/libslic3r/SLA/SLAAutoSupports.hpp +++ b/src/libslic3r/SLA/SLAAutoSupports.hpp @@ -5,6 +5,8 @@ #include #include +#include + // #define SLA_AUTOSUPPORTS_DEBUG namespace Slic3r { @@ -24,22 +26,20 @@ public: const std::vector& heights, const Config& config, std::function throw_on_cancel); const std::vector& output() { return m_output; } -private: - std::vector m_output; - - SLAAutoSupports::Config m_config; + struct MyLayer; struct Structure { - Structure(const ExPolygon& poly, const BoundingBox &bbox, const Vec2f ¢roid, float area, float h) : polygon(&poly), bbox(bbox), centroid(centroid), area(area), height(h) + Structure(MyLayer &layer, const ExPolygon& poly, const BoundingBox &bbox, const Vec2f ¢roid, float area, float h) : + layer(&layer), polygon(&poly), bbox(bbox), centroid(centroid), area(area), height(h) #ifdef SLA_AUTOSUPPORTS_DEBUG , unique_id(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch())) #endif /* SLA_AUTOSUPPORTS_DEBUG */ {} + MyLayer *layer; const ExPolygon* polygon = nullptr; const BoundingBox bbox; const Vec2f centroid = Vec2f::Zero(); const float area = 0.f; - std::vector structures_below; float height = 0; // How well is this ExPolygon held to the print base? // Positive number, the higher the better. @@ -48,20 +48,26 @@ private: std::chrono::milliseconds unique_id; #endif /* SLA_AUTOSUPPORTS_DEBUG */ + boost::container::small_vector islands_above; + boost::container::small_vector islands_below; + ExPolygons dangling_areas; + ExPolygons overhangs; + float overhangs_area; + bool overlaps(const Structure &rhs) const { return this->bbox.overlap(rhs.bbox) && (this->polygon->overlaps(*rhs.polygon) || rhs.polygon->overlaps(*this->polygon)); } float area_below() const { float area = 0.f; - for (const Structure *below : this->structures_below) + for (const Structure *below : this->islands_below) area += below->area; return area; } Polygons polygons_below() const { size_t cnt = 0; - for (const Structure *below : this->structures_below) + for (const Structure *below : this->islands_below) cnt += 1 + below->polygon->holes.size(); Polygons out; out.reserve(cnt); - for (const Structure *below : this->structures_below) { + for (const Structure *below : this->islands_below) { out.emplace_back(below->polygon->contour); append(out, below->polygon->holes); } @@ -69,17 +75,87 @@ private: } ExPolygons expolygons_below() const { ExPolygons out; - out.reserve(this->structures_below.size()); - for (const Structure *below : this->structures_below) + out.reserve(this->islands_below.size()); + for (const Structure *below : this->islands_below) out.emplace_back(*below->polygon); return out; } }; + struct MyLayer { + MyLayer(const size_t layer_id, coordf_t print_z) : layer_id(layer_id), print_z(print_z) {} + size_t layer_id; + coordf_t print_z; + std::vector islands; + }; + + struct RichSupportPoint { + Vec3f position; + Structure *island; + }; + + struct PointGrid3D { + struct GridHash { + std::size_t operator()(const Vec3i &cell_id) { + 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, Structure *island, float radius) { + Vec3f pos3d(pos.x(), pos.y(), float(island->layer->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; + } + }; + +private: + std::vector m_output; + + SLAAutoSupports::Config m_config; + float m_supports_force_total = 0.f; void process(const std::vector& slices, const std::vector& heights); - void uniformly_cover(const ExPolygon& island, Structure& structure, bool is_new_island = false, bool just_one = false); + void uniformly_cover(const ExPolygon& island, Structure& structure, PointGrid3D &grid3d, bool is_new_island = false, bool just_one = false); void project_onto_mesh(std::vector& points) const; #ifdef SLA_AUTOSUPPORTS_DEBUG