1 #include "IndexedMesh.hpp"
2 #include "Concurrency.hpp"
3
4 #include <libslic3r/AABBTreeIndirect.hpp>
5 #include <libslic3r/TriangleMesh.hpp>
6
7 #include <numeric>
8
9 #ifdef SLIC3R_HOLE_RAYCASTER
10 #include <libslic3r/SLA/Hollowing.hpp>
11 #endif
12
13 namespace Slic3r { namespace sla {
14
15 class IndexedMesh::AABBImpl {
16 private:
17 AABBTreeIndirect::Tree3f m_tree;
18
19 public:
init(const TriangleMesh & tm)20 void init(const TriangleMesh& tm)
21 {
22 m_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
23 tm.its.vertices, tm.its.indices);
24 }
25
intersect_ray(const TriangleMesh & tm,const Vec3d & s,const Vec3d & dir,igl::Hit & hit)26 void intersect_ray(const TriangleMesh& tm,
27 const Vec3d& s, const Vec3d& dir, igl::Hit& hit)
28 {
29 AABBTreeIndirect::intersect_ray_first_hit(tm.its.vertices,
30 tm.its.indices,
31 m_tree,
32 s, dir, hit);
33 }
34
intersect_ray(const TriangleMesh & tm,const Vec3d & s,const Vec3d & dir,std::vector<igl::Hit> & hits)35 void intersect_ray(const TriangleMesh& tm,
36 const Vec3d& s, const Vec3d& dir, std::vector<igl::Hit>& hits)
37 {
38 AABBTreeIndirect::intersect_ray_all_hits(tm.its.vertices,
39 tm.its.indices,
40 m_tree,
41 s, dir, hits);
42 }
43
squared_distance(const TriangleMesh & tm,const Vec3d & point,int & i,Eigen::Matrix<double,1,3> & closest)44 double squared_distance(const TriangleMesh& tm,
45 const Vec3d& point, int& i, Eigen::Matrix<double, 1, 3>& closest) {
46 size_t idx_unsigned = 0;
47 Vec3d closest_vec3d(closest);
48 double dist = AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
49 tm.its.vertices,
50 tm.its.indices,
51 m_tree, point, idx_unsigned, closest_vec3d);
52 i = int(idx_unsigned);
53 closest = closest_vec3d;
54 return dist;
55 }
56 };
57
58 static const constexpr double MESH_EPS = 1e-6;
59
IndexedMesh(const TriangleMesh & tmesh)60 IndexedMesh::IndexedMesh(const TriangleMesh& tmesh)
61 : m_aabb(new AABBImpl()), m_tm(&tmesh)
62 {
63 auto&& bb = tmesh.bounding_box();
64 m_ground_level += bb.min(Z);
65
66 // Build the AABB accelaration tree
67 m_aabb->init(tmesh);
68 }
69
~IndexedMesh()70 IndexedMesh::~IndexedMesh() {}
71
IndexedMesh(const IndexedMesh & other)72 IndexedMesh::IndexedMesh(const IndexedMesh &other):
73 m_tm(other.m_tm), m_ground_level(other.m_ground_level),
74 m_aabb( new AABBImpl(*other.m_aabb) ) {}
75
76
operator =(const IndexedMesh & other)77 IndexedMesh &IndexedMesh::operator=(const IndexedMesh &other)
78 {
79 m_tm = other.m_tm;
80 m_ground_level = other.m_ground_level;
81 m_aabb.reset(new AABBImpl(*other.m_aabb)); return *this;
82 }
83
84 IndexedMesh &IndexedMesh::operator=(IndexedMesh &&other) = default;
85
86 IndexedMesh::IndexedMesh(IndexedMesh &&other) = default;
87
88
89
vertices() const90 const std::vector<Vec3f>& IndexedMesh::vertices() const
91 {
92 return m_tm->its.vertices;
93 }
94
95
96
indices() const97 const std::vector<Vec3i>& IndexedMesh::indices() const
98 {
99 return m_tm->its.indices;
100 }
101
102
103
vertices(size_t idx) const104 const Vec3f& IndexedMesh::vertices(size_t idx) const
105 {
106 return m_tm->its.vertices[idx];
107 }
108
109
110
indices(size_t idx) const111 const Vec3i& IndexedMesh::indices(size_t idx) const
112 {
113 return m_tm->its.indices[idx];
114 }
115
116
117
normal_by_face_id(int face_id) const118 Vec3d IndexedMesh::normal_by_face_id(int face_id) const {
119 return m_tm->stl.facet_start[face_id].normal.cast<double>();
120 }
121
122
123 IndexedMesh::hit_result
query_ray_hit(const Vec3d & s,const Vec3d & dir) const124 IndexedMesh::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
125 {
126 assert(is_approx(dir.norm(), 1.));
127 igl::Hit hit;
128 hit.t = std::numeric_limits<float>::infinity();
129
130 #ifdef SLIC3R_HOLE_RAYCASTER
131 if (! m_holes.empty()) {
132
133 // If there are holes, the hit_results will be made by
134 // query_ray_hits (object) and filter_hits (holes):
135 return filter_hits(query_ray_hits(s, dir));
136 }
137 #endif
138
139 m_aabb->intersect_ray(*m_tm, s, dir, hit);
140 hit_result ret(*this);
141 ret.m_t = double(hit.t);
142 ret.m_dir = dir;
143 ret.m_source = s;
144 if(!std::isinf(hit.t) && !std::isnan(hit.t)) {
145 ret.m_normal = this->normal_by_face_id(hit.id);
146 ret.m_face_id = hit.id;
147 }
148
149 return ret;
150 }
151
152 std::vector<IndexedMesh::hit_result>
query_ray_hits(const Vec3d & s,const Vec3d & dir) const153 IndexedMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const
154 {
155 std::vector<IndexedMesh::hit_result> outs;
156 std::vector<igl::Hit> hits;
157 m_aabb->intersect_ray(*m_tm, s, dir, hits);
158
159 // The sort is necessary, the hits are not always sorted.
160 std::sort(hits.begin(), hits.end(),
161 [](const igl::Hit& a, const igl::Hit& b) { return a.t < b.t; });
162
163 // Remove duplicates. They sometimes appear, for example when the ray is cast
164 // along an axis of a cube due to floating-point approximations in igl (?)
165 hits.erase(std::unique(hits.begin(), hits.end(),
166 [](const igl::Hit& a, const igl::Hit& b)
167 { return a.t == b.t; }),
168 hits.end());
169
170 // Convert the igl::Hit into hit_result
171 outs.reserve(hits.size());
172 for (const igl::Hit& hit : hits) {
173 outs.emplace_back(IndexedMesh::hit_result(*this));
174 outs.back().m_t = double(hit.t);
175 outs.back().m_dir = dir;
176 outs.back().m_source = s;
177 if(!std::isinf(hit.t) && !std::isnan(hit.t)) {
178 outs.back().m_normal = this->normal_by_face_id(hit.id);
179 outs.back().m_face_id = hit.id;
180 }
181 }
182
183 return outs;
184 }
185
186
187 #ifdef SLIC3R_HOLE_RAYCASTER
filter_hits(const std::vector<IndexedMesh::hit_result> & object_hits) const188 IndexedMesh::hit_result IndexedMesh::filter_hits(
189 const std::vector<IndexedMesh::hit_result>& object_hits) const
190 {
191 assert(! m_holes.empty());
192 hit_result out(*this);
193
194 if (object_hits.empty())
195 return out;
196
197 const Vec3d& s = object_hits.front().source();
198 const Vec3d& dir = object_hits.front().direction();
199
200 // A helper struct to save an intersetion with a hole
201 struct HoleHit {
202 HoleHit(float t_p, const Vec3d& normal_p, bool entry_p) :
203 t(t_p), normal(normal_p), entry(entry_p) {}
204 float t;
205 Vec3d normal;
206 bool entry;
207 };
208 std::vector<HoleHit> hole_isects;
209 hole_isects.reserve(m_holes.size());
210
211 auto sf = s.cast<float>();
212 auto dirf = dir.cast<float>();
213
214 // Collect hits on all holes, preserve information about entry/exit
215 for (const sla::DrainHole& hole : m_holes) {
216 std::array<std::pair<float, Vec3d>, 2> isects;
217 if (hole.get_intersections(sf, dirf, isects)) {
218 // Ignore hole hits behind the source
219 if (isects[0].first > 0.f) hole_isects.emplace_back(isects[0].first, isects[0].second, true);
220 if (isects[1].first > 0.f) hole_isects.emplace_back(isects[1].first, isects[1].second, false);
221 }
222 }
223
224 // Holes can intersect each other, sort the hits by t
225 std::sort(hole_isects.begin(), hole_isects.end(),
226 [](const HoleHit& a, const HoleHit& b) { return a.t < b.t; });
227
228 // Now inspect the intersections with object and holes, in the order of
229 // increasing distance. Keep track how deep are we nested in mesh/holes and
230 // pick the correct intersection.
231 // This needs to be done twice - first to find out how deep in the structure
232 // the source is, then to pick the correct intersection.
233 int hole_nested = 0;
234 int object_nested = 0;
235 for (int dry_run=1; dry_run>=0; --dry_run) {
236 hole_nested = -hole_nested;
237 object_nested = -object_nested;
238
239 bool is_hole = false;
240 bool is_entry = false;
241 const HoleHit* next_hole_hit = hole_isects.empty() ? nullptr : &hole_isects.front();
242 const hit_result* next_mesh_hit = &object_hits.front();
243
244 while (next_hole_hit || next_mesh_hit) {
245 if (next_hole_hit && next_mesh_hit) // still have hole and obj hits
246 is_hole = (next_hole_hit->t < next_mesh_hit->m_t);
247 else
248 is_hole = next_hole_hit; // one or the other ran out
249
250 // Is this entry or exit hit?
251 is_entry = is_hole ? next_hole_hit->entry : ! next_mesh_hit->is_inside();
252
253 if (! dry_run) {
254 if (! is_hole && hole_nested == 0) {
255 // This is a valid object hit
256 return *next_mesh_hit;
257 }
258 if (is_hole && ! is_entry && object_nested != 0) {
259 // This holehit is the one we seek
260 out.m_t = next_hole_hit->t;
261 out.m_normal = next_hole_hit->normal;
262 out.m_source = s;
263 out.m_dir = dir;
264 return out;
265 }
266 }
267
268 // Increase/decrease the counter
269 (is_hole ? hole_nested : object_nested) += (is_entry ? 1 : -1);
270
271 // Advance the respective pointer
272 if (is_hole && next_hole_hit++ == &hole_isects.back())
273 next_hole_hit = nullptr;
274 if (! is_hole && next_mesh_hit++ == &object_hits.back())
275 next_mesh_hit = nullptr;
276 }
277 }
278
279 // if we got here, the ray ended up in infinity
280 return out;
281 }
282 #endif
283
284
squared_distance(const Vec3d & p,int & i,Vec3d & c) const285 double IndexedMesh::squared_distance(const Vec3d &p, int& i, Vec3d& c) const {
286 double sqdst = 0;
287 Eigen::Matrix<double, 1, 3> pp = p;
288 Eigen::Matrix<double, 1, 3> cc;
289 sqdst = m_aabb->squared_distance(*m_tm, pp, i, cc);
290 c = cc;
291 return sqdst;
292 }
293
294
point_on_edge(const Vec3d & p,const Vec3d & e1,const Vec3d & e2,double eps=0.05)295 static bool point_on_edge(const Vec3d& p, const Vec3d& e1, const Vec3d& e2,
296 double eps = 0.05)
297 {
298 using Line3D = Eigen::ParametrizedLine<double, 3>;
299
300 auto line = Line3D::Through(e1, e2);
301 double d = line.distance(p);
302 return std::abs(d) < eps;
303 }
304
normals(const PointSet & points,const IndexedMesh & mesh,double eps,std::function<void ()> thr,const std::vector<unsigned> & pt_indices)305 PointSet normals(const PointSet& points,
306 const IndexedMesh& mesh,
307 double eps,
308 std::function<void()> thr, // throw on cancel
309 const std::vector<unsigned>& pt_indices)
310 {
311 if (points.rows() == 0 || mesh.vertices().empty() || mesh.indices().empty())
312 return {};
313
314 std::vector<unsigned> range = pt_indices;
315 if (range.empty()) {
316 range.resize(size_t(points.rows()), 0);
317 std::iota(range.begin(), range.end(), 0);
318 }
319
320 PointSet ret(range.size(), 3);
321
322 // for (size_t ridx = 0; ridx < range.size(); ++ridx)
323 ccr::for_each(size_t(0), range.size(),
324 [&ret, &mesh, &points, thr, eps, &range](size_t ridx) {
325 thr();
326 unsigned el = range[ridx];
327 auto eidx = Eigen::Index(el);
328 int faceid = 0;
329 Vec3d p;
330
331 mesh.squared_distance(points.row(eidx), faceid, p);
332
333 auto trindex = mesh.indices(faceid);
334
335 const Vec3d &p1 = mesh.vertices(trindex(0)).cast<double>();
336 const Vec3d &p2 = mesh.vertices(trindex(1)).cast<double>();
337 const Vec3d &p3 = mesh.vertices(trindex(2)).cast<double>();
338
339 // We should check if the point lies on an edge of the hosting
340 // triangle. If it does then all the other triangles using the
341 // same two points have to be searched and the final normal should
342 // be some kind of aggregation of the participating triangle
343 // normals. We should also consider the cases where the support
344 // point lies right on a vertex of its triangle. The procedure is
345 // the same, get the neighbor triangles and calculate an average
346 // normal.
347
348 // mark the vertex indices of the edge. ia and ib marks and edge
349 // ic will mark a single vertex.
350 int ia = -1, ib = -1, ic = -1;
351
352 if (std::abs((p - p1).norm()) < eps) {
353 ic = trindex(0);
354 } else if (std::abs((p - p2).norm()) < eps) {
355 ic = trindex(1);
356 } else if (std::abs((p - p3).norm()) < eps) {
357 ic = trindex(2);
358 } else if (point_on_edge(p, p1, p2, eps)) {
359 ia = trindex(0);
360 ib = trindex(1);
361 } else if (point_on_edge(p, p2, p3, eps)) {
362 ia = trindex(1);
363 ib = trindex(2);
364 } else if (point_on_edge(p, p1, p3, eps)) {
365 ia = trindex(0);
366 ib = trindex(2);
367 }
368
369 // vector for the neigboring triangles including the detected one.
370 std::vector<size_t> neigh;
371 if (ic >= 0) { // The point is right on a vertex of the triangle
372 for (size_t n = 0; n < mesh.indices().size(); ++n) {
373 thr();
374 Vec3i ni = mesh.indices(n);
375 if ((ni(X) == ic || ni(Y) == ic || ni(Z) == ic))
376 neigh.emplace_back(n);
377 }
378 } else if (ia >= 0 && ib >= 0) { // the point is on and edge
379 // now get all the neigboring triangles
380 for (size_t n = 0; n < mesh.indices().size(); ++n) {
381 thr();
382 Vec3i ni = mesh.indices(n);
383 if ((ni(X) == ia || ni(Y) == ia || ni(Z) == ia) &&
384 (ni(X) == ib || ni(Y) == ib || ni(Z) == ib))
385 neigh.emplace_back(n);
386 }
387 }
388
389 // Calculate the normals for the neighboring triangles
390 std::vector<Vec3d> neighnorms;
391 neighnorms.reserve(neigh.size());
392 for (size_t &tri_id : neigh)
393 neighnorms.emplace_back(mesh.normal_by_face_id(tri_id));
394
395 // Throw out duplicates. They would cause trouble with summing. We
396 // will use std::unique which works on sorted ranges. We will sort
397 // by the coefficient-wise sum of the normals. It should force the
398 // same elements to be consecutive.
399 std::sort(neighnorms.begin(), neighnorms.end(),
400 [](const Vec3d &v1, const Vec3d &v2) {
401 return v1.sum() < v2.sum();
402 });
403
404 auto lend = std::unique(neighnorms.begin(), neighnorms.end(),
405 [](const Vec3d &n1, const Vec3d &n2) {
406 // Compare normals for equivalence.
407 // This is controvers stuff.
408 auto deq = [](double a, double b) {
409 return std::abs(a - b) < 1e-3;
410 };
411 return deq(n1(X), n2(X)) &&
412 deq(n1(Y), n2(Y)) &&
413 deq(n1(Z), n2(Z));
414 });
415
416 if (!neighnorms.empty()) { // there were neighbors to count with
417 // sum up the normals and then normalize the result again.
418 // This unification seems to be enough.
419 Vec3d sumnorm(0, 0, 0);
420 sumnorm = std::accumulate(neighnorms.begin(), lend, sumnorm);
421 sumnorm.normalize();
422 ret.row(long(ridx)) = sumnorm;
423 } else { // point lies safely within its triangle
424 Eigen::Vector3d U = p2 - p1;
425 Eigen::Vector3d V = p3 - p1;
426 ret.row(long(ridx)) = U.cross(V).normalized();
427 }
428 });
429
430 return ret;
431 }
432
433 }} // namespace Slic3r::sla
434