1 //#include "igl/random_points_on_mesh.h"
2 //#include "igl/AABB.h"
3
4 #include <tbb/parallel_for.h>
5
6 #include "SupportPointGenerator.hpp"
7 #include "Concurrency.hpp"
8 #include "Model.hpp"
9 #include "ExPolygon.hpp"
10 #include "SVG.hpp"
11 #include "Point.hpp"
12 #include "ClipperUtils.hpp"
13 #include "Tesselate.hpp"
14 #include "ExPolygonCollection.hpp"
15 #include "libslic3r.h"
16
17 #include "libnest2d/backends/clipper/geometries.hpp"
18 #include "libnest2d/utils/rotcalipers.hpp"
19
20 #include <iostream>
21 #include <random>
22
23 namespace Slic3r {
24 namespace sla {
25
26 /*float SupportPointGenerator::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2)
27 {
28 n1.normalize();
29 n2.normalize();
30
31 Vec3d v = (p2-p1);
32 v.normalize();
33
34 float c1 = n1.dot(v);
35 float c2 = n2.dot(v);
36 float result = pow(p1(0)-p2(0), 2) + pow(p1(1)-p2(1), 2) + pow(p1(2)-p2(2), 2);
37 // Check for division by zero:
38 if(fabs(c1 - c2) > 0.0001)
39 result *= (asin(c1) - asin(c2)) / (c1 - c2);
40 return result;
41 }
42
43
44 float SupportPointGenerator::get_required_density(float angle) const
45 {
46 // calculation would be density_0 * cos(angle). To provide one more degree of freedom, we will scale the angle
47 // to get the user-set density for 45 deg. So it ends up as density_0 * cos(K * angle).
48 float K = 4.f * float(acos(m_config.density_at_45/m_config.density_at_horizontal) / M_PI);
49 return std::max(0.f, float(m_config.density_at_horizontal * cos(K*angle)));
50 }
51
52 float SupportPointGenerator::distance_limit(float angle) const
53 {
54 return 1./(2.4*get_required_density(angle));
55 }*/
56
SupportPointGenerator(const sla::IndexedMesh & emesh,const std::vector<ExPolygons> & slices,const std::vector<float> & heights,const Config & config,std::function<void (void)> throw_on_cancel,std::function<void (int)> statusfn)57 SupportPointGenerator::SupportPointGenerator(
58 const sla::IndexedMesh &emesh,
59 const std::vector<ExPolygons> &slices,
60 const std::vector<float> & heights,
61 const Config & config,
62 std::function<void(void)> throw_on_cancel,
63 std::function<void(int)> statusfn)
64 : SupportPointGenerator(emesh, config, throw_on_cancel, statusfn)
65 {
66 std::random_device rd;
67 m_rng.seed(rd());
68 execute(slices, heights);
69 }
70
SupportPointGenerator(const IndexedMesh & emesh,const SupportPointGenerator::Config & config,std::function<void ()> throw_on_cancel,std::function<void (int)> statusfn)71 SupportPointGenerator::SupportPointGenerator(
72 const IndexedMesh &emesh,
73 const SupportPointGenerator::Config &config,
74 std::function<void ()> throw_on_cancel,
75 std::function<void (int)> statusfn)
76 : m_config(config)
77 , m_emesh(emesh)
78 , m_throw_on_cancel(throw_on_cancel)
79 , m_statusfn(statusfn)
80 {
81 }
82
execute(const std::vector<ExPolygons> & slices,const std::vector<float> & heights)83 void SupportPointGenerator::execute(const std::vector<ExPolygons> &slices,
84 const std::vector<float> & heights)
85 {
86 process(slices, heights);
87 project_onto_mesh(m_output);
88 }
89
project_onto_mesh(std::vector<sla::SupportPoint> & points) const90 void SupportPointGenerator::project_onto_mesh(std::vector<sla::SupportPoint>& points) const
91 {
92 // The function makes sure that all the points are really exactly placed on the mesh.
93
94 // Use a reasonable granularity to account for the worker thread synchronization cost.
95 static constexpr size_t gransize = 64;
96
97 ccr_par::for_each(size_t(0), points.size(), [this, &points](size_t idx)
98 {
99 if ((idx % 16) == 0)
100 // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
101 m_throw_on_cancel();
102
103 Vec3f& p = points[idx].pos;
104 // Project the point upward and downward and choose the closer intersection with the mesh.
105 sla::IndexedMesh::hit_result hit_up = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., 1.));
106 sla::IndexedMesh::hit_result hit_down = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., -1.));
107
108 bool up = hit_up.is_hit();
109 bool down = hit_down.is_hit();
110
111 if (!up && !down)
112 return;
113
114 sla::IndexedMesh::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down;
115 p = p + (hit.distance() * hit.direction()).cast<float>();
116 }, gransize);
117 }
118
make_layers(const std::vector<ExPolygons> & slices,const std::vector<float> & heights,std::function<void (void)> throw_on_cancel)119 static std::vector<SupportPointGenerator::MyLayer> make_layers(
120 const std::vector<ExPolygons>& slices, const std::vector<float>& heights,
121 std::function<void(void)> throw_on_cancel)
122 {
123 assert(slices.size() == heights.size());
124
125 // Allocate empty layers.
126 std::vector<SupportPointGenerator::MyLayer> layers;
127 layers.reserve(slices.size());
128 for (size_t i = 0; i < slices.size(); ++ i)
129 layers.emplace_back(i, heights[i]);
130
131 // FIXME: calculate actual pixel area from printer config:
132 //const float pixel_area = pow(wxGetApp().preset_bundle->project_config.option<ConfigOptionFloat>("display_width") / wxGetApp().preset_bundle->project_config.option<ConfigOptionInt>("display_pixels_x"), 2.f); //
133 const float pixel_area = pow(0.047f, 2.f);
134
135 ccr_par::for_each(size_t(0), layers.size(),
136 [&layers, &slices, &heights, pixel_area, throw_on_cancel](size_t layer_id)
137 {
138 if ((layer_id % 8) == 0)
139 // Don't call the following function too often as it flushes
140 // CPU write caches due to synchronization primitves.
141 throw_on_cancel();
142
143 SupportPointGenerator::MyLayer &layer = layers[layer_id];
144 const ExPolygons & islands = slices[layer_id];
145 // FIXME WTF?
146 const float height = (layer_id > 2 ?
147 heights[layer_id - 3] :
148 heights[0] - (heights[1] - heights[0]));
149 layer.islands.reserve(islands.size());
150 for (const ExPolygon &island : islands) {
151 float area = float(island.area() * SCALING_FACTOR * SCALING_FACTOR);
152 if (area >= pixel_area)
153 // FIXME this is not a correct centroid of a polygon with holes.
154 layer.islands.emplace_back(layer, island, get_extents(island.contour),
155 unscaled<float>(island.contour.centroid()), area, height);
156 }
157 }, 32 /*gransize*/);
158
159 // Calculate overlap of successive layers. Link overlapping islands.
160 ccr_par::for_each(size_t(1), layers.size(),
161 [&layers, &heights, throw_on_cancel] (size_t layer_id)
162 {
163 if ((layer_id % 2) == 0)
164 // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
165 throw_on_cancel();
166 SupportPointGenerator::MyLayer &layer_above = layers[layer_id];
167 SupportPointGenerator::MyLayer &layer_below = layers[layer_id - 1];
168 //FIXME WTF?
169 const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]);
170 const float safe_angle = 35.f * (float(M_PI)/180.f); // smaller number - less supports
171 const float between_layers_offset = scaled<float>(layer_height * std::tan(safe_angle));
172 const float slope_angle = 75.f * (float(M_PI)/180.f); // smaller number - less supports
173 const float slope_offset = scaled<float>(layer_height * std::tan(slope_angle));
174 //FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands.
175 for (SupportPointGenerator::Structure &top : layer_above.islands) {
176 for (SupportPointGenerator::Structure &bottom : layer_below.islands) {
177 float overlap_area = top.overlap_area(bottom);
178 if (overlap_area > 0) {
179 top.islands_below.emplace_back(&bottom, overlap_area);
180 bottom.islands_above.emplace_back(&top, overlap_area);
181 }
182 }
183 if (! top.islands_below.empty()) {
184 Polygons top_polygons = to_polygons(*top.polygon);
185 Polygons bottom_polygons = top.polygons_below();
186 top.overhangs = diff_ex(top_polygons, bottom_polygons);
187 if (! top.overhangs.empty()) {
188
189 // Produce 2 bands around the island, a safe band for dangling overhangs
190 // and an unsafe band for sloped overhangs.
191 // These masks include the original island
192 auto dangl_mask = offset(bottom_polygons, between_layers_offset, ClipperLib::jtSquare);
193 auto overh_mask = offset(bottom_polygons, slope_offset, ClipperLib::jtSquare);
194
195 // Absolutely hopeless overhangs are those outside the unsafe band
196 top.overhangs = diff_ex(top_polygons, overh_mask);
197
198 // Now cut out the supported core from the safe band
199 // and cut the safe band from the unsafe band to get distinct
200 // zones.
201 overh_mask = diff(overh_mask, dangl_mask);
202 dangl_mask = diff(dangl_mask, bottom_polygons);
203
204 top.dangling_areas = intersection_ex(top_polygons, dangl_mask);
205 top.overhangs_slopes = intersection_ex(top_polygons, overh_mask);
206
207 top.overhangs_area = 0.f;
208 std::vector<std::pair<ExPolygon*, float>> expolys_with_areas;
209 for (ExPolygon &ex : top.overhangs) {
210 float area = float(ex.area());
211 expolys_with_areas.emplace_back(&ex, area);
212 top.overhangs_area += area;
213 }
214 std::sort(expolys_with_areas.begin(), expolys_with_areas.end(),
215 [](const std::pair<ExPolygon*, float> &p1, const std::pair<ExPolygon*, float> &p2)
216 { return p1.second > p2.second; });
217 ExPolygons overhangs_sorted;
218 for (auto &p : expolys_with_areas)
219 overhangs_sorted.emplace_back(std::move(*p.first));
220 top.overhangs = std::move(overhangs_sorted);
221 top.overhangs_area *= float(SCALING_FACTOR * SCALING_FACTOR);
222 }
223 }
224 }
225 }, 8 /* gransize */);
226
227 return layers;
228 }
229
process(const std::vector<ExPolygons> & slices,const std::vector<float> & heights)230 void SupportPointGenerator::process(const std::vector<ExPolygons>& slices, const std::vector<float>& heights)
231 {
232 #ifdef SLA_SUPPORTPOINTGEN_DEBUG
233 std::vector<std::pair<ExPolygon, coord_t>> islands;
234 #endif /* SLA_SUPPORTPOINTGEN_DEBUG */
235
236 std::vector<SupportPointGenerator::MyLayer> layers = make_layers(slices, heights, m_throw_on_cancel);
237
238 PointGrid3D point_grid;
239 point_grid.cell_size = Vec3f(10.f, 10.f, 10.f);
240
241 double increment = 100.0 / layers.size();
242 double status = 0;
243
244 for (unsigned int layer_id = 0; layer_id < layers.size(); ++ layer_id) {
245 SupportPointGenerator::MyLayer *layer_top = &layers[layer_id];
246 SupportPointGenerator::MyLayer *layer_bottom = (layer_id > 0) ? &layers[layer_id - 1] : nullptr;
247 std::vector<float> support_force_bottom;
248 if (layer_bottom != nullptr) {
249 support_force_bottom.assign(layer_bottom->islands.size(), 0.f);
250 for (size_t i = 0; i < layer_bottom->islands.size(); ++ i)
251 support_force_bottom[i] = layer_bottom->islands[i].supports_force_total();
252 }
253 for (Structure &top : layer_top->islands)
254 for (Structure::Link &bottom_link : top.islands_below) {
255 Structure &bottom = *bottom_link.island;
256 //float centroids_dist = (bottom.centroid - top.centroid).norm();
257 // Penalization resulting from centroid offset:
258 // bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, (1600.f * layer_height) * centroids_dist * centroids_dist / bottom.area));
259 float &support_force = support_force_bottom[&bottom - layer_bottom->islands.data()];
260 //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.
261 // One should rather work with the overlap area vs overhang area.
262 // support_force *= std::min(1.f, 1.f - std::min(1.f, 0.1f * centroids_dist * centroids_dist / bottom.area));
263 // Penalization resulting from increasing polygon area:
264 support_force *= std::min(1.f, 20.f * bottom.area / top.area);
265 }
266 // Let's assign proper support force to each of them:
267 if (layer_id > 0) {
268 for (Structure &below : layer_bottom->islands) {
269 float below_support_force = support_force_bottom[&below - layer_bottom->islands.data()];
270 float above_overlap_area = 0.f;
271 for (Structure::Link &above_link : below.islands_above)
272 above_overlap_area += above_link.overlap_area;
273 for (Structure::Link &above_link : below.islands_above)
274 above_link.island->supports_force_inherited += below_support_force * above_link.overlap_area / above_overlap_area;
275 }
276 }
277 // Now iterate over all polygons and append new points if needed.
278 for (Structure &s : layer_top->islands) {
279 // Penalization resulting from large diff from the last layer:
280 s.supports_force_inherited /= std::max(1.f, 0.17f * (s.overhangs_area) / s.area);
281
282 add_support_points(s, point_grid);
283 }
284
285 m_throw_on_cancel();
286
287 status += increment;
288 m_statusfn(int(std::round(status)));
289
290 #ifdef SLA_SUPPORTPOINTGEN_DEBUG
291 /*std::string layer_num_str = std::string((i<10 ? "0" : "")) + std::string((i<100 ? "0" : "")) + std::to_string(i);
292 output_expolygons(expolys_top, "top" + layer_num_str + ".svg");
293 output_expolygons(diff, "diff" + layer_num_str + ".svg");
294 if (!islands.empty())
295 output_expolygons(islands, "islands" + layer_num_str + ".svg");*/
296 #endif /* SLA_SUPPORTPOINTGEN_DEBUG */
297 }
298 }
299
add_support_points(SupportPointGenerator::Structure & s,SupportPointGenerator::PointGrid3D & grid3d)300 void SupportPointGenerator::add_support_points(SupportPointGenerator::Structure &s, SupportPointGenerator::PointGrid3D &grid3d)
301 {
302 // Select each type of surface (overrhang, dangling, slope), derive the support
303 // force deficit for it and call uniformly conver with the right params
304
305 float tp = m_config.tear_pressure();
306 float current = s.supports_force_total();
307 static constexpr float DANGL_DAMPING = .5f;
308 static constexpr float SLOPE_DAMPING = .1f;
309
310 if (s.islands_below.empty()) {
311 // completely new island - needs support no doubt
312 // deficit is full, there is nothing below that would hold this island
313 uniformly_cover({ *s.polygon }, s, s.area * tp, grid3d, IslandCoverageFlags(icfIsNew | icfWithBoundary) );
314 return;
315 }
316
317 if (! s.overhangs.empty()) {
318 uniformly_cover(s.overhangs, s, s.overhangs_area * tp, grid3d);
319 }
320
321 auto areafn = [](double sum, auto &p) { return sum + p.area() * SCALING_FACTOR * SCALING_FACTOR; };
322
323 current = s.supports_force_total();
324 if (! s.dangling_areas.empty()) {
325 // Let's see if there's anything that overlaps enough to need supports:
326 // What we now have in polygons needs support, regardless of what the forces are, so we can add them.
327
328 double a = std::accumulate(s.dangling_areas.begin(), s.dangling_areas.end(), 0., areafn);
329 uniformly_cover(s.dangling_areas, s, a * tp - a * current * s.area, grid3d, icfWithBoundary);
330 }
331
332 current = s.supports_force_total();
333 if (! s.overhangs_slopes.empty()) {
334 double a = std::accumulate(s.overhangs_slopes.begin(), s.overhangs_slopes.end(), 0., areafn);
335 uniformly_cover(s.overhangs_slopes, s, a * tp - a * current / s.area, grid3d, icfWithBoundary);
336 }
337 }
338
sample_expolygon(const ExPolygon & expoly,float samples_per_mm2,std::mt19937 & rng)339 std::vector<Vec2f> sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng)
340 {
341 // Triangulate the polygon with holes into triplets of 3D points.
342 std::vector<Vec2f> triangles = Slic3r::triangulate_expolygon_2f(expoly);
343
344 std::vector<Vec2f> out;
345 if (! triangles.empty())
346 {
347 // Calculate area of each triangle.
348 auto areas = reserve_vector<float>(triangles.size() / 3);
349 double aback = 0.;
350 for (size_t i = 0; i < triangles.size(); ) {
351 const Vec2f &a = triangles[i ++];
352 const Vec2f v1 = triangles[i ++] - a;
353 const Vec2f v2 = triangles[i ++] - a;
354
355 // Prefix sum of the areas.
356 areas.emplace_back(aback + 0.5f * std::abs(cross2(v1, v2)));
357 aback = areas.back();
358 }
359
360 size_t num_samples = size_t(ceil(areas.back() * samples_per_mm2));
361 std::uniform_real_distribution<> random_triangle(0., double(areas.back()));
362 std::uniform_real_distribution<> random_float(0., 1.);
363 for (size_t i = 0; i < num_samples; ++ i) {
364 double r = random_triangle(rng);
365 size_t idx_triangle = std::min<size_t>(std::upper_bound(areas.begin(), areas.end(), (float)r) - areas.begin(), areas.size() - 1) * 3;
366 // Select a random point on the triangle.
367 const Vec2f &a = triangles[idx_triangle ++];
368 const Vec2f &b = triangles[idx_triangle++];
369 const Vec2f &c = triangles[idx_triangle];
370 #if 1
371 // https://www.cs.princeton.edu/~funk/tog02.pdf
372 // page 814, formula 1.
373 double u = float(std::sqrt(random_float(rng)));
374 double v = float(random_float(rng));
375 out.emplace_back(a * (1.f - u) + b * (u * (1.f - v)) + c * (v * u));
376 #else
377 // Greg Turk, Graphics Gems
378 // https://devsplorer.wordpress.com/2019/08/07/find-a-random-point-on-a-plane-using-barycentric-coordinates-in-unity/
379 double u = float(random_float(rng));
380 double v = float(random_float(rng));
381 if (u + v >= 1.f) {
382 u = 1.f - u;
383 v = 1.f - v;
384 }
385 out.emplace_back(a + u * (b - a) + v * (c - a));
386 #endif
387 }
388 }
389 return out;
390 }
391
392
sample_expolygon(const ExPolygons & expolys,float samples_per_mm2,std::mt19937 & rng)393 std::vector<Vec2f> sample_expolygon(const ExPolygons &expolys, float samples_per_mm2, std::mt19937 &rng)
394 {
395 std::vector<Vec2f> out;
396 for (const ExPolygon &expoly : expolys)
397 append(out, sample_expolygon(expoly, samples_per_mm2, rng));
398
399 return out;
400 }
401
sample_expolygon_boundary(const ExPolygon & expoly,float samples_per_mm,std::vector<Vec2f> & out,std::mt19937 & rng)402 void sample_expolygon_boundary(const ExPolygon & expoly,
403 float samples_per_mm,
404 std::vector<Vec2f> &out,
405 std::mt19937 & rng)
406 {
407 double point_stepping_scaled = scale_(1.f) / samples_per_mm;
408 for (size_t i_contour = 0; i_contour <= expoly.holes.size(); ++ i_contour) {
409 const Polygon &contour = (i_contour == 0) ? expoly.contour :
410 expoly.holes[i_contour - 1];
411
412 const Points pts = contour.equally_spaced_points(point_stepping_scaled);
413 for (size_t i = 0; i < pts.size(); ++ i)
414 out.emplace_back(unscale<float>(pts[i].x()),
415 unscale<float>(pts[i].y()));
416 }
417 }
418
sample_expolygon_with_boundary(const ExPolygon & expoly,float samples_per_mm2,float samples_per_mm_boundary,std::mt19937 & rng)419 std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygon &expoly, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
420 {
421 std::vector<Vec2f> out = sample_expolygon(expoly, samples_per_mm2, rng);
422 sample_expolygon_boundary(expoly, samples_per_mm_boundary, out, rng);
423 return out;
424 }
425
sample_expolygon_with_boundary(const ExPolygons & expolys,float samples_per_mm2,float samples_per_mm_boundary,std::mt19937 & rng)426 std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygons &expolys, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
427 {
428 std::vector<Vec2f> out;
429 for (const ExPolygon &expoly : expolys)
430 append(out, sample_expolygon_with_boundary(expoly, samples_per_mm2, samples_per_mm_boundary, rng));
431 return out;
432 }
433
434 template<typename REFUSE_FUNCTION>
poisson_disk_from_samples(const std::vector<Vec2f> & raw_samples,float radius,REFUSE_FUNCTION refuse_function)435 static inline std::vector<Vec2f> poisson_disk_from_samples(const std::vector<Vec2f> &raw_samples, float radius, REFUSE_FUNCTION refuse_function)
436 {
437 Vec2f corner_min(std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
438 for (const Vec2f &pt : raw_samples) {
439 corner_min.x() = std::min(corner_min.x(), pt.x());
440 corner_min.y() = std::min(corner_min.y(), pt.y());
441 }
442
443 // Assign the raw samples to grid cells, sort the grid cells lexicographically.
444 struct RawSample
445 {
446 Vec2f coord;
447 Vec2i cell_id;
448 RawSample(const Vec2f &crd = {}, const Vec2i &id = {}): coord{crd}, cell_id{id} {}
449 };
450
451 auto raw_samples_sorted = reserve_vector<RawSample>(raw_samples.size());
452 for (const Vec2f &pt : raw_samples)
453 raw_samples_sorted.emplace_back(pt, ((pt - corner_min) / radius).cast<int>());
454
455 std::sort(raw_samples_sorted.begin(), raw_samples_sorted.end(), [](const RawSample &lhs, const RawSample &rhs)
456 { 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()); });
457
458 struct PoissonDiskGridEntry {
459 // Resulting output sample points for this cell:
460 enum {
461 max_positions = 4
462 };
463 Vec2f poisson_samples[max_positions];
464 int num_poisson_samples = 0;
465
466 // Index into raw_samples:
467 int first_sample_idx;
468 int sample_cnt;
469 };
470
471 struct CellIDHash {
472 std::size_t operator()(const Vec2i &cell_id) const {
473 return std::hash<int>()(cell_id.x()) ^ std::hash<int>()(cell_id.y() * 593);
474 }
475 };
476
477 // Map from cell IDs to hash_data. Each hash_data points to the range in raw_samples corresponding to that cell.
478 // (We could just store the samples in hash_data. This implementation is an artifact of the reference paper, which
479 // is optimizing for GPU acceleration that we haven't implemented currently.)
480 typedef std::unordered_map<Vec2i, PoissonDiskGridEntry, CellIDHash> Cells;
481 Cells cells;
482 {
483 typename Cells::iterator last_cell_id_it;
484 Vec2i last_cell_id(-1, -1);
485 for (size_t i = 0; i < raw_samples_sorted.size(); ++ i) {
486 const RawSample &sample = raw_samples_sorted[i];
487 if (sample.cell_id == last_cell_id) {
488 // This sample is in the same cell as the previous, so just increase the count. Cells are
489 // always contiguous, since we've sorted raw_samples_sorted by cell ID.
490 ++ last_cell_id_it->second.sample_cnt;
491 } else {
492 // This is a new cell.
493 PoissonDiskGridEntry data;
494 data.first_sample_idx = int(i);
495 data.sample_cnt = 1;
496 auto result = cells.insert({sample.cell_id, data});
497 last_cell_id = sample.cell_id;
498 last_cell_id_it = result.first;
499 }
500 }
501 }
502
503 const int max_trials = 5;
504 const float radius_squared = radius * radius;
505 for (int trial = 0; trial < max_trials; ++ trial) {
506 // Create sample points for each entry in cells.
507 for (auto &it : cells) {
508 const Vec2i &cell_id = it.first;
509 PoissonDiskGridEntry &cell_data = it.second;
510 // 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.
511 int next_sample_idx = cell_data.first_sample_idx + trial;
512 if (trial >= cell_data.sample_cnt)
513 // There are no more points to try for this cell.
514 continue;
515 const RawSample &candidate = raw_samples_sorted[next_sample_idx];
516 // See if this point conflicts with any other points in this cell, or with any points in
517 // neighboring cells. Note that it's possible to have more than one point in the same cell.
518 bool conflict = refuse_function(candidate.coord);
519 for (int i = -1; i < 2 && ! conflict; ++ i) {
520 for (int j = -1; j < 2; ++ j) {
521 const auto &it_neighbor = cells.find(cell_id + Vec2i(i, j));
522 if (it_neighbor != cells.end()) {
523 const PoissonDiskGridEntry &neighbor = it_neighbor->second;
524 for (int i_sample = 0; i_sample < neighbor.num_poisson_samples; ++ i_sample)
525 if ((neighbor.poisson_samples[i_sample] - candidate.coord).squaredNorm() < radius_squared) {
526 conflict = true;
527 break;
528 }
529 }
530 }
531 }
532 if (! conflict) {
533 // Store the new sample.
534 assert(cell_data.num_poisson_samples < cell_data.max_positions);
535 if (cell_data.num_poisson_samples < cell_data.max_positions)
536 cell_data.poisson_samples[cell_data.num_poisson_samples ++] = candidate.coord;
537 }
538 }
539 }
540
541 // Copy the results to the output.
542 std::vector<Vec2f> out;
543 for (const auto& it : cells)
544 for (int i = 0; i < it.second.num_poisson_samples; ++ i)
545 out.emplace_back(it.second.poisson_samples[i]);
546 return out;
547 }
548
549
uniformly_cover(const ExPolygons & islands,Structure & structure,float deficit,PointGrid3D & grid3d,IslandCoverageFlags flags)550 void SupportPointGenerator::uniformly_cover(const ExPolygons& islands, Structure& structure, float deficit, PointGrid3D &grid3d, IslandCoverageFlags flags)
551 {
552 //int num_of_points = std::max(1, (int)((island.area()*pow(SCALING_FACTOR, 2) * m_config.tear_pressure)/m_config.support_force));
553
554 float support_force_deficit = deficit;
555 // auto bb = get_extents(islands);
556
557 if (flags & icfIsNew) {
558 auto chull_ex = ExPolygonCollection{islands}.convex_hull();
559 auto chull = Slic3rMultiPoint_to_ClipperPath(chull_ex);
560 auto rotbox = libnest2d::minAreaBoundingBox(chull);
561 Vec2d bbdim = {unscaled(rotbox.width()), unscaled(rotbox.height())};
562
563 if (bbdim.x() > bbdim.y()) std::swap(bbdim.x(), bbdim.y());
564 double aspectr = bbdim.y() / bbdim.x();
565
566 support_force_deficit *= (1 + aspectr / 2.);
567 }
568
569 if (support_force_deficit < 0)
570 return;
571
572 // Number of newly added points.
573 const size_t poisson_samples_target = size_t(ceil(support_force_deficit / m_config.support_force()));
574
575 const float density_horizontal = m_config.tear_pressure() / m_config.support_force();
576 //FIXME why?
577 float poisson_radius = std::max(m_config.minimal_distance, 1.f / (5.f * density_horizontal));
578 // const float poisson_radius = 1.f / (15.f * density_horizontal);
579 const float samples_per_mm2 = 30.f / (float(M_PI) * poisson_radius * poisson_radius);
580 // Minimum distance between samples, in 3D space.
581 // float min_spacing = poisson_radius / 3.f;
582 float min_spacing = poisson_radius;
583
584 //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.
585
586 std::vector<Vec2f> raw_samples =
587 flags & icfWithBoundary ?
588 sample_expolygon_with_boundary(islands, samples_per_mm2,
589 5.f / poisson_radius, m_rng) :
590 sample_expolygon(islands, samples_per_mm2, m_rng);
591
592 std::vector<Vec2f> poisson_samples;
593 for (size_t iter = 0; iter < 4; ++ iter) {
594 poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius,
595 [&structure, &grid3d, min_spacing](const Vec2f &pos) {
596 return grid3d.collides_with(pos, structure.layer->print_z, min_spacing);
597 });
598 if (poisson_samples.size() >= poisson_samples_target || m_config.minimal_distance > poisson_radius-EPSILON)
599 break;
600 float coeff = 0.5f;
601 if (poisson_samples.size() * 2 > poisson_samples_target)
602 coeff = float(poisson_samples.size()) / float(poisson_samples_target);
603 poisson_radius = std::max(m_config.minimal_distance, poisson_radius * coeff);
604 min_spacing = std::max(m_config.minimal_distance, min_spacing * coeff);
605 }
606
607 #ifdef SLA_SUPPORTPOINTGEN_DEBUG
608 {
609 static int irun = 0;
610 Slic3r::SVG svg(debug_out_path("SLA_supports-uniformly_cover-%d.svg", irun ++), get_extents(islands));
611 for (const ExPolygon &island : islands)
612 svg.draw(island);
613 for (const Vec2f &pt : raw_samples)
614 svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "red");
615 for (const Vec2f &pt : poisson_samples)
616 svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "blue");
617 }
618 #endif /* NDEBUG */
619
620 // assert(! poisson_samples.empty());
621 if (poisson_samples_target < poisson_samples.size()) {
622 std::shuffle(poisson_samples.begin(), poisson_samples.end(), m_rng);
623 poisson_samples.erase(poisson_samples.begin() + poisson_samples_target, poisson_samples.end());
624 }
625 for (const Vec2f &pt : poisson_samples) {
626 m_output.emplace_back(float(pt(0)), float(pt(1)), structure.zlevel, m_config.head_diameter/2.f, flags & icfIsNew);
627 structure.supports_force_this_layer += m_config.support_force();
628 grid3d.insert(pt, &structure);
629 }
630 }
631
632
remove_bottom_points(std::vector<SupportPoint> & pts,float lvl)633 void remove_bottom_points(std::vector<SupportPoint> &pts, float lvl)
634 {
635 // get iterator to the reorganized vector end
636 auto endit = std::remove_if(pts.begin(), pts.end(), [lvl]
637 (const sla::SupportPoint &sp) {
638 return sp.pos.z() <= lvl;
639 });
640
641 // erase all elements after the new end
642 pts.erase(endit, pts.end());
643 }
644
645 #ifdef SLA_SUPPORTPOINTGEN_DEBUG
output_structures(const std::vector<Structure> & structures)646 void SupportPointGenerator::output_structures(const std::vector<Structure>& structures)
647 {
648 for (unsigned int i=0 ; i<structures.size(); ++i) {
649 std::stringstream ss;
650 ss << structures[i].unique_id.count() << "_" << std::setw(10) << std::setfill('0') << 1000 + (int)structures[i].height/1000 << ".png";
651 output_expolygons(std::vector<ExPolygon>{*structures[i].polygon}, ss.str());
652 }
653 }
654
output_expolygons(const ExPolygons & expolys,const std::string & filename)655 void SupportPointGenerator::output_expolygons(const ExPolygons& expolys, const std::string &filename)
656 {
657 BoundingBox bb(Point(-30000000, -30000000), Point(30000000, 30000000));
658 Slic3r::SVG svg_cummulative(filename, bb);
659 for (size_t i = 0; i < expolys.size(); ++ i) {
660 /*Slic3r::SVG svg("single"+std::to_string(i)+".svg", bb);
661 svg.draw(expolys[i]);
662 svg.draw_outline(expolys[i].contour, "black", scale_(0.05));
663 svg.draw_outline(expolys[i].holes, "blue", scale_(0.05));
664 svg.Close();*/
665
666 svg_cummulative.draw(expolys[i]);
667 svg_cummulative.draw_outline(expolys[i].contour, "black", scale_(0.05));
668 svg_cummulative.draw_outline(expolys[i].holes, "blue", scale_(0.05));
669 }
670 }
671 #endif
672
673 } // namespace sla
674 } // namespace Slic3r
675