1 #include "../ClipperUtils.hpp"
2 #include "../ExPolygon.hpp"
3 #include "../Surface.hpp"
4 #include "../Geometry.hpp"
5 #include "../Layer.hpp"
6 #include "../Print.hpp"
7 #include "../ShortestPath.hpp"
8
9 #include "FillAdaptive.hpp"
10
11 // for indexed_triangle_set
12 #include <admesh/stl.h>
13
14 #include <cstdlib>
15 #include <cmath>
16 #include <algorithm>
17 #include <numeric>
18
19 // Boost pool: Don't use mutexes to synchronize memory allocation.
20 #define BOOST_POOL_NO_MT
21 #include <boost/pool/object_pool.hpp>
22
23 #include <boost/geometry.hpp>
24 #include <boost/geometry/geometries/point.hpp>
25 #include <boost/geometry/geometries/segment.hpp>
26 #include <boost/geometry/index/rtree.hpp>
27
28
29 namespace Slic3r {
30 namespace FillAdaptive {
31
32 // Derived from https://github.com/juj/MathGeoLib/blob/master/src/Geometry/Triangle.cpp
33 // The AABB-Triangle test implementation is based on the pseudo-code in
34 // Christer Ericson's Real-Time Collision Detection, pp. 169-172. It is
35 // practically a standard SAT test.
36 //
37 // Original MathGeoLib benchmark:
38 // Best: 17.282 nsecs / 46.496 ticks, Avg: 17.804 nsecs, Worst: 18.434 nsecs
39 //
40 //FIXME Vojtech: The MathGeoLib contains a vectorized implementation.
41 template<typename Vector>
triangle_AABB_intersects(const Vector & a,const Vector & b,const Vector & c,const BoundingBoxBase<Vector> & aabb)42 bool triangle_AABB_intersects(const Vector &a, const Vector &b, const Vector &c, const BoundingBoxBase<Vector> &aabb)
43 {
44 using Scalar = typename Vector::Scalar;
45
46 Vector tMin = a.cwiseMin(b.cwiseMin(c));
47 Vector tMax = a.cwiseMax(b.cwiseMax(c));
48
49 if (tMin.x() >= aabb.max.x() || tMax.x() <= aabb.min.x()
50 || tMin.y() >= aabb.max.y() || tMax.y() <= aabb.min.y()
51 || tMin.z() >= aabb.max.z() || tMax.z() <= aabb.min.z())
52 return false;
53
54 Vector center = (aabb.min + aabb.max) * 0.5f;
55 Vector h = aabb.max - center;
56
57 const Vector t[3] { b-a, c-a, c-b };
58
59 Vector ac = a - center;
60
61 Vector n = t[0].cross(t[1]);
62 Scalar s = n.dot(ac);
63 Scalar r = std::abs(h.dot(n.cwiseAbs()));
64 if (abs(s) >= r)
65 return false;
66
67 const Vector at[3] = { t[0].cwiseAbs(), t[1].cwiseAbs(), t[2].cwiseAbs() };
68
69 Vector bc = b - center;
70 Vector cc = c - center;
71
72 // SAT test all cross-axes.
73 // The following is a fully unrolled loop of this code, stored here for reference:
74 /*
75 Scalar d1, d2, a1, a2;
76 const Vector e[3] = { DIR_VEC(1, 0, 0), DIR_VEC(0, 1, 0), DIR_VEC(0, 0, 1) };
77 for(int i = 0; i < 3; ++i)
78 for(int j = 0; j < 3; ++j)
79 {
80 Vector axis = Cross(e[i], t[j]);
81 ProjectToAxis(axis, d1, d2);
82 aabb.ProjectToAxis(axis, a1, a2);
83 if (d2 <= a1 || d1 >= a2) return false;
84 }
85 */
86
87 // eX <cross> t[0]
88 Scalar d1 = t[0].y() * ac.z() - t[0].z() * ac.y();
89 Scalar d2 = t[0].y() * cc.z() - t[0].z() * cc.y();
90 Scalar tc = (d1 + d2) * 0.5f;
91 r = std::abs(h.y() * at[0].z() + h.z() * at[0].y());
92 if (r + std::abs(tc - d1) < std::abs(tc))
93 return false;
94
95 // eX <cross> t[1]
96 d1 = t[1].y() * ac.z() - t[1].z() * ac.y();
97 d2 = t[1].y() * bc.z() - t[1].z() * bc.y();
98 tc = (d1 + d2) * 0.5f;
99 r = std::abs(h.y() * at[1].z() + h.z() * at[1].y());
100 if (r + std::abs(tc - d1) < std::abs(tc))
101 return false;
102
103 // eX <cross> t[2]
104 d1 = t[2].y() * ac.z() - t[2].z() * ac.y();
105 d2 = t[2].y() * bc.z() - t[2].z() * bc.y();
106 tc = (d1 + d2) * 0.5f;
107 r = std::abs(h.y() * at[2].z() + h.z() * at[2].y());
108 if (r + std::abs(tc - d1) < std::abs(tc))
109 return false;
110
111 // eY <cross> t[0]
112 d1 = t[0].z() * ac.x() - t[0].x() * ac.z();
113 d2 = t[0].z() * cc.x() - t[0].x() * cc.z();
114 tc = (d1 + d2) * 0.5f;
115 r = std::abs(h.x() * at[0].z() + h.z() * at[0].x());
116 if (r + std::abs(tc - d1) < std::abs(tc))
117 return false;
118
119 // eY <cross> t[1]
120 d1 = t[1].z() * ac.x() - t[1].x() * ac.z();
121 d2 = t[1].z() * bc.x() - t[1].x() * bc.z();
122 tc = (d1 + d2) * 0.5f;
123 r = std::abs(h.x() * at[1].z() + h.z() * at[1].x());
124 if (r + std::abs(tc - d1) < std::abs(tc))
125 return false;
126
127 // eY <cross> t[2]
128 d1 = t[2].z() * ac.x() - t[2].x() * ac.z();
129 d2 = t[2].z() * bc.x() - t[2].x() * bc.z();
130 tc = (d1 + d2) * 0.5f;
131 r = std::abs(h.x() * at[2].z() + h.z() * at[2].x());
132 if (r + std::abs(tc - d1) < std::abs(tc))
133 return false;
134
135 // eZ <cross> t[0]
136 d1 = t[0].x() * ac.y() - t[0].y() * ac.x();
137 d2 = t[0].x() * cc.y() - t[0].y() * cc.x();
138 tc = (d1 + d2) * 0.5f;
139 r = std::abs(h.y() * at[0].x() + h.x() * at[0].y());
140 if (r + std::abs(tc - d1) < std::abs(tc))
141 return false;
142
143 // eZ <cross> t[1]
144 d1 = t[1].x() * ac.y() - t[1].y() * ac.x();
145 d2 = t[1].x() * bc.y() - t[1].y() * bc.x();
146 tc = (d1 + d2) * 0.5f;
147 r = std::abs(h.y() * at[1].x() + h.x() * at[1].y());
148 if (r + std::abs(tc - d1) < std::abs(tc))
149 return false;
150
151 // eZ <cross> t[2]
152 d1 = t[2].x() * ac.y() - t[2].y() * ac.x();
153 d2 = t[2].x() * bc.y() - t[2].y() * bc.x();
154 tc = (d1 + d2) * 0.5f;
155 r = std::abs(h.y() * at[2].x() + h.x() * at[2].y());
156 if (r + std::abs(tc - d1) < std::abs(tc))
157 return false;
158
159 // No separating axis exists, the AABB and triangle intersect.
160 return true;
161 }
162
dist2_to_triangle(const Vec3d & a,const Vec3d & b,const Vec3d & c,const Vec3d & p)163 static double dist2_to_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p)
164 {
165 double out = std::numeric_limits<double>::max();
166 const Vec3d v1 = b - a;
167 auto l1 = v1.squaredNorm();
168 const Vec3d v2 = c - b;
169 auto l2 = v2.squaredNorm();
170 const Vec3d v3 = a - c;
171 auto l3 = v3.squaredNorm();
172
173 // Is the triangle valid?
174 if (l1 > 0. && l2 > 0. && l3 > 0.)
175 {
176 // 1) Project point into the plane of the triangle.
177 const Vec3d n = v1.cross(v2);
178 double d = (p - a).dot(n);
179 const Vec3d foot_pt = p - n * d / n.squaredNorm();
180
181 // 2) Maximum projection of n.
182 int proj_axis;
183 n.array().cwiseAbs().maxCoeff(&proj_axis);
184
185 // 3) Test whether the foot_pt is inside the triangle.
186 {
187 auto inside_triangle = [](const Vec2d& v1, const Vec2d& v2, const Vec2d& v3, const Vec2d& pt) {
188 const double d1 = cross2(v1, pt);
189 const double d2 = cross2(v2, pt);
190 const double d3 = cross2(v3, pt);
191 // Testing both CCW and CW orientations.
192 return (d1 >= 0. && d2 >= 0. && d3 >= 0.) || (d1 <= 0. && d2 <= 0. && d3 <= 0.);
193 };
194 bool inside;
195 switch (proj_axis) {
196 case 0:
197 inside = inside_triangle({v1.y(), v1.z()}, {v2.y(), v2.z()}, {v3.y(), v3.z()}, {foot_pt.y(), foot_pt.z()}); break;
198 case 1:
199 inside = inside_triangle({v1.z(), v1.x()}, {v2.z(), v2.x()}, {v3.z(), v3.x()}, {foot_pt.z(), foot_pt.x()}); break;
200 default:
201 assert(proj_axis == 2);
202 inside = inside_triangle({v1.x(), v1.y()}, {v2.x(), v2.y()}, {v3.x(), v3.y()}, {foot_pt.x(), foot_pt.y()}); break;
203 }
204 if (inside)
205 return (p - foot_pt).squaredNorm();
206 }
207
208 // 4) Find minimum distance to triangle vertices and edges.
209 out = std::min((p - a).squaredNorm(), std::min((p - b).squaredNorm(), (p - c).squaredNorm()));
210 auto t = (p - a).dot(v1);
211 if (t > 0. && t < l1)
212 out = std::min(out, (a + v1 * (t / l1) - p).squaredNorm());
213 t = (p - b).dot(v2);
214 if (t > 0. && t < l2)
215 out = std::min(out, (b + v2 * (t / l2) - p).squaredNorm());
216 t = (p - c).dot(v3);
217 if (t > 0. && t < l3)
218 out = std::min(out, (c + v3 * (t / l3) - p).squaredNorm());
219 }
220
221 return out;
222 }
223
224 // Ordering of children cubes.
225 static const std::array<Vec3d, 8> child_centers {
226 Vec3d(-1, -1, -1), Vec3d( 1, -1, -1), Vec3d(-1, 1, -1), Vec3d( 1, 1, -1),
227 Vec3d(-1, -1, 1), Vec3d( 1, -1, 1), Vec3d(-1, 1, 1), Vec3d( 1, 1, 1)
228 };
229
230 // Traversal order of octree children cells for three infill directions,
231 // so that a single line will be discretized in a strictly monotonic order.
232 static constexpr std::array<std::array<int, 8>, 3> child_traversal_order {
233 std::array<int, 8>{ 2, 3, 0, 1, 6, 7, 4, 5 },
234 std::array<int, 8>{ 4, 0, 6, 2, 5, 1, 7, 3 },
235 std::array<int, 8>{ 1, 5, 0, 4, 3, 7, 2, 6 },
236 };
237
238 struct Cube
239 {
240 Vec3d center;
241 #ifndef NDEBUG
242 Vec3d center_octree;
243 #endif // NDEBUG
244 std::array<Cube*, 8> children {}; // initialized to nullptrs
CubeSlic3r::FillAdaptive::Cube245 Cube(const Vec3d ¢er) : center(center) {}
246 };
247
248 struct CubeProperties
249 {
250 double edge_length; // Lenght of edge of a cube
251 double height; // Height of rotated cube (standing on the corner)
252 double diagonal_length; // Length of diagonal of a cube a face
253 double line_z_distance; // Defines maximal distance from a center of a cube on Z axis on which lines will be created
254 double line_xy_distance;// Defines maximal distance from a center of a cube on X and Y axis on which lines will be created
255 };
256
257 struct Octree
258 {
259 // Octree will allocate its Cubes from the pool. The pool only supports deletion of the complete pool,
260 // perfect for building up our octree.
261 boost::object_pool<Cube> pool;
262 Cube* root_cube { nullptr };
263 Vec3d origin;
264 std::vector<CubeProperties> cubes_properties;
265
OctreeSlic3r::FillAdaptive::Octree266 Octree(const Vec3d &origin, const std::vector<CubeProperties> &cubes_properties)
267 : root_cube(pool.construct(origin)), origin(origin), cubes_properties(cubes_properties) {}
268
269 void insert_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, Cube *current_cube, const BoundingBoxf3 ¤t_bbox, int depth);
270 };
271
operator ()(Octree * p)272 void OctreeDeleter::operator()(Octree *p) {
273 delete p;
274 }
275
adaptive_fill_line_spacing(const PrintObject & print_object)276 std::pair<double, double> adaptive_fill_line_spacing(const PrintObject &print_object)
277 {
278 // Output, spacing for icAdaptiveCubic and icSupportCubic
279 double adaptive_line_spacing = 0.;
280 double support_line_spacing = 0.;
281
282 enum class Tristate {
283 Yes,
284 No,
285 Maybe
286 };
287 struct RegionFillData {
288 Tristate has_adaptive_infill;
289 Tristate has_support_infill;
290 double density;
291 double extrusion_width;
292 };
293 std::vector<RegionFillData> region_fill_data;
294 region_fill_data.reserve(print_object.print()->regions().size());
295 bool build_octree = false;
296 const std::vector<double> &nozzle_diameters = print_object.print()->config().nozzle_diameter.values;
297 double max_nozzle_diameter = *std::max_element(nozzle_diameters.begin(), nozzle_diameters.end());
298 double default_infill_extrusion_width = Flow::auto_extrusion_width(FlowRole::frInfill, float(max_nozzle_diameter));
299 for (const PrintRegion *region : print_object.print()->regions()) {
300 const PrintRegionConfig &config = region->config();
301 bool nonempty = config.fill_density > 0;
302 bool has_adaptive_infill = nonempty && config.fill_pattern == ipAdaptiveCubic;
303 bool has_support_infill = nonempty && config.fill_pattern == ipSupportCubic;
304 double infill_extrusion_width = config.infill_extrusion_width.percent ? default_infill_extrusion_width * 0.01 * config.infill_extrusion_width : config.infill_extrusion_width;
305 region_fill_data.push_back(RegionFillData({
306 has_adaptive_infill ? Tristate::Maybe : Tristate::No,
307 has_support_infill ? Tristate::Maybe : Tristate::No,
308 config.fill_density,
309 infill_extrusion_width != 0. ? infill_extrusion_width : default_infill_extrusion_width
310 }));
311 build_octree |= has_adaptive_infill || has_support_infill;
312 }
313
314 if (build_octree) {
315 // Compute the average of above parameters over all layers
316 for (const Layer *layer : print_object.layers())
317 for (size_t region_id = 0; region_id < layer->regions().size(); ++ region_id) {
318 RegionFillData &rd = region_fill_data[region_id];
319 if (rd.has_adaptive_infill == Tristate::Maybe && ! layer->regions()[region_id]->fill_surfaces.empty())
320 rd.has_adaptive_infill = Tristate::Yes;
321 if (rd.has_support_infill == Tristate::Maybe && ! layer->regions()[region_id]->fill_surfaces.empty())
322 rd.has_support_infill = Tristate::Yes;
323 }
324
325 double adaptive_fill_density = 0.;
326 double adaptive_infill_extrusion_width = 0.;
327 int adaptive_cnt = 0;
328 double support_fill_density = 0.;
329 double support_infill_extrusion_width = 0.;
330 int support_cnt = 0;
331
332 for (const RegionFillData &rd : region_fill_data) {
333 if (rd.has_adaptive_infill == Tristate::Yes) {
334 adaptive_fill_density += rd.density;
335 adaptive_infill_extrusion_width += rd.extrusion_width;
336 ++ adaptive_cnt;
337 } else if (rd.has_support_infill == Tristate::Yes) {
338 support_fill_density += rd.density;
339 support_infill_extrusion_width += rd.extrusion_width;
340 ++ support_cnt;
341 }
342 }
343
344 auto to_line_spacing = [](int cnt, double density, double extrusion_width) {
345 if (cnt) {
346 density /= double(cnt);
347 extrusion_width /= double(cnt);
348 return extrusion_width / ((density / 100.0f) * 0.333333333f);
349 } else
350 return 0.;
351 };
352 adaptive_line_spacing = to_line_spacing(adaptive_cnt, adaptive_fill_density, adaptive_infill_extrusion_width);
353 support_line_spacing = to_line_spacing(support_cnt, support_fill_density, support_infill_extrusion_width);
354 }
355
356 return std::make_pair(adaptive_line_spacing, support_line_spacing);
357 }
358
359 // Context used by generate_infill_lines() when recursively traversing an octree in a DDA fashion
360 // (Digital Differential Analyzer).
361 struct FillContext
362 {
363 // The angles have to agree with child_traversal_order.
364 static constexpr double direction_angles[3] {
365 0.,
366 (2.0 * M_PI) / 3.0,
367 -(2.0 * M_PI) / 3.0
368 };
369
FillContextSlic3r::FillAdaptive::FillContext370 FillContext(const Octree &octree, double z_position, int direction_idx) :
371 cubes_properties(octree.cubes_properties),
372 z_position(z_position),
373 traversal_order(child_traversal_order[direction_idx]),
374 cos_a(cos(direction_angles[direction_idx])),
375 sin_a(sin(direction_angles[direction_idx]))
376 {
377 static constexpr auto unused = std::numeric_limits<coord_t>::max();
378 temp_lines.assign((1 << octree.cubes_properties.size()) - 1, Line(Point(unused, unused), Point(unused, unused)));
379 }
380
381 // Rotate the point, uses the same convention as Point::rotate().
rotateSlic3r::FillAdaptive::FillContext382 Vec2d rotate(const Vec2d& v) { return Vec2d(this->cos_a * v.x() - this->sin_a * v.y(), this->sin_a * v.x() + this->cos_a * v.y()); }
383
384 const std::vector<CubeProperties> &cubes_properties;
385 // Top of the current layer.
386 const double z_position;
387 // Order of traversal for this line direction.
388 const std::array<int, 8> traversal_order;
389 // Rotation of the generated line for this line direction.
390 const double cos_a;
391 const double sin_a;
392
393 // Linearized tree spanning a single Octree wall, used to connect lines spanning
394 // neighboring Octree cells. Unused lines have the Line::a::x set to infinity.
395 std::vector<Line> temp_lines;
396 // Final output
397 std::vector<Line> output_lines;
398 };
399
400 static constexpr double octree_rot[3] = { 5.0 * M_PI / 4.0, Geometry::deg2rad(215.264), M_PI / 6.0 };
401
transform_to_world()402 Eigen::Quaterniond transform_to_world()
403 {
404 return Eigen::AngleAxisd(octree_rot[2], Vec3d::UnitZ()) * Eigen::AngleAxisd(octree_rot[1], Vec3d::UnitY()) * Eigen::AngleAxisd(octree_rot[0], Vec3d::UnitX());
405 }
406
transform_to_octree()407 Eigen::Quaterniond transform_to_octree()
408 {
409 return Eigen::AngleAxisd(- octree_rot[0], Vec3d::UnitX()) * Eigen::AngleAxisd(- octree_rot[1], Vec3d::UnitY()) * Eigen::AngleAxisd(- octree_rot[2], Vec3d::UnitZ());
410 }
411
412 #ifndef NDEBUG
413 // Verify that the traversal order of the octree children matches the line direction,
414 // therefore the infill line may get extended with O(1) time & space complexity.
verify_traversal_order(FillContext & context,const Cube * cube,int depth,const Vec2d & line_from,const Vec2d & line_to)415 static bool verify_traversal_order(
416 FillContext &context,
417 const Cube *cube,
418 int depth,
419 const Vec2d &line_from,
420 const Vec2d &line_to)
421 {
422 std::array<Vec3d, 8> c;
423 Eigen::Quaterniond to_world = transform_to_world();
424 for (int i = 0; i < 8; ++i) {
425 int j = context.traversal_order[i];
426 Vec3d cntr = to_world * (cube->center_octree + (child_centers[j] * (context.cubes_properties[depth].edge_length / 4.)));
427 assert(!cube->children[j] || cube->children[j]->center.isApprox(cntr));
428 c[i] = cntr;
429 }
430 std::array<Vec3d, 10> dirs = {
431 c[1] - c[0], c[2] - c[0], c[3] - c[1], c[3] - c[2], c[3] - c[0],
432 c[5] - c[4], c[6] - c[4], c[7] - c[5], c[7] - c[6], c[7] - c[4]
433 };
434 assert(std::abs(dirs[4].z()) < 0.005);
435 assert(std::abs(dirs[9].z()) < 0.005);
436 assert(dirs[0].isApprox(dirs[3]));
437 assert(dirs[1].isApprox(dirs[2]));
438 assert(dirs[5].isApprox(dirs[8]));
439 assert(dirs[6].isApprox(dirs[7]));
440 Vec3d line_dir = Vec3d(line_to.x() - line_from.x(), line_to.y() - line_from.y(), 0.).normalized();
441 for (auto& dir : dirs) {
442 double d = dir.normalized().dot(line_dir);
443 assert(d > 0.7);
444 }
445 return true;
446 }
447 #endif // NDEBUG
448
generate_infill_lines_recursive(FillContext & context,const Cube * cube,int address,int depth)449 static void generate_infill_lines_recursive(
450 FillContext &context,
451 const Cube *cube,
452 // Address of this wall in the octree, used to address context.temp_lines.
453 int address,
454 int depth)
455 {
456 assert(cube != nullptr);
457
458 const std::vector<CubeProperties> &cubes_properties = context.cubes_properties;
459 const double z_diff = context.z_position - cube->center.z();
460 const double z_diff_abs = std::abs(z_diff);
461
462 if (z_diff_abs > cubes_properties[depth].height / 2.)
463 return;
464
465 if (z_diff_abs < cubes_properties[depth].line_z_distance) {
466 // Discretize a single wall splitting the cube into two.
467 const double zdist = cubes_properties[depth].line_z_distance;
468 Vec2d from(
469 0.5 * cubes_properties[depth].diagonal_length * (zdist - z_diff_abs) / zdist,
470 cubes_properties[depth].line_xy_distance - (zdist + z_diff) / sqrt(2.));
471 Vec2d to(-from.x(), from.y());
472 from = context.rotate(from);
473 to = context.rotate(to);
474 // Relative to cube center
475 const Vec2d offset(cube->center.x(), cube->center.y());
476 from += offset;
477 to += offset;
478 // Verify that the traversal order of the octree children matches the line direction,
479 // therefore the infill line may get extended with O(1) time & space complexity.
480 assert(verify_traversal_order(context, cube, depth, from, to));
481 // Either extend an existing line or start a new one.
482 Line &last_line = context.temp_lines[address];
483 Line new_line(Point::new_scale(from), Point::new_scale(to));
484 if (last_line.a.x() == std::numeric_limits<coord_t>::max()) {
485 last_line.a = new_line.a;
486 } else if ((new_line.a - last_line.b).cwiseAbs().maxCoeff() > 1000) { // SCALED_EPSILON is 100 and it is not enough
487 context.output_lines.emplace_back(last_line);
488 last_line.a = new_line.a;
489 }
490 last_line.b = new_line.b;
491 }
492
493 // left child index
494 address = address * 2 + 1;
495 -- depth;
496 size_t i = 0;
497 for (const int child_idx : context.traversal_order) {
498 const Cube *child = cube->children[child_idx];
499 if (child != nullptr)
500 generate_infill_lines_recursive(context, child, address, depth);
501 if (++ i == 4)
502 // right child index
503 ++ address;
504 }
505 }
506
507 #ifndef NDEBUG
508 // #define ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
509 #endif
510
511 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
export_infill_lines_to_svg(const ExPolygon & expoly,const Polylines & polylines,const std::string & path,const Points & pts=Points ())512 static void export_infill_lines_to_svg(const ExPolygon &expoly, const Polylines &polylines, const std::string &path, const Points &pts = Points())
513 {
514 BoundingBox bbox = get_extents(expoly);
515 bbox.offset(scale_(3.));
516
517 ::Slic3r::SVG svg(path, bbox);
518 svg.draw(expoly);
519 svg.draw_outline(expoly, "green");
520 svg.draw(polylines, "red");
521 static constexpr double trim_length = scale_(0.4);
522 for (Polyline polyline : polylines)
523 if (! polyline.empty()) {
524 Vec2d a = polyline.points.front().cast<double>();
525 Vec2d d = polyline.points.back().cast<double>();
526 if (polyline.size() == 2) {
527 Vec2d v = d - a;
528 double l = v.norm();
529 if (l > 2. * trim_length) {
530 a += v * trim_length / l;
531 d -= v * trim_length / l;
532 polyline.points.front() = a.cast<coord_t>();
533 polyline.points.back() = d.cast<coord_t>();
534 } else
535 polyline.points.clear();
536 } else if (polyline.size() > 2) {
537 Vec2d b = polyline.points[1].cast<double>();
538 Vec2d c = polyline.points[polyline.points.size() - 2].cast<double>();
539 Vec2d v = b - a;
540 double l = v.norm();
541 if (l > trim_length) {
542 a += v * trim_length / l;
543 polyline.points.front() = a.cast<coord_t>();
544 } else
545 polyline.points.erase(polyline.points.begin());
546 v = d - c;
547 l = v.norm();
548 if (l > trim_length)
549 polyline.points.back() = (d - v * trim_length / l).cast<coord_t>();
550 else
551 polyline.points.pop_back();
552 }
553 svg.draw(polyline, "black");
554 }
555 svg.draw(pts, "magenta");
556 }
557 #endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
558
559 // Representing a T-joint (in general case) between two infill lines
560 // (between one end point of intersect_pl/intersect_line and
561 struct Intersection
562 {
563 // Closest line to intersect_point.
564 const Line *closest_line;
565
566 // The line for which is computed closest line from intersect_point to closest_line
567 const Line *intersect_line;
568 // Pointer to the polyline from which is computed closest_line
569 Polyline *intersect_pl;
570 // Point for which is computed closest line (closest_line)
571 Point intersect_point;
572 // Indicate if intersect_point is the first or the last point of intersect_pl
573 bool front;
574 // Signum of intersect_line_dir.cross(closest_line.dir()):
575 bool left;
576
577 // Indication if this intersection has been proceed
578 bool used = false;
579
freshSlic3r::FillAdaptive::Intersection580 bool fresh() const throw() { return ! used && ! intersect_pl->empty(); }
581
IntersectionSlic3r::FillAdaptive::Intersection582 Intersection(const Line &closest_line, const Line &intersect_line, Polyline *intersect_pl, const Point &intersect_point, bool front) :
583 closest_line(&closest_line), intersect_line(&intersect_line), intersect_pl(intersect_pl), intersect_point(intersect_point), front(front)
584 {
585 // Calculate side of this intersection line of the closest line.
586 Vec2d v1((this->closest_line->b - this->closest_line->a).cast<double>());
587 Vec2d v2(this->intersect_line_dir());
588 #ifndef NDEBUG
589 {
590 Vec2d v1n = v1.normalized();
591 Vec2d v2n = v2.normalized();
592 double c = cross2(v1n, v2n);
593 assert(std::abs(c) > sin(M_PI / 12.));
594 }
595 #endif // NDEBUG
596 this->left = cross2(v1, v2) > 0.;
597 }
598
other_hookSlic3r::FillAdaptive::Intersection599 std::optional<Line> other_hook() const {
600 std::optional<Line> out;
601 const Points &pts = intersect_pl->points;
602 if (pts.size() >= 3)
603 out = this->front ? Line(pts[1], pts[2]) : Line(pts[pts.size() - 2], pts[pts.size() - 3]);
604 return out;
605 }
606
other_hook_intersectsSlic3r::FillAdaptive::Intersection607 bool other_hook_intersects(const Line &l, Point &pt) {
608 std::optional<Line> h = other_hook();
609 return h && h->intersection(l, &pt);
610 }
other_hook_intersectsSlic3r::FillAdaptive::Intersection611 bool other_hook_intersects(const Line &l) { Point pt; return this->other_hook_intersects(l, pt); }
612
613 // Direction to intersect_point.
intersect_line_dirSlic3r::FillAdaptive::Intersection614 Vec2d intersect_line_dir() const throw() {
615 return (this->intersect_point == intersect_line->a ? intersect_line->b - intersect_line->a : intersect_line->a - intersect_line->b).cast<double>();
616 }
617 };
618
get_nearest_intersection(std::vector<std::pair<Intersection *,double>> & intersect_line,const size_t first_idx)619 static inline Intersection* get_nearest_intersection(std::vector<std::pair<Intersection*, double>>& intersect_line, const size_t first_idx)
620 {
621 assert(intersect_line.size() >= 2);
622 bool take_next = false;
623 if (first_idx == 0)
624 take_next = true;
625 else if (first_idx + 1 == intersect_line.size())
626 take_next = false;
627 else {
628 // Has both prev and next.
629 const std::pair<Intersection*, double> &ithis = intersect_line[first_idx];
630 const std::pair<Intersection*, double> &iprev = intersect_line[first_idx - 1];
631 const std::pair<Intersection*, double> &inext = intersect_line[first_idx + 1];
632 take_next = iprev.first->fresh() && inext.first->fresh() ?
633 inext.second - ithis.second < ithis.second - iprev.second :
634 inext.first->fresh();
635 }
636 return intersect_line[take_next ? first_idx + 1 : first_idx - 1].first;
637 }
638
639 // Create a line representing the anchor aka hook extrusion based on line_to_offset
640 // translated in the direction of the intersection line (intersection.intersect_line).
create_offset_line(Line offset_line,const Intersection & intersection,const double scaled_offset)641 static Line create_offset_line(Line offset_line, const Intersection &intersection, const double scaled_offset)
642 {
643 offset_line.translate((perp(intersection.closest_line->vector().cast<double>().normalized()) * (intersection.left ? scaled_offset : - scaled_offset)).cast<coord_t>());
644 // Extend the line by a small value to guarantee a collision with adjacent lines
645 offset_line.extend(coord_t(scaled_offset * 1.16)); // / cos(PI/6)
646 return offset_line;
647 }
648
649 namespace bg = boost::geometry;
650 namespace bgm = boost::geometry::model;
651 namespace bgi = boost::geometry::index;
652
653 // float is needed because for coord_t bgi::intersects throws "bad numeric conversion: positive overflow"
654 using rtree_point_t = bgm::point<float, 2, boost::geometry::cs::cartesian>;
655 using rtree_segment_t = bgm::segment<rtree_point_t>;
656 using rtree_t = bgi::rtree<std::pair<rtree_segment_t, size_t>, bgi::rstar<16, 4>>;
657
mk_rtree_point(const Point & pt)658 static inline rtree_point_t mk_rtree_point(const Point &pt) {
659 return rtree_point_t(float(pt.x()), float(pt.y()));
660 }
mk_rtree_seg(const Point & a,const Point & b)661 static inline rtree_segment_t mk_rtree_seg(const Point &a, const Point &b) {
662 return { mk_rtree_point(a), mk_rtree_point(b) };
663 }
mk_rtree_seg(const Line & l)664 static inline rtree_segment_t mk_rtree_seg(const Line &l) {
665 return mk_rtree_seg(l.a, l.b);
666 }
667
668 // Create a hook based on hook_line and append it to the begin or end of the polyline in the intersection
add_hook(const Intersection & intersection,const double scaled_offset,const coordf_t hook_length,double scaled_trim_distance,const rtree_t & rtree,const Lines & lines_src)669 static void add_hook(
670 const Intersection &intersection, const double scaled_offset,
671 const coordf_t hook_length, double scaled_trim_distance,
672 const rtree_t &rtree, const Lines &lines_src)
673 {
674 if (hook_length < SCALED_EPSILON)
675 // Ignore open hooks.
676 return;
677
678 #ifndef NDEBUG
679 {
680 const Vec2d v = (intersection.closest_line->b - intersection.closest_line->a).cast<double>();
681 const Vec2d va = (intersection.intersect_point - intersection.closest_line->a).cast<double>();
682 const double l2 = v.squaredNorm(); // avoid a sqrt
683 assert(l2 > 0.);
684 const double t = va.dot(v) / l2;
685 assert(t > 0. && t < 1.);
686 const double d = (t * v - va).norm();
687 assert(d < 1000.);
688 }
689 #endif // NDEBUG
690
691 // Trim the hook start by the infill line it will connect to.
692 Point hook_start;
693 bool intersection_found = intersection.intersect_line->intersection(
694 create_offset_line(*intersection.closest_line, intersection, scaled_offset),
695 &hook_start);
696 assert(intersection_found);
697
698 std::optional<Line> other_hook = intersection.other_hook();
699
700 Vec2d hook_vector_norm = intersection.closest_line->vector().cast<double>().normalized();
701 // hook_vector is extended by the thickness of the infill line, so that a collision is found against
702 // the infill centerline to be later trimmed by the thickened line.
703 Vector hook_vector = ((hook_length + 1.16 * scaled_trim_distance) * hook_vector_norm).cast<coord_t>();
704 Line hook_forward(hook_start, hook_start + hook_vector);
705
706 auto filter_itself = [&intersection, &lines_src](const auto &item) { return item.second != intersection.intersect_line - lines_src.data(); };
707
708 std::vector<std::pair<rtree_segment_t, size_t>> hook_intersections;
709 rtree.query(bgi::intersects(mk_rtree_seg(hook_forward)) && bgi::satisfies(filter_itself), std::back_inserter(hook_intersections));
710 Point self_intersection_point;
711 bool self_intersection = other_hook && other_hook->intersection(hook_forward, &self_intersection_point);
712
713 // Find closest intersection of a line segment starting with pt pointing in dir
714 // with any of the hook_intersections, returns Euclidian distance.
715 // dir is normalized.
716 auto max_hook_length = [hook_length, scaled_trim_distance, &lines_src](
717 const Vec2d &pt, const Vec2d &dir,
718 const std::vector<std::pair<rtree_segment_t, size_t>> &hook_intersections,
719 bool self_intersection, const std::optional<Line> &self_intersection_line, const Point &self_intersection_point) {
720 // No hook is longer than hook_length, there shouldn't be any intersection closer than that.
721 auto max_length = hook_length;
722 auto update_max_length = [&max_length](double d) {
723 if (d < max_length)
724 max_length = d;
725 };
726 // Shift the trimming point away from the colliding thick line.
727 auto shift_from_thick_line = [&dir, scaled_trim_distance](const Vec2d& dir2) {
728 return scaled_trim_distance * std::abs(cross2(dir, dir2.normalized()));
729 };
730
731 for (const auto &hook_intersection : hook_intersections) {
732 const rtree_segment_t &segment = hook_intersection.first;
733 // Segment start and end points, segment vector.
734 Vec2d pt2(bg::get<0, 0>(segment), bg::get<0, 1>(segment));
735 Vec2d dir2 = Vec2d(bg::get<1, 0>(segment), bg::get<1, 1>(segment)) - pt2;
736 // Find intersection of (pt, dir) with (pt2, dir2), where dir is normalized.
737 double denom = cross2(dir, dir2);
738 assert(std::abs(denom) > EPSILON);
739 double t = cross2(pt2 - pt, dir2) / denom;
740 if (hook_intersection.second < lines_src.size())
741 // Trimming by another infill line. Reduce overlap.
742 t -= shift_from_thick_line(dir2);
743 update_max_length(t);
744 }
745 if (self_intersection) {
746 double t = (self_intersection_point.cast<double>() - pt).dot(dir) - shift_from_thick_line((*self_intersection_line).vector().cast<double>());
747 max_length = std::min(max_length, t);
748 }
749 return std::max(0., max_length);
750 };
751
752 Vec2d hook_startf = hook_start.cast<double>();
753 double hook_forward_max_length = max_hook_length(hook_startf, hook_vector_norm, hook_intersections, self_intersection, other_hook, self_intersection_point);
754 double hook_backward_max_length = 0.;
755 if (hook_forward_max_length < hook_length - SCALED_EPSILON) {
756 // Try the other side.
757 hook_intersections.clear();
758 Line hook_backward(hook_start, hook_start - hook_vector);
759 rtree.query(bgi::intersects(mk_rtree_seg(hook_backward)) && bgi::satisfies(filter_itself), std::back_inserter(hook_intersections));
760 self_intersection = other_hook && other_hook->intersection(hook_backward, &self_intersection_point);
761 hook_backward_max_length = max_hook_length(hook_startf, - hook_vector_norm, hook_intersections, self_intersection, other_hook, self_intersection_point);
762 }
763
764 // Take the longer hook.
765 Vec2d hook_dir = (hook_forward_max_length > hook_backward_max_length ? hook_forward_max_length : - hook_backward_max_length) * hook_vector_norm;
766 Point hook_end = hook_start + hook_dir.cast<coord_t>();
767
768 Points &pl = intersection.intersect_pl->points;
769 if (intersection.front) {
770 pl.front() = hook_start;
771 pl.emplace(pl.begin(), hook_end);
772 } else {
773 pl.back() = hook_start;
774 pl.emplace_back(hook_end);
775 }
776 }
777
778 #ifndef NDEBUG
validate_intersection_t_joint(const Intersection & intersection)779 bool validate_intersection_t_joint(const Intersection &intersection)
780 {
781 const Vec2d v = (intersection.closest_line->b - intersection.closest_line->a).cast<double>();
782 const Vec2d va = (intersection.intersect_point - intersection.closest_line->a).cast<double>();
783 const double l2 = v.squaredNorm(); // avoid a sqrt
784 assert(l2 > 0.);
785 const double t = va.dot(v);
786 assert(t > SCALED_EPSILON && t < l2 - SCALED_EPSILON);
787 const double d = ((t / l2) * v - va).norm();
788 assert(d < 1000.);
789 return true;
790 }
validate_intersections(const std::vector<Intersection> & intersections)791 bool validate_intersections(const std::vector<Intersection> &intersections)
792 {
793 for (const Intersection& intersection : intersections)
794 assert(validate_intersection_t_joint(intersection));
795 return true;
796 }
797 #endif // NDEBUG
798
connect_lines_using_hooks(Polylines && lines,const ExPolygon & boundary,const double spacing,const coordf_t hook_length,const coordf_t hook_length_max)799 static Polylines connect_lines_using_hooks(Polylines &&lines, const ExPolygon &boundary, const double spacing, const coordf_t hook_length, const coordf_t hook_length_max)
800 {
801 rtree_t rtree;
802 size_t poly_idx = 0;
803
804 // 19% overlap, slightly lower than the allowed overlap in Fill::connect_infill()
805 const float scaled_offset = float(scale_(spacing) * 0.81);
806 // 25% overlap
807 const float scaled_trim_distance = float(scale_(spacing) * 0.5 * 0.75);
808
809 // Keeping the vector of closest points outside the loop, so the vector does not need to be reallocated.
810 std::vector<std::pair<rtree_segment_t, size_t>> closest;
811 // Pairs of lines touching at one end point. The pair is sorted to make the end point connection test symmetric.
812 std::vector<std::pair<const Polyline*, const Polyline*>> lines_touching_at_endpoints;
813 {
814 // Insert infill lines into rtree, merge close collinear segments split by the infill boundary,
815 // collect lines_touching_at_endpoints.
816 double r2_close = Slic3r::sqr(1200.);
817 for (Polyline &poly : lines) {
818 assert(poly.points.size() == 2);
819 if (&poly != lines.data()) {
820 // Join collinear segments separated by a tiny gap. These gaps were likely created by clipping the infill lines with a concave dent in an infill boundary.
821 auto collinear_segment = [&rtree, &closest, &lines, &lines_touching_at_endpoints, r2_close](const Point& pt, const Point& pt_other, const Polyline* polyline) -> std::pair<Polyline*, bool> {
822 closest.clear();
823 rtree.query(bgi::nearest(mk_rtree_point(pt), 1), std::back_inserter(closest));
824 const Polyline *other = &lines[closest.front().second];
825 double dist2_front = (other->points.front() - pt).cast<double>().squaredNorm();
826 double dist2_back = (other->points.back() - pt).cast<double>().squaredNorm();
827 double dist2_min = std::min(dist2_front, dist2_back);
828 if (dist2_min < r2_close) {
829 // Don't connect the segments in an opposite direction.
830 double dist2_min_other = std::min((other->points.front() - pt_other).cast<double>().squaredNorm(), (other->points.back() - pt_other).cast<double>().squaredNorm());
831 if (dist2_min_other > dist2_min) {
832 // End points of the two lines are very close, they should have been merged together if they are collinear.
833 Vec2d v1 = (pt_other - pt).cast<double>();
834 Vec2d v2 = (other->points.back() - other->points.front()).cast<double>();
835 Vec2d v1n = v1.normalized();
836 Vec2d v2n = v2.normalized();
837 // The vectors must not be collinear.
838 double d = v1n.dot(v2n);
839 if (std::abs(d) > 0.99f) {
840 // Lines are collinear, merge them.
841 rtree.remove(closest.front());
842 return std::make_pair(const_cast<Polyline*>(other), dist2_min == dist2_front);
843 } else {
844 if (polyline > other)
845 std::swap(polyline, other);
846 lines_touching_at_endpoints.emplace_back(polyline, other);
847 }
848 }
849 }
850 return std::make_pair(static_cast<Polyline*>(nullptr), false);
851 };
852 auto collinear_front = collinear_segment(poly.points.front(), poly.points.back(), &poly);
853 auto collinear_back = collinear_segment(poly.points.back(), poly.points.front(), &poly);
854 assert(! collinear_front.first || ! collinear_back.first || collinear_front.first != collinear_back.first);
855 if (collinear_front.first) {
856 Polyline &other = *collinear_front.first;
857 assert(&other != &poly);
858 poly.points.front() = collinear_front.second ? other.points.back() : other.points.front();
859 other.points.clear();
860 }
861 if (collinear_back.first) {
862 Polyline &other = *collinear_back.first;
863 assert(&other != &poly);
864 poly.points.back() = collinear_back.second ? other.points.back() : other.points.front();
865 other.points.clear();
866 }
867 }
868 rtree.insert(std::make_pair(mk_rtree_seg(poly.points.front(), poly.points.back()), poly_idx++));
869 }
870 }
871
872 // Convert input polylines to lines_src after the colinear segments were merged.
873 Lines lines_src;
874 lines_src.reserve(lines.size());
875 std::transform(lines.begin(), lines.end(), std::back_inserter(lines_src), [](const Polyline &pl) {
876 return pl.empty() ? Line(Point(0, 0), Point(0, 0)) : Line(pl.points.front(), pl.points.back()); });
877
878 sort_remove_duplicates(lines_touching_at_endpoints);
879
880 std::vector<Intersection> intersections;
881 {
882 // Minimum lenght of an infill line to anchor. Very short lines cannot be trimmed from both sides,
883 // it does not help to anchor extremely short infill lines, it consumes too much plastic while not adding
884 // to the object rigidity.
885 assert(scaled_offset > scaled_trim_distance);
886 const double line_len_threshold_drop_both_sides = scaled_offset * (2. / cos(PI / 6.) + 0.5) + SCALED_EPSILON;
887 const double line_len_threshold_anchor_both_sides = line_len_threshold_drop_both_sides + scaled_offset;
888 const double line_len_threshold_drop_single_side = scaled_offset * (1. / cos(PI / 6.) + 1.5) + SCALED_EPSILON;
889 const double line_len_threshold_anchor_single_side = line_len_threshold_drop_single_side + scaled_offset;
890 for (size_t line_idx = 0; line_idx < lines.size(); ++ line_idx) {
891 Polyline &line = lines[line_idx];
892 if (line.points.empty())
893 continue;
894
895 Point &front_point = line.points.front();
896 Point &back_point = line.points.back();
897
898 // Find the nearest line from the start point of the line.
899 std::optional<size_t> tjoint_front, tjoint_back;
900 {
901 auto has_tjoint = [&closest, line_idx, &rtree, &lines, &lines_src](const Point &pt) {
902 auto filter_t_joint = [line_idx, &lines_src, pt](const auto &item) {
903 if (item.second != line_idx) {
904 // Verify that the point projects onto the line.
905 const Line &line = lines_src[item.second];
906 const Vec2d v = (line.b - line.a).cast<double>();
907 const Vec2d va = (pt - line.a).cast<double>();
908 const double l2 = v.squaredNorm(); // avoid a sqrt
909 if (l2 > 0.) {
910 const double t = va.dot(v);
911 return t > SCALED_EPSILON && t < l2 - SCALED_EPSILON;
912 }
913 }
914 return false;
915 };
916 closest.clear();
917 rtree.query(bgi::nearest(mk_rtree_point(pt), 1) && bgi::satisfies(filter_t_joint), std::back_inserter(closest));
918 std::optional<size_t> out;
919 if (! closest.empty()) {
920 const Polyline &pl = lines[closest.front().second];
921 if (pl.points.empty()) {
922 // The closest infill line was already dropped as it was too short.
923 // Such an infill line should not make a T-joint anyways.
924 #if 0 // #ifndef NDEBUG
925 const auto &seg = closest.front().first;
926 struct Linef { Vec2d a; Vec2d b; };
927 Linef l { { bg::get<0, 0>(seg), bg::get<0, 1>(seg) }, { bg::get<1, 0>(seg), bg::get<1, 1>(seg) } };
928 assert(line_alg::distance_to_squared(l, Vec2d(pt.cast<double>())) > 1000 * 1000);
929 #endif // NDEBUG
930 } else if (((Line)pl).distance_to_squared(pt) <= 1000 * 1000)
931 out = closest.front().second;
932 }
933 return out;
934 };
935 // Refuse to create a T-joint if the infill lines touch at their ends.
936 auto filter_end_point_connections = [&lines_touching_at_endpoints, &lines, &line](std::optional<size_t> in) {
937 std::optional<size_t> out;
938 if (in) {
939 const Polyline *lo = &line;
940 const Polyline *hi = &lines[*in];
941 if (lo > hi)
942 std::swap(lo, hi);
943 if (! std::binary_search(lines_touching_at_endpoints.begin(), lines_touching_at_endpoints.end(), std::make_pair(lo, hi)))
944 // Not an end-point connection, it is a valid T-joint.
945 out = in;
946 }
947 return out;
948 };
949 tjoint_front = filter_end_point_connections(has_tjoint(front_point));
950 tjoint_back = filter_end_point_connections(has_tjoint(back_point));
951 }
952
953 int num_tjoints = int(tjoint_front.has_value()) + int(tjoint_back.has_value());
954 if (num_tjoints > 0) {
955 double line_len = line.length();
956 bool drop = false;
957 bool anchor = false;
958 if (num_tjoints == 1) {
959 // Connected to perimeters on a single side only, connected to another infill line on the other side.
960 drop = line_len < line_len_threshold_drop_single_side;
961 anchor = line_len > line_len_threshold_anchor_single_side;
962 } else {
963 // Not connected to perimeters at all, connected to two infill lines.
964 assert(num_tjoints == 2);
965 drop = line_len < line_len_threshold_drop_both_sides;
966 anchor = line_len > line_len_threshold_anchor_both_sides;
967 }
968 if (drop) {
969 // Drop a very short line if connected to another infill line.
970 // Lines shorter than spacing are skipped because it is needed to shrink a line by the value of spacing.
971 // A shorter line than spacing could produce a degenerate polyline.
972 line.points.clear();
973 } else if (anchor) {
974 if (tjoint_front) {
975 // T-joint of line's front point with the 'closest' line.
976 intersections.emplace_back(lines_src[*tjoint_front], lines_src[line_idx], &line, front_point, true);
977 assert(validate_intersection_t_joint(intersections.back()));
978 }
979 if (tjoint_back) {
980 // T-joint of line's back point with the 'closest' line.
981 intersections.emplace_back(lines_src[*tjoint_back], lines_src[line_idx], &line, back_point, false);
982 assert(validate_intersection_t_joint(intersections.back()));
983 }
984 } else {
985 if (tjoint_front)
986 // T joint at the front at a 60 degree angle, the line is very short.
987 // Trim the front side.
988 front_point += ((scaled_trim_distance * 1.155) * (back_point - front_point).cast<double>().normalized()).cast<coord_t>();
989 if (tjoint_back)
990 // T joint at the front at a 60 degree angle, the line is very short.
991 // Trim the front side.
992 back_point += ((scaled_trim_distance * 1.155) * (front_point - back_point).cast<double>().normalized()).cast<coord_t>();
993 }
994 }
995 }
996 // Remove those intersections, that point to a dropped line.
997 for (auto it = intersections.begin(); it != intersections.end(); ) {
998 assert(! lines[it->intersect_line - lines_src.data()].points.empty());
999 if (lines[it->closest_line - lines_src.data()].points.empty()) {
1000 *it = intersections.back();
1001 intersections.pop_back();
1002 } else
1003 ++ it;
1004 }
1005 }
1006 assert(validate_intersections(intersections));
1007
1008 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1009 static int iRun = 0;
1010 int iStep = 0;
1011 {
1012 Points pts;
1013 for (const Intersection &i : intersections)
1014 pts.emplace_back(i.intersect_point);
1015 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-Tjoints-%d.svg", iRun++), pts);
1016 }
1017 #endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
1018
1019 // Sort lexicographically by closest_line_idx and left/right orientation.
1020 std::sort(intersections.begin(), intersections.end(),
1021 [](const Intersection &i1, const Intersection &i2) {
1022 return (i1.closest_line == i2.closest_line) ?
1023 int(i1.left) < int(i2.left) :
1024 i1.closest_line < i2.closest_line;
1025 });
1026
1027 std::vector<size_t> merged_with(lines.size());
1028 std::iota(merged_with.begin(), merged_with.end(), 0);
1029
1030 // Appends the boundary polygon with all holes to rtree for detection to check whether hooks are not crossing the boundary
1031 {
1032 Point prev = boundary.contour.points.back();
1033 for (const Point &point : boundary.contour.points) {
1034 rtree.insert(std::make_pair(mk_rtree_seg(prev, point), poly_idx++));
1035 prev = point;
1036 }
1037 for (const Polygon &polygon : boundary.holes) {
1038 Point prev = polygon.points.back();
1039 for (const Point &point : polygon.points) {
1040 rtree.insert(std::make_pair(mk_rtree_seg(prev, point), poly_idx++));
1041 prev = point;
1042 }
1043 }
1044 }
1045
1046 auto update_merged_polyline_idx = [&merged_with](size_t pl_idx) {
1047 // Update the polyline index to index which is merged
1048 for (size_t last = pl_idx;;) {
1049 size_t lower = merged_with[last];
1050 if (lower == last) {
1051 merged_with[pl_idx] = lower;
1052 return lower;
1053 }
1054 last = lower;
1055 }
1056 assert(false);
1057 return size_t(0);
1058 };
1059 auto update_merged_polyline = [&lines, update_merged_polyline_idx](Intersection& intersection) {
1060 // Update the polyline index to index which is merged
1061 size_t intersect_pl_idx = update_merged_polyline_idx(intersection.intersect_pl - lines.data());
1062 intersection.intersect_pl = &lines[intersect_pl_idx];
1063 // After polylines are merged, it is necessary to update "forward" based on if intersect_point is the first or the last point of intersect_pl.
1064 if (intersection.fresh()) {
1065 assert(intersection.intersect_pl->points.front() == intersection.intersect_point ||
1066 intersection.intersect_pl->points.back() == intersection.intersect_point);
1067 intersection.front = intersection.intersect_pl->points.front() == intersection.intersect_point;
1068 }
1069 };
1070
1071 // Merge polylines touching at their ends. This should be a very rare case, but it happens surprisingly often.
1072 for (auto it = lines_touching_at_endpoints.rbegin(); it != lines_touching_at_endpoints.rend(); ++ it) {
1073 Polyline *pl1 = const_cast<Polyline*>(it->first);
1074 Polyline *pl2 = const_cast<Polyline*>(it->second);
1075 assert(pl1 < pl2);
1076 // pl1 was visited for the 1st time.
1077 // pl2 may have alread been merged with another polyline, even with this one.
1078 pl2 = &lines[update_merged_polyline_idx(pl2 - lines.data())];
1079 assert(pl1 <= pl2);
1080 // Avoid closing a loop, ignore dropped infill lines.
1081 if (pl1 != pl2 && ! pl1->points.empty() && ! pl2->points.empty()) {
1082 // Merge the polylines.
1083 assert(pl1 < pl2);
1084 assert(pl1->points.size() >= 2);
1085 assert(pl2->points.size() >= 2);
1086 double d11 = (pl1->points.front() - pl2->points.front()).cast<double>().squaredNorm();
1087 double d12 = (pl1->points.front() - pl2->points.back()) .cast<double>().squaredNorm();
1088 double d21 = (pl1->points.back() - pl2->points.front()).cast<double>().squaredNorm();
1089 double d22 = (pl1->points.back() - pl2->points.back()) .cast<double>().squaredNorm();
1090 double d1min = std::min(d11, d12);
1091 double d2min = std::min(d21, d22);
1092 if (d1min < d2min) {
1093 pl1->reverse();
1094 if (d12 == d1min)
1095 pl2->reverse();
1096 } else if (d22 == d2min)
1097 pl2->reverse();
1098 pl1->points.back() = (pl1->points.back() + pl2->points.front()) / 2;
1099 pl1->append(pl2->points.begin() + 1, pl2->points.end());
1100 pl2->points.clear();
1101 merged_with[pl2 - lines.data()] = pl1 - lines.data();
1102 }
1103 }
1104
1105 // Keep intersect_line outside the loop, so it does not get reallocated.
1106 std::vector<std::pair<Intersection*, double>> intersect_line;
1107 for (size_t min_idx = 0; min_idx < intersections.size();) {
1108 intersect_line.clear();
1109 // All the nearest points (T-joints) ending at the same line are projected onto this line. Because of it, it can easily find the nearest point.
1110 {
1111 const Vec2d line_dir = intersections[min_idx].closest_line->vector().cast<double>();
1112 size_t max_idx = min_idx;
1113 for (; max_idx < intersections.size() &&
1114 intersections[min_idx].closest_line == intersections[max_idx].closest_line &&
1115 intersections[min_idx].left == intersections[max_idx].left;
1116 ++ max_idx)
1117 intersect_line.emplace_back(&intersections[max_idx], line_dir.dot(intersections[max_idx].intersect_point.cast<double>()));
1118 min_idx = max_idx;
1119 assert(intersect_line.size() > 0);
1120 // Sort the intersections along line_dir.
1121 std::sort(intersect_line.begin(), intersect_line.end(), [](const auto &i1, const auto &i2) { return i1.second < i2.second; });
1122 }
1123
1124 if (intersect_line.size() == 1) {
1125 // Simple case: The current intersection is the only one touching its adjacent line.
1126 Intersection &first_i = *intersect_line.front().first;
1127 update_merged_polyline(first_i);
1128 if (first_i.fresh()) {
1129 // Try to connect left or right. If not enough space for hook_length, take the longer side.
1130 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1131 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-add_hook0-pre-%d-%d.svg", iRun, iStep), { first_i.intersect_point });
1132 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1133 add_hook(first_i, scaled_offset, hook_length, scaled_trim_distance, rtree, lines_src);
1134 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1135 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-add_hook0-pre-%d-%d.svg", iRun, iStep), { first_i.intersect_point });
1136 ++ iStep;
1137 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1138 first_i.used = true;
1139 }
1140 continue;
1141 }
1142
1143 for (size_t first_idx = 0; first_idx < intersect_line.size(); ++ first_idx) {
1144 Intersection &first_i = *intersect_line[first_idx].first;
1145 update_merged_polyline(first_i);
1146 if (! first_i.fresh())
1147 // The intersection has been processed, or the polyline has been merged to another polyline.
1148 continue;
1149
1150 // Get the previous or next intersection on the same line, pick the closer one.
1151 if (first_idx > 0)
1152 update_merged_polyline(*intersect_line[first_idx - 1].first);
1153 if (first_idx + 1 < intersect_line.size())
1154 update_merged_polyline(*intersect_line[first_idx + 1].first);
1155 Intersection &nearest_i = *get_nearest_intersection(intersect_line, first_idx);
1156 assert(first_i.closest_line == nearest_i.closest_line);
1157 assert(first_i.intersect_line != nearest_i.intersect_line);
1158 assert(first_i.intersect_line != first_i.closest_line);
1159 assert(nearest_i.intersect_line != first_i.closest_line);
1160 // A line between two intersections points
1161 Line offset_line = create_offset_line(Line(first_i.intersect_point, nearest_i.intersect_point), first_i, scaled_offset);
1162 // Check if both intersections lie on the offset_line and simultaneously get their points of intersecting.
1163 // These points are used as start and end of the hook
1164 Point first_i_point, nearest_i_point;
1165 bool could_connect = false;
1166 if (nearest_i.fresh()) {
1167 could_connect = first_i.intersect_line->intersection(offset_line, &first_i_point) &&
1168 nearest_i.intersect_line->intersection(offset_line, &nearest_i_point);
1169 assert(could_connect);
1170 }
1171 Points &first_points = first_i.intersect_pl->points;
1172 Points &second_points = nearest_i.intersect_pl->points;
1173 could_connect &= (nearest_i_point - first_i_point).cast<double>().squaredNorm() <= Slic3r::sqr(hook_length_max);
1174 if (could_connect) {
1175 // Both intersections are so close that their polylines can be connected.
1176 // Verify that no other infill line intersects this anchor line.
1177 closest.clear();
1178 rtree.query(
1179 bgi::intersects(mk_rtree_seg(first_i_point, nearest_i_point)) &&
1180 bgi::satisfies([&first_i, &nearest_i, &lines_src](const auto &item)
1181 { return item.second != first_i.intersect_line - lines_src.data() && item.second != nearest_i.intersect_line - lines_src.data(); }),
1182 std::back_inserter(closest));
1183 could_connect = closest.empty();
1184 #if 0
1185 // Avoid self intersections. Maybe it is better to trim the self intersection after the connection?
1186 if (could_connect && first_i.intersect_pl != nearest_i.intersect_pl) {
1187 Line l(first_i_point, nearest_i_point);
1188 could_connect = ! first_i.other_hook_intersects(l) && ! nearest_i.other_hook_intersects(l);
1189 }
1190 #endif
1191 }
1192 bool connected = false;
1193 if (could_connect) {
1194 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1195 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-connecting-pre-%d-%d.svg", iRun, iStep), { first_i.intersect_point, nearest_i.intersect_point });
1196 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1197 // No other infill line intersects this anchor line. Extrude it as a whole.
1198 if (first_i.intersect_pl == nearest_i.intersect_pl) {
1199 // Both intersections are on the same polyline, that means a loop is being closed.
1200 assert(first_i.front != nearest_i.front);
1201 if (! first_i.front)
1202 std::swap(first_i_point, nearest_i_point);
1203 first_points.front() = first_i_point;
1204 first_points.back() = nearest_i_point;
1205 //FIXME trim the end of a closed loop a bit?
1206 first_points.emplace(first_points.begin(), nearest_i_point);
1207 } else {
1208 // Both intersections are on different polylines
1209 Line l(first_i_point, nearest_i_point);
1210 l.translate((perp(first_i.closest_line->vector().cast<double>().normalized()) * (first_i.left ? scaled_trim_distance : - scaled_trim_distance)).cast<coord_t>());
1211 Point pt_start, pt_end;
1212 bool trim_start = first_i .intersect_pl->points.size() == 3 && first_i .other_hook_intersects(l, pt_start);
1213 bool trim_end = nearest_i.intersect_pl->points.size() == 3 && nearest_i.other_hook_intersects(l, pt_end);
1214 first_points.reserve(first_points.size() + second_points.size());
1215 if (first_i.front)
1216 std::reverse(first_points.begin(), first_points.end());
1217 if (trim_start)
1218 first_points.front() = pt_start;
1219 first_points.back() = first_i_point;
1220 first_points.emplace_back(nearest_i_point);
1221 if (nearest_i.front)
1222 first_points.insert(first_points.end(), second_points.begin() + 1, second_points.end());
1223 else
1224 first_points.insert(first_points.end(), second_points.rbegin() + 1, second_points.rend());
1225 if (trim_end)
1226 first_points.back() = pt_end;
1227 // Keep the polyline at the lower index slot.
1228 if (first_i.intersect_pl < nearest_i.intersect_pl) {
1229 second_points.clear();
1230 merged_with[nearest_i.intersect_pl - lines.data()] = first_i.intersect_pl - lines.data();
1231 } else {
1232 second_points = std::move(first_points);
1233 first_points.clear();
1234 merged_with[first_i.intersect_pl - lines.data()] = nearest_i.intersect_pl - lines.data();
1235 }
1236 }
1237 nearest_i.used = true;
1238 connected = true;
1239 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1240 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-connecting-post-%d-%d.svg", iRun, iStep), { first_i.intersect_point, nearest_i.intersect_point });
1241 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1242 }
1243 if (! connected) {
1244 // Try to connect left or right. If not enough space for hook_length, take the longer side.
1245 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1246 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-add_hook-pre-%d-%d.svg", iRun, iStep), { first_i.intersect_point });
1247 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1248 add_hook(first_i, scaled_offset, hook_length, scaled_trim_distance, rtree, lines_src);
1249 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1250 export_infill_lines_to_svg(boundary, lines, debug_out_path("FillAdaptive-add_hook-post-%d-%d.svg", iRun, iStep), { first_i.intersect_point });
1251 #endif // ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1252 }
1253 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1254 ++ iStep;
1255 #endif ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1256 first_i.used = true;
1257 }
1258 }
1259
1260 Polylines polylines_out;
1261 polylines_out.reserve(polylines_out.size() + std::count_if(lines.begin(), lines.end(), [](const Polyline &pl) { return !pl.empty(); }));
1262 for (Polyline &pl : lines)
1263 if (!pl.empty()) polylines_out.emplace_back(std::move(pl));
1264 return polylines_out;
1265 }
1266
1267 #ifndef NDEBUG
has_no_collinear_lines(const Polylines & polylines)1268 bool has_no_collinear_lines(const Polylines &polylines)
1269 {
1270 // Create line end point lookup.
1271 struct LineEnd {
1272 LineEnd(const Polyline *line, bool start) : line(line), start(start) {}
1273 const Polyline *line;
1274 // Is it the start or end point?
1275 bool start;
1276 const Point& point() const { return start ? line->points.front() : line->points.back(); }
1277 const Point& other_point() const { return start ? line->points.back() : line->points.front(); }
1278 LineEnd other_end() const { return LineEnd(line, !start); }
1279 Vec2d vec() const { return Vec2d((this->other_point() - this->point()).cast<double>()); }
1280 bool operator==(const LineEnd &rhs) const { return this->line == rhs.line && this->start == rhs.start; }
1281 };
1282 struct LineEndAccessor {
1283 const Point* operator()(const LineEnd &pt) const { return &pt.point(); }
1284 };
1285 typedef ClosestPointInRadiusLookup<LineEnd, LineEndAccessor> ClosestPointLookupType;
1286 ClosestPointLookupType closest_end_point_lookup(coord_t(1001. * sqrt(2.)));
1287 for (const Polyline& pl : polylines) {
1288 // assert(pl.points.size() == 2);
1289 auto line_start = LineEnd(&pl, true);
1290 auto line_end = LineEnd(&pl, false);
1291
1292 auto assert_not_collinear = [&closest_end_point_lookup](const LineEnd &line_start) {
1293 std::vector<std::pair<const LineEnd*, double>> hits = closest_end_point_lookup.find_all(line_start.point());
1294 for (const std::pair<const LineEnd*, double> &hit : hits)
1295 if ((line_start.point() - hit.first->point()).cwiseAbs().maxCoeff() <= 1000) {
1296 // End points of the two lines are very close, they should have been merged together if they are collinear.
1297 Vec2d v1 = line_start.vec();
1298 Vec2d v2 = hit.first->vec();
1299 Vec2d v1n = v1.normalized();
1300 Vec2d v2n = v2.normalized();
1301 // The vectors must not be collinear.
1302 assert(std::abs(v1n.dot(v2n)) < cos(M_PI / 12.));
1303 }
1304 };
1305 assert_not_collinear(line_start);
1306 assert_not_collinear(line_end);
1307
1308 closest_end_point_lookup.insert(line_start);
1309 closest_end_point_lookup.insert(line_end);
1310 }
1311
1312 return true;
1313 }
1314 #endif
1315
_fill_surface_single(const FillParams & params,unsigned int thickness_layers,const std::pair<float,Point> & direction,ExPolygon expolygon,Polylines & polylines_out)1316 void Filler::_fill_surface_single(
1317 const FillParams ¶ms,
1318 unsigned int thickness_layers,
1319 const std::pair<float, Point> &direction,
1320 ExPolygon expolygon,
1321 Polylines &polylines_out)
1322 {
1323 assert (this->adapt_fill_octree);
1324
1325 Polylines all_polylines;
1326 {
1327 // 3 contexts for three directions of infill lines
1328 std::array<FillContext, 3> contexts {
1329 FillContext { *adapt_fill_octree, this->z, 0 },
1330 FillContext { *adapt_fill_octree, this->z, 1 },
1331 FillContext { *adapt_fill_octree, this->z, 2 }
1332 };
1333 // Generate the infill lines along the octree cells, merge touching lines of the same direction.
1334 size_t num_lines = 0;
1335 for (auto &context : contexts) {
1336 generate_infill_lines_recursive(context, adapt_fill_octree->root_cube, 0, int(adapt_fill_octree->cubes_properties.size()) - 1);
1337 num_lines += context.output_lines.size() + context.temp_lines.size();
1338 }
1339
1340 #if 0
1341 // Collect the lines, trim them by the expolygon.
1342 all_polylines.reserve(num_lines);
1343 auto boundary = to_polygons(expolygon);
1344 for (auto &context : contexts) {
1345 Polylines lines;
1346 lines.reserve(context.output_lines.size() + context.temp_lines.size());
1347 std::transform(context.output_lines.begin(), context.output_lines.end(), std::back_inserter(lines), [](const Line& l) { return Polyline{ l.a, l.b }; });
1348 for (const Line &l : context.temp_lines)
1349 if (l.a.x() != std::numeric_limits<coord_t>::max())
1350 lines.push_back({ l.a, l.b });
1351 // Crop all polylines
1352 append(all_polylines, intersection_pl(std::move(lines), boundary));
1353 }
1354 // assert(has_no_collinear_lines(all_polylines));
1355 #else
1356 // Collect the lines.
1357 std::vector<Line> lines;
1358 lines.reserve(num_lines);
1359 for (auto &context : contexts) {
1360 append(lines, context.output_lines);
1361 for (const Line &line : context.temp_lines)
1362 if (line.a.x() != std::numeric_limits<coord_t>::max())
1363 lines.emplace_back(line);
1364 }
1365 // Convert lines to polylines.
1366 all_polylines.reserve(lines.size());
1367 std::transform(lines.begin(), lines.end(), std::back_inserter(all_polylines), [](const Line& l) { return Polyline{ l.a, l.b }; });
1368 // Crop all polylines
1369 all_polylines = intersection_pl(std::move(all_polylines), to_polygons(expolygon));
1370 #endif
1371 }
1372
1373 // After intersection_pl some polylines with only one line are split into more lines
1374 for (Polyline &polyline : all_polylines) {
1375 //FIXME assert that all the points are collinear and in between the start and end point.
1376 if (polyline.points.size() > 2)
1377 polyline.points.erase(polyline.points.begin() + 1, polyline.points.end() - 1);
1378 }
1379 // assert(has_no_collinear_lines(all_polylines));
1380
1381 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1382 {
1383 static int iRun = 0;
1384 export_infill_lines_to_svg(expolygon, all_polylines, debug_out_path("FillAdaptive-initial-%d.svg", iRun++));
1385 }
1386 #endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
1387
1388 const auto hook_length = coordf_t(std::min<float>(std::numeric_limits<coord_t>::max(), scale_(params.anchor_length)));
1389 const auto hook_length_max = coordf_t(std::min<float>(std::numeric_limits<coord_t>::max(), scale_(params.anchor_length_max)));
1390
1391 Polylines all_polylines_with_hooks = all_polylines.size() > 1 ? connect_lines_using_hooks(std::move(all_polylines), expolygon, this->spacing, hook_length, hook_length_max) : std::move(all_polylines);
1392
1393 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1394 {
1395 static int iRun = 0;
1396 export_infill_lines_to_svg(expolygon, all_polylines_with_hooks, debug_out_path("FillAdaptive-hooks-%d.svg", iRun++));
1397 }
1398 #endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
1399
1400 if (params.dont_connect() || all_polylines_with_hooks.size() <= 1)
1401 append(polylines_out, chain_polylines(std::move(all_polylines_with_hooks)));
1402 else
1403 connect_infill(std::move(all_polylines_with_hooks), expolygon, polylines_out, this->spacing, params);
1404
1405 #ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
1406 {
1407 static int iRun = 0;
1408 export_infill_lines_to_svg(expolygon, polylines_out, debug_out_path("FillAdaptive-final-%d.svg", iRun ++));
1409 }
1410 #endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
1411 }
1412
bbox_max_radius(const BoundingBoxf3 & bbox,const Vec3d & center)1413 static double bbox_max_radius(const BoundingBoxf3 &bbox, const Vec3d ¢er)
1414 {
1415 const auto p = (bbox.min - center);
1416 const auto s = bbox.size();
1417 double r2max = 0.;
1418 for (int i = 0; i < 8; ++ i)
1419 r2max = std::max(r2max, (p + Vec3d(s.x() * double(i & 1), s.y() * double(i & 2), s.z() * double(i & 4))).squaredNorm());
1420 return sqrt(r2max);
1421 }
1422
make_cubes_properties(double max_cube_edge_length,double line_spacing)1423 static std::vector<CubeProperties> make_cubes_properties(double max_cube_edge_length, double line_spacing)
1424 {
1425 max_cube_edge_length += EPSILON;
1426
1427 std::vector<CubeProperties> cubes_properties;
1428 for (double edge_length = line_spacing * 2.;; edge_length *= 2.)
1429 {
1430 CubeProperties props{};
1431 props.edge_length = edge_length;
1432 props.height = edge_length * sqrt(3);
1433 props.diagonal_length = edge_length * sqrt(2);
1434 props.line_z_distance = edge_length / sqrt(3);
1435 props.line_xy_distance = edge_length / sqrt(6);
1436 cubes_properties.emplace_back(props);
1437 if (edge_length > max_cube_edge_length)
1438 break;
1439 }
1440 return cubes_properties;
1441 }
1442
is_overhang_triangle(const Vec3d & a,const Vec3d & b,const Vec3d & c,const Vec3d & up)1443 static inline bool is_overhang_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &up)
1444 {
1445 // Calculate triangle normal.
1446 auto n = (b - a).cross(c - b);
1447 return n.dot(up) > 0.707 * n.norm();
1448 }
1449
transform_center(Cube * current_cube,const Eigen::Matrix3d & rot)1450 static void transform_center(Cube *current_cube, const Eigen::Matrix3d &rot)
1451 {
1452 #ifndef NDEBUG
1453 current_cube->center_octree = current_cube->center;
1454 #endif // NDEBUG
1455 current_cube->center = rot * current_cube->center;
1456 for (auto *child : current_cube->children)
1457 if (child)
1458 transform_center(child, rot);
1459 }
1460
build_octree(const indexed_triangle_set & triangle_mesh,const std::vector<Vec3d> & overhang_triangles,coordf_t line_spacing,bool support_overhangs_only)1461 OctreePtr build_octree(
1462 // Mesh is rotated to the coordinate system of the octree.
1463 const indexed_triangle_set &triangle_mesh,
1464 // Overhang triangles extracted from fill surfaces with stInternalBridge type,
1465 // rotated to the coordinate system of the octree.
1466 const std::vector<Vec3d> &overhang_triangles,
1467 coordf_t line_spacing,
1468 bool support_overhangs_only)
1469 {
1470 assert(line_spacing > 0);
1471 assert(! std::isnan(line_spacing));
1472
1473 BoundingBox3Base<Vec3f> bbox(triangle_mesh.vertices);
1474 Vec3d cube_center = bbox.center().cast<double>();
1475 std::vector<CubeProperties> cubes_properties = make_cubes_properties(double(bbox.size().maxCoeff()), line_spacing);
1476 auto octree = OctreePtr(new Octree(cube_center, cubes_properties));
1477
1478 if (cubes_properties.size() > 1) {
1479 Octree *octree_ptr = octree.get();
1480 double edge_length_half = 0.5 * cubes_properties.back().edge_length;
1481 Vec3d diag_half(edge_length_half, edge_length_half, edge_length_half);
1482 int max_depth = int(cubes_properties.size()) - 1;
1483 auto process_triangle = [octree_ptr, max_depth, diag_half](const Vec3d &a, const Vec3d &b, const Vec3d &c) {
1484 octree_ptr->insert_triangle(
1485 a, b, c,
1486 octree_ptr->root_cube,
1487 BoundingBoxf3(octree_ptr->root_cube->center - diag_half, octree_ptr->root_cube->center + diag_half),
1488 max_depth);
1489 };
1490 auto up_vector = support_overhangs_only ? Vec3d(transform_to_octree() * Vec3d(0., 0., 1.)) : Vec3d();
1491 for (auto &tri : triangle_mesh.indices) {
1492 auto a = triangle_mesh.vertices[tri[0]].cast<double>();
1493 auto b = triangle_mesh.vertices[tri[1]].cast<double>();
1494 auto c = triangle_mesh.vertices[tri[2]].cast<double>();
1495 if (! support_overhangs_only || is_overhang_triangle(a, b, c, up_vector))
1496 process_triangle(a, b, c);
1497 }
1498 for (size_t i = 0; i < overhang_triangles.size(); i += 3)
1499 process_triangle(overhang_triangles[i], overhang_triangles[i + 1], overhang_triangles[i + 2]);
1500 {
1501 // Transform the octree to world coordinates to reduce computation when extracting infill lines.
1502 auto rot = transform_to_world().toRotationMatrix();
1503 transform_center(octree->root_cube, rot);
1504 octree->origin = rot * octree->origin;
1505 }
1506 }
1507
1508 return octree;
1509 }
1510
insert_triangle(const Vec3d & a,const Vec3d & b,const Vec3d & c,Cube * current_cube,const BoundingBoxf3 & current_bbox,int depth)1511 void Octree::insert_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, Cube *current_cube, const BoundingBoxf3 ¤t_bbox, int depth)
1512 {
1513 assert(current_cube);
1514 assert(depth > 0);
1515
1516 // Squared radius of a sphere around the child cube.
1517 const double r2_cube = Slic3r::sqr(0.5 * this->cubes_properties[-- depth].height + EPSILON);
1518
1519 for (size_t i = 0; i < 8; ++ i) {
1520 const Vec3d &child_center_dir = child_centers[i];
1521 // Calculate a slightly expanded bounding box of a child cube to cope with triangles touching a cube wall and other numeric errors.
1522 // We will rather densify the octree a bit more than necessary instead of missing a triangle.
1523 BoundingBoxf3 bbox;
1524 for (int k = 0; k < 3; ++ k) {
1525 if (child_center_dir[k] == -1.) {
1526 bbox.min[k] = current_bbox.min[k];
1527 bbox.max[k] = current_cube->center[k] + EPSILON;
1528 } else {
1529 bbox.min[k] = current_cube->center[k] - EPSILON;
1530 bbox.max[k] = current_bbox.max[k];
1531 }
1532 }
1533 Vec3d child_center = current_cube->center + (child_center_dir * (this->cubes_properties[depth].edge_length / 2.));
1534 //if (dist2_to_triangle(a, b, c, child_center) < r2_cube) {
1535 if (triangle_AABB_intersects(a, b, c, bbox)) {
1536 if (! current_cube->children[i])
1537 current_cube->children[i] = this->pool.construct(child_center);
1538 if (depth > 0)
1539 this->insert_triangle(a, b, c, current_cube->children[i], bbox, depth);
1540 }
1541 }
1542 }
1543
1544 } // namespace FillAdaptive
1545 } // namespace Slic3r
1546