1 /*
2 * This program is free software; you can redistribute it and/or
3 * modify it under the terms of the GNU General Public License
4 * as published by the Free Software Foundation; either version 2
5 * of the License, or (at your option) any later version.
6 *
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
11 *
12 * You should have received a copy of the GNU General Public License
13 * along with this program; if not, write to the Free Software Foundation,
14 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15 */
16
17 /** \file
18 * \ingroup bli
19 */
20
21 /* The #blender::meshintersect API needs GMP. */
22 #ifdef WITH_GMP
23
24 # include <algorithm>
25 # include <fstream>
26 # include <iostream>
27
28 # include "BLI_allocator.hh"
29 # include "BLI_array.hh"
30 # include "BLI_assert.h"
31 # include "BLI_delaunay_2d.h"
32 # include "BLI_double3.hh"
33 # include "BLI_float3.hh"
34 # include "BLI_hash.hh"
35 # include "BLI_kdopbvh.h"
36 # include "BLI_map.hh"
37 # include "BLI_math_boolean.hh"
38 # include "BLI_math_mpq.hh"
39 # include "BLI_mpq2.hh"
40 # include "BLI_mpq3.hh"
41 # include "BLI_span.hh"
42 # include "BLI_task.h"
43 # include "BLI_threads.h"
44 # include "BLI_vector.hh"
45 # include "BLI_vector_set.hh"
46
47 # include "PIL_time.h"
48
49 # include "BLI_mesh_intersect.hh"
50
51 // # define PERFDEBUG
52
53 namespace blender::meshintersect {
54
55 # ifdef PERFDEBUG
56 static void perfdata_init(void);
57 static void incperfcount(int countnum);
58 static void bumpperfcount(int countnum, int amt);
59 static void doperfmax(int maxnum, int val);
60 static void dump_perfdata(void);
61 # endif
62
63 /** For debugging, can disable threading in intersect code with this static constant. */
64 static constexpr bool intersect_use_threading = true;
65
Vert(const mpq3 & mco,const double3 & dco,int id,int orig)66 Vert::Vert(const mpq3 &mco, const double3 &dco, int id, int orig)
67 : co_exact(mco), co(dco), id(id), orig(orig)
68 {
69 }
70
operator ==(const Vert & other) const71 bool Vert::operator==(const Vert &other) const
72 {
73 return this->co_exact == other.co_exact;
74 }
75
hash() const76 uint64_t Vert::hash() const
77 {
78 return co_exact.hash();
79 }
80
operator <<(std::ostream & os,const Vert * v)81 std::ostream &operator<<(std::ostream &os, const Vert *v)
82 {
83 constexpr int dbg_level = 0;
84 os << "v" << v->id;
85 if (v->orig != NO_INDEX) {
86 os << "o" << v->orig;
87 }
88 os << v->co;
89 if (dbg_level > 0) {
90 os << "=" << v->co_exact;
91 }
92 return os;
93 }
94
operator ==(const Plane & other) const95 bool Plane::operator==(const Plane &other) const
96 {
97 return norm_exact == other.norm_exact && d_exact == other.d_exact;
98 }
99
make_canonical()100 void Plane::make_canonical()
101 {
102 if (norm_exact[0] != 0) {
103 mpq_class den = norm_exact[0];
104 norm_exact = mpq3(1, norm_exact[1] / den, norm_exact[2] / den);
105 d_exact = d_exact / den;
106 }
107 else if (norm_exact[1] != 0) {
108 mpq_class den = norm_exact[1];
109 norm_exact = mpq3(0, 1, norm_exact[2] / den);
110 d_exact = d_exact / den;
111 }
112 else {
113 if (norm_exact[2] != 0) {
114 mpq_class den = norm_exact[2];
115 norm_exact = mpq3(0, 0, 1);
116 d_exact = d_exact / den;
117 }
118 else {
119 /* A degenerate plane. */
120 d_exact = 0;
121 }
122 }
123 norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d());
124 d = d_exact.get_d();
125 }
126
Plane(const mpq3 & norm_exact,const mpq_class & d_exact)127 Plane::Plane(const mpq3 &norm_exact, const mpq_class &d_exact)
128 : norm_exact(norm_exact), d_exact(d_exact)
129 {
130 norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d());
131 d = d_exact.get_d();
132 }
133
Plane(const double3 & norm,const double d)134 Plane::Plane(const double3 &norm, const double d) : norm(norm), d(d)
135 {
136 norm_exact = mpq3(0, 0, 0); /* Marks as "exact not yet populated". */
137 }
138
139 /** This is wrong for degenerate planes, but we don't expect to call it on those. */
exact_populated() const140 bool Plane::exact_populated() const
141 {
142 return norm_exact[0] != 0 || norm_exact[1] != 0 || norm_exact[2] != 0;
143 }
144
hash() const145 uint64_t Plane::hash() const
146 {
147 constexpr uint64_t h1 = 33;
148 constexpr uint64_t h2 = 37;
149 constexpr uint64_t h3 = 39;
150 uint64_t hashx = hash_mpq_class(this->norm_exact.x);
151 uint64_t hashy = hash_mpq_class(this->norm_exact.y);
152 uint64_t hashz = hash_mpq_class(this->norm_exact.z);
153 uint64_t hashd = hash_mpq_class(this->d_exact);
154 uint64_t ans = hashx ^ (hashy * h1) ^ (hashz * h1 * h2) ^ (hashd * h1 * h2 * h3);
155 return ans;
156 }
157
operator <<(std::ostream & os,const Plane * plane)158 std::ostream &operator<<(std::ostream &os, const Plane *plane)
159 {
160 os << "[" << plane->norm << ";" << plane->d << "]";
161 return os;
162 }
163
Face(Span<const Vert * > verts,int id,int orig,Span<int> edge_origs,Span<bool> is_intersect)164 Face::Face(
165 Span<const Vert *> verts, int id, int orig, Span<int> edge_origs, Span<bool> is_intersect)
166 : vert(verts), edge_orig(edge_origs), is_intersect(is_intersect), id(id), orig(orig)
167 {
168 }
169
Face(Span<const Vert * > verts,int id,int orig)170 Face::Face(Span<const Vert *> verts, int id, int orig) : vert(verts), id(id), orig(orig)
171 {
172 }
173
populate_plane(bool need_exact)174 void Face::populate_plane(bool need_exact)
175 {
176 if (plane != nullptr) {
177 if (!need_exact || plane->exact_populated()) {
178 return;
179 }
180 }
181 if (need_exact) {
182 mpq3 normal_exact;
183 if (vert.size() > 3) {
184 Array<mpq3> co(vert.size());
185 for (int i : index_range()) {
186 co[i] = vert[i]->co_exact;
187 }
188 normal_exact = mpq3::cross_poly(co);
189 }
190 else {
191 mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact;
192 mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact;
193 normal_exact = mpq3::cross(tr02, tr12);
194 }
195 mpq_class d_exact = -mpq3::dot(normal_exact, vert[0]->co_exact);
196 plane = new Plane(normal_exact, d_exact);
197 }
198 else {
199 double3 normal;
200 if (vert.size() > 3) {
201 Array<double3> co(vert.size());
202 for (int i : index_range()) {
203 co[i] = vert[i]->co;
204 }
205 normal = double3::cross_poly(co);
206 }
207 else {
208 double3 tr02 = vert[0]->co - vert[2]->co;
209 double3 tr12 = vert[1]->co - vert[2]->co;
210 normal = double3::cross_high_precision(tr02, tr12);
211 }
212 double d = -double3::dot(normal, vert[0]->co);
213 plane = new Plane(normal, d);
214 }
215 }
216
~Face()217 Face::~Face()
218 {
219 delete plane;
220 }
221
operator ==(const Face & other) const222 bool Face::operator==(const Face &other) const
223 {
224 if (this->size() != other.size()) {
225 return false;
226 }
227 for (FacePos i : index_range()) {
228 /* Can test pointer equality since we will have
229 * unique vert pointers for unique co_equal's. */
230 if (this->vert[i] != other.vert[i]) {
231 return false;
232 }
233 }
234 return true;
235 }
236
cyclic_equal(const Face & other) const237 bool Face::cyclic_equal(const Face &other) const
238 {
239 if (this->size() != other.size()) {
240 return false;
241 }
242 int flen = this->size();
243 for (FacePos start : index_range()) {
244 for (FacePos start_other : index_range()) {
245 bool ok = true;
246 for (int i = 0; ok && i < flen; ++i) {
247 FacePos p = (start + i) % flen;
248 FacePos p_other = (start_other + i) % flen;
249 if (this->vert[p] != other.vert[p_other]) {
250 ok = false;
251 }
252 }
253 if (ok) {
254 return true;
255 }
256 }
257 }
258 return false;
259 }
260
operator <<(std::ostream & os,const Face * f)261 std::ostream &operator<<(std::ostream &os, const Face *f)
262 {
263 os << "f" << f->id << "o" << f->orig << "[";
264 for (const Vert *v : *f) {
265 os << v;
266 if (v != f->vert[f->size() - 1]) {
267 os << " ";
268 }
269 }
270 os << "]";
271 if (f->orig != NO_INDEX) {
272 os << "o" << f->orig;
273 }
274 os << " e_orig[";
275 for (int i : f->index_range()) {
276 os << f->edge_orig[i];
277 if (f->is_intersect[i]) {
278 os << "#";
279 }
280 if (i != f->size() - 1) {
281 os << " ";
282 }
283 }
284 os << "]";
285 return os;
286 }
287
288 /**
289 * Un-comment the following to try using a spin-lock instead of
290 * a mutex in the arena allocation routines.
291 * Initial tests showed that it doesn't seem to help very much,
292 * if at all, to use a spin-lock.
293 */
294 // #define USE_SPINLOCK
295
296 /**
297 * #IMeshArena is the owner of the Vert and Face resources used
298 * during a run of one of the mesh-intersect main functions.
299 * It also keeps has a hash table of all Verts created so that it can
300 * ensure that only one instance of a Vert with a given co_exact will
301 * exist. I.e., it de-duplicates the vertices.
302 */
303 class IMeshArena::IMeshArenaImpl : NonCopyable, NonMovable {
304
305 /**
306 * Don't use Vert itself as key since resizing may move
307 * pointers to the Vert around, and we need to have those pointers
308 * stay the same throughout the lifetime of the #IMeshArena.
309 */
310 struct VSetKey {
311 Vert *vert;
312
VSetKeyblender::meshintersect::IMeshArena::IMeshArenaImpl::VSetKey313 VSetKey(Vert *p) : vert(p)
314 {
315 }
316
hashblender::meshintersect::IMeshArena::IMeshArenaImpl::VSetKey317 uint32_t hash() const
318 {
319 return vert->hash();
320 }
321
operator ==blender::meshintersect::IMeshArena::IMeshArenaImpl::VSetKey322 bool operator==(const VSetKey &other) const
323 {
324 return *this->vert == *other.vert;
325 }
326 };
327
328 VectorSet<VSetKey> vset_; /* TODO: replace with Set */
329
330 /**
331 * Ownership of the Vert memory is here, so destroying this reclaims that memory.
332 *
333 * TODO: replace these with pooled allocation, and just destroy the pools at the end.
334 */
335 Vector<std::unique_ptr<Vert>> allocated_verts_;
336 Vector<std::unique_ptr<Face>> allocated_faces_;
337
338 /* Use these to allocate ids when Verts and Faces are allocated. */
339 int next_vert_id_ = 0;
340 int next_face_id_ = 0;
341
342 /* Need a lock when multi-threading to protect allocation of new elements. */
343 # ifdef USE_SPINLOCK
344 SpinLock lock_;
345 # else
346 ThreadMutex *mutex_;
347 # endif
348
349 public:
IMeshArenaImpl()350 IMeshArenaImpl()
351 {
352 if (intersect_use_threading) {
353 # ifdef USE_SPINLOCK
354 BLI_spin_init(&lock_);
355 # else
356 mutex_ = BLI_mutex_alloc();
357 # endif
358 }
359 }
~IMeshArenaImpl()360 ~IMeshArenaImpl()
361 {
362 if (intersect_use_threading) {
363 # ifdef USE_SPINLOCK
364 BLI_spin_end(&lock_);
365 # else
366 BLI_mutex_free(mutex_);
367 # endif
368 }
369 }
370
reserve(int vert_num_hint,int face_num_hint)371 void reserve(int vert_num_hint, int face_num_hint)
372 {
373 vset_.reserve(vert_num_hint);
374 allocated_verts_.reserve(vert_num_hint);
375 allocated_faces_.reserve(face_num_hint);
376 }
377
tot_allocated_verts() const378 int tot_allocated_verts() const
379 {
380 return allocated_verts_.size();
381 }
382
tot_allocated_faces() const383 int tot_allocated_faces() const
384 {
385 return allocated_faces_.size();
386 }
387
add_or_find_vert(const mpq3 & co,int orig)388 const Vert *add_or_find_vert(const mpq3 &co, int orig)
389 {
390 double3 dco(co[0].get_d(), co[1].get_d(), co[2].get_d());
391 return add_or_find_vert(co, dco, orig);
392 }
393
add_or_find_vert(const double3 & co,int orig)394 const Vert *add_or_find_vert(const double3 &co, int orig)
395 {
396 mpq3 mco(co[0], co[1], co[2]);
397 return add_or_find_vert(mco, co, orig);
398 }
399
add_face(Span<const Vert * > verts,int orig,Span<int> edge_origs,Span<bool> is_intersect)400 Face *add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs, Span<bool> is_intersect)
401 {
402 Face *f = new Face(verts, next_face_id_++, orig, edge_origs, is_intersect);
403 if (intersect_use_threading) {
404 # ifdef USE_SPINLOCK
405 BLI_spin_lock(&lock_);
406 # else
407 BLI_mutex_lock(mutex_);
408 # endif
409 }
410 allocated_faces_.append(std::unique_ptr<Face>(f));
411 if (intersect_use_threading) {
412 # ifdef USE_SPINLOCK
413 BLI_spin_unlock(&lock_);
414 # else
415 BLI_mutex_unlock(mutex_);
416 # endif
417 }
418 return f;
419 }
420
add_face(Span<const Vert * > verts,int orig,Span<int> edge_origs)421 Face *add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs)
422 {
423 Array<bool> is_intersect(verts.size(), false);
424 return add_face(verts, orig, edge_origs, is_intersect);
425 }
426
add_face(Span<const Vert * > verts,int orig)427 Face *add_face(Span<const Vert *> verts, int orig)
428 {
429 Array<int> edge_origs(verts.size(), NO_INDEX);
430 Array<bool> is_intersect(verts.size(), false);
431 return add_face(verts, orig, edge_origs, is_intersect);
432 }
433
find_vert(const mpq3 & co)434 const Vert *find_vert(const mpq3 &co)
435 {
436 const Vert *ans;
437 Vert vtry(co, double3(), NO_INDEX, NO_INDEX);
438 VSetKey vskey(&vtry);
439 if (intersect_use_threading) {
440 # ifdef USE_SPINLOCK
441 BLI_spin_lock(&lock_);
442 # else
443 BLI_mutex_lock(mutex_);
444 # endif
445 }
446 int i = vset_.index_of_try(vskey);
447 if (i == -1) {
448 ans = nullptr;
449 }
450 else {
451 ans = vset_[i].vert;
452 }
453 if (intersect_use_threading) {
454 # ifdef USE_SPINLOCK
455 BLI_spin_unlock(&lock_);
456 # else
457 BLI_mutex_unlock(mutex_);
458 # endif
459 }
460 return ans;
461 }
462
463 /**
464 * This is slow. Only used for unit tests right now.
465 * Since it is only used for that purpose, access is not lock-protected.
466 * The argument vs can be a cyclic shift of the actual stored Face.
467 */
find_face(Span<const Vert * > vs)468 const Face *find_face(Span<const Vert *> vs)
469 {
470 Array<int> eorig(vs.size(), NO_INDEX);
471 Array<bool> is_intersect(vs.size(), false);
472 Face ftry(vs, NO_INDEX, NO_INDEX, eorig, is_intersect);
473 for (const int i : allocated_faces_.index_range()) {
474 if (ftry.cyclic_equal(*allocated_faces_[i])) {
475 return allocated_faces_[i].get();
476 }
477 }
478 return nullptr;
479 }
480
481 private:
add_or_find_vert(const mpq3 & mco,const double3 & dco,int orig)482 const Vert *add_or_find_vert(const mpq3 &mco, const double3 &dco, int orig)
483 {
484 /* Don't allocate Vert yet, in case it is already there. */
485 Vert vtry(mco, dco, NO_INDEX, NO_INDEX);
486 const Vert *ans;
487 VSetKey vskey(&vtry);
488 if (intersect_use_threading) {
489 # ifdef USE_SPINLOCK
490 BLI_spin_lock(&lock_);
491 # else
492 BLI_mutex_lock(mutex_);
493 # endif
494 }
495 int i = vset_.index_of_try(vskey);
496 if (i == -1) {
497 vskey.vert = new Vert(mco, dco, next_vert_id_++, orig);
498 vset_.add_new(vskey);
499 allocated_verts_.append(std::unique_ptr<Vert>(vskey.vert));
500 ans = vskey.vert;
501 }
502 else {
503 /* It was a duplicate, so return the existing one.
504 * Note that the returned Vert may have a different orig.
505 * This is the intended semantics: if the Vert already
506 * exists then we are merging verts and using the first-seen
507 * one as the canonical one. */
508 ans = vset_[i].vert;
509 }
510 if (intersect_use_threading) {
511 # ifdef USE_SPINLOCK
512 BLI_spin_unlock(&lock_);
513 # else
514 BLI_mutex_unlock(mutex_);
515 # endif
516 }
517 return ans;
518 };
519 };
520
IMeshArena()521 IMeshArena::IMeshArena()
522 {
523 pimpl_ = std::unique_ptr<IMeshArenaImpl>(new IMeshArenaImpl());
524 }
525
~IMeshArena()526 IMeshArena::~IMeshArena()
527 {
528 }
529
reserve(int vert_num_hint,int face_num_hint)530 void IMeshArena::reserve(int vert_num_hint, int face_num_hint)
531 {
532 pimpl_->reserve(vert_num_hint, face_num_hint);
533 }
534
tot_allocated_verts() const535 int IMeshArena::tot_allocated_verts() const
536 {
537 return pimpl_->tot_allocated_verts();
538 }
539
tot_allocated_faces() const540 int IMeshArena::tot_allocated_faces() const
541 {
542 return pimpl_->tot_allocated_faces();
543 }
544
add_or_find_vert(const mpq3 & co,int orig)545 const Vert *IMeshArena::add_or_find_vert(const mpq3 &co, int orig)
546 {
547 return pimpl_->add_or_find_vert(co, orig);
548 }
549
add_face(Span<const Vert * > verts,int orig,Span<int> edge_origs,Span<bool> is_intersect)550 Face *IMeshArena::add_face(Span<const Vert *> verts,
551 int orig,
552 Span<int> edge_origs,
553 Span<bool> is_intersect)
554 {
555 return pimpl_->add_face(verts, orig, edge_origs, is_intersect);
556 }
557
add_face(Span<const Vert * > verts,int orig,Span<int> edge_origs)558 Face *IMeshArena::add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs)
559 {
560 return pimpl_->add_face(verts, orig, edge_origs);
561 }
562
add_face(Span<const Vert * > verts,int orig)563 Face *IMeshArena::add_face(Span<const Vert *> verts, int orig)
564 {
565 return pimpl_->add_face(verts, orig);
566 }
567
add_or_find_vert(const double3 & co,int orig)568 const Vert *IMeshArena::add_or_find_vert(const double3 &co, int orig)
569 {
570 return pimpl_->add_or_find_vert(co, orig);
571 }
572
find_vert(const mpq3 & co) const573 const Vert *IMeshArena::find_vert(const mpq3 &co) const
574 {
575 return pimpl_->find_vert(co);
576 }
577
find_face(Span<const Vert * > verts) const578 const Face *IMeshArena::find_face(Span<const Vert *> verts) const
579 {
580 return pimpl_->find_face(verts);
581 }
582
set_faces(Span<Face * > faces)583 void IMesh::set_faces(Span<Face *> faces)
584 {
585 face_ = faces;
586 }
587
lookup_vert(const Vert * v) const588 int IMesh::lookup_vert(const Vert *v) const
589 {
590 BLI_assert(vert_populated_);
591 return vert_to_index_.lookup_default(v, NO_INDEX);
592 }
593
populate_vert()594 void IMesh::populate_vert()
595 {
596 /* This is likely an overestimate, since verts are shared between
597 * faces. It is ok if estimate is over or even under. */
598 constexpr int ESTIMATE_VERTS_PER_FACE = 4;
599 int estimate_num_verts = ESTIMATE_VERTS_PER_FACE * face_.size();
600 populate_vert(estimate_num_verts);
601 }
602
populate_vert(int max_verts)603 void IMesh::populate_vert(int max_verts)
604 {
605 if (vert_populated_) {
606 return;
607 }
608 vert_to_index_.reserve(max_verts);
609 int next_allocate_index = 0;
610 for (const Face *f : face_) {
611 for (const Vert *v : *f) {
612 if (v->id == 1) {
613 }
614 int index = vert_to_index_.lookup_default(v, NO_INDEX);
615 if (index == NO_INDEX) {
616 BLI_assert(next_allocate_index < UINT_MAX - 2);
617 vert_to_index_.add(v, next_allocate_index++);
618 }
619 }
620 }
621 int tot_v = next_allocate_index;
622 vert_ = Array<const Vert *>(tot_v);
623 for (auto item : vert_to_index_.items()) {
624 int index = item.value;
625 BLI_assert(index < tot_v);
626 vert_[index] = item.key;
627 }
628 /* Easier debugging (at least when there are no merged input verts)
629 * if output vert order is same as input, with new verts at the end.
630 * TODO: when all debugged, set fix_order = false. */
631 const bool fix_order = true;
632 if (fix_order) {
633 std::sort(vert_.begin(), vert_.end(), [](const Vert *a, const Vert *b) {
634 if (a->orig != NO_INDEX && b->orig != NO_INDEX) {
635 return a->orig < b->orig;
636 }
637 if (a->orig != NO_INDEX) {
638 return true;
639 }
640 if (b->orig != NO_INDEX) {
641 return false;
642 }
643 return a->id < b->id;
644 });
645 for (int i : vert_.index_range()) {
646 const Vert *v = vert_[i];
647 vert_to_index_.add_overwrite(v, i);
648 }
649 }
650 vert_populated_ = true;
651 }
652
erase_face_positions(int f_index,Span<bool> face_pos_erase,IMeshArena * arena)653 bool IMesh::erase_face_positions(int f_index, Span<bool> face_pos_erase, IMeshArena *arena)
654 {
655 const Face *cur_f = this->face(f_index);
656 int cur_len = cur_f->size();
657 int num_to_erase = 0;
658 for (int i : cur_f->index_range()) {
659 if (face_pos_erase[i]) {
660 ++num_to_erase;
661 }
662 }
663 if (num_to_erase == 0) {
664 return false;
665 }
666 int new_len = cur_len - num_to_erase;
667 if (new_len < 3) {
668 /* This erase causes removal of whole face.
669 * Because this may be called from a loop over the face array,
670 * we don't want to compress that array right here; instead will
671 * mark with null pointer and caller should call remove_null_faces().
672 * the loop is done.
673 */
674 this->face_[f_index] = NULL;
675 return true;
676 }
677 Array<const Vert *> new_vert(new_len);
678 Array<int> new_edge_orig(new_len);
679 Array<bool> new_is_intersect(new_len);
680 int new_index = 0;
681 for (int i : cur_f->index_range()) {
682 if (!face_pos_erase[i]) {
683 new_vert[new_index] = (*cur_f)[i];
684 new_edge_orig[new_index] = cur_f->edge_orig[i];
685 new_is_intersect[new_index] = cur_f->is_intersect[i];
686 ++new_index;
687 }
688 }
689 BLI_assert(new_index == new_len);
690 this->face_[f_index] = arena->add_face(new_vert, cur_f->orig, new_edge_orig, new_is_intersect);
691 return false;
692 }
693
remove_null_faces()694 void IMesh::remove_null_faces()
695 {
696 int64_t nullcount = 0;
697 for (Face *f : this->face_) {
698 if (f == NULL) {
699 ++nullcount;
700 }
701 }
702 if (nullcount == 0) {
703 return;
704 }
705 int64_t new_size = this->face_.size() - nullcount;
706 int64_t copy_to_index = 0;
707 int64_t copy_from_index = 0;
708 Array<Face *> new_face(new_size);
709 while (copy_from_index < face_.size()) {
710 Face *f_from = face_[copy_from_index++];
711 if (f_from) {
712 new_face[copy_to_index++] = f_from;
713 }
714 }
715 this->face_ = new_face;
716 }
717
operator <<(std::ostream & os,const IMesh & mesh)718 std::ostream &operator<<(std::ostream &os, const IMesh &mesh)
719 {
720 if (mesh.has_verts()) {
721 os << "Verts:\n";
722 int i = 0;
723 for (const Vert *v : mesh.vertices()) {
724 os << i << ": " << v << "\n";
725 ++i;
726 }
727 }
728 os << "\nFaces:\n";
729 int i = 0;
730 for (const Face *f : mesh.faces()) {
731 os << i << ": " << f << "\n";
732 if (f->plane != nullptr) {
733 os << " plane=" << f->plane << " eorig=[";
734 for (Face::FacePos p = 0; p < f->size(); ++p) {
735 os << f->edge_orig[p] << " ";
736 }
737 os << "]\n";
738 }
739 ++i;
740 }
741 return os;
742 }
743
744 struct BoundingBox {
745 float3 min{FLT_MAX, FLT_MAX, FLT_MAX};
746 float3 max{-FLT_MAX, -FLT_MAX, -FLT_MAX};
747
748 BoundingBox() = default;
BoundingBoxblender::meshintersect::BoundingBox749 BoundingBox(const float3 &min, const float3 &max) : min(min), max(max)
750 {
751 }
BoundingBoxblender::meshintersect::BoundingBox752 BoundingBox(const BoundingBox &other) : min(other.min), max(other.max)
753 {
754 }
BoundingBoxblender::meshintersect::BoundingBox755 BoundingBox(BoundingBox &&other) noexcept : min(std::move(other.min)), max(std::move(other.max))
756 {
757 }
758 ~BoundingBox() = default;
operator =blender::meshintersect::BoundingBox759 BoundingBox operator=(const BoundingBox &other)
760 {
761 if (this != &other) {
762 min = other.min;
763 max = other.max;
764 }
765 return *this;
766 }
operator =blender::meshintersect::BoundingBox767 BoundingBox operator=(BoundingBox &&other) noexcept
768 {
769 min = std::move(other.min);
770 max = std::move(other.max);
771 return *this;
772 }
773
combineblender::meshintersect::BoundingBox774 void combine(const float3 &p)
775 {
776 min.x = min_ff(min.x, p.x);
777 min.y = min_ff(min.y, p.y);
778 min.z = min_ff(min.z, p.z);
779 max.x = max_ff(max.x, p.x);
780 max.y = max_ff(max.y, p.y);
781 max.z = max_ff(max.z, p.z);
782 }
783
combineblender::meshintersect::BoundingBox784 void combine(const double3 &p)
785 {
786 min.x = min_ff(min.x, static_cast<float>(p.x));
787 min.y = min_ff(min.y, static_cast<float>(p.y));
788 min.z = min_ff(min.z, static_cast<float>(p.z));
789 max.x = max_ff(max.x, static_cast<float>(p.x));
790 max.y = max_ff(max.y, static_cast<float>(p.y));
791 max.z = max_ff(max.z, static_cast<float>(p.z));
792 }
793
combineblender::meshintersect::BoundingBox794 void combine(const BoundingBox &bb)
795 {
796 min.x = min_ff(min.x, bb.min.x);
797 min.y = min_ff(min.y, bb.min.y);
798 min.z = min_ff(min.z, bb.min.z);
799 max.x = max_ff(max.x, bb.max.x);
800 max.y = max_ff(max.y, bb.max.y);
801 max.z = max_ff(max.z, bb.max.z);
802 }
803
expandblender::meshintersect::BoundingBox804 void expand(float pad)
805 {
806 min.x -= pad;
807 min.y -= pad;
808 min.z -= pad;
809 max.x += pad;
810 max.y += pad;
811 max.z += pad;
812 }
813 };
814
815 /**
816 * Assume bounding boxes have been expanded by a sufficient epsilon on all sides
817 * so that the comparisons against the bb bounds are sufficient to guarantee that
818 * if an overlap or even touching could happen, this will return true.
819 */
bbs_might_intersect(const BoundingBox & bb_a,const BoundingBox & bb_b)820 static bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b)
821 {
822 return isect_aabb_aabb_v3(bb_a.min, bb_a.max, bb_b.min, bb_b.max);
823 }
824
825 /**
826 * Data and functions to calculate bounding boxes and pad them, in parallel.
827 * The bounding box calculation has the additional task of calculating the maximum
828 * absolute value of any coordinate in the mesh, which will be used to calculate
829 * the pad value.
830 */
831 struct BBChunkData {
832 double max_abs_val = 0.0;
833 };
834
835 struct BBCalcData {
836 const IMesh &im;
837 Array<BoundingBox> *face_bounding_box;
838
BBCalcDatablender::meshintersect::BBCalcData839 BBCalcData(const IMesh &im, Array<BoundingBox> *fbb) : im(im), face_bounding_box(fbb){};
840 };
841
calc_face_bb_range_func(void * __restrict userdata,const int iter,const TaskParallelTLS * __restrict tls)842 static void calc_face_bb_range_func(void *__restrict userdata,
843 const int iter,
844 const TaskParallelTLS *__restrict tls)
845 {
846 BBCalcData *bbdata = static_cast<BBCalcData *>(userdata);
847 double max_abs = 0.0;
848 const Face &face = *bbdata->im.face(iter);
849 BoundingBox &bb = (*bbdata->face_bounding_box)[iter];
850 for (const Vert *v : face) {
851 bb.combine(v->co);
852 for (int i = 0; i < 3; ++i) {
853 max_abs = max_dd(max_abs, fabs(v->co[i]));
854 }
855 }
856 BBChunkData *chunk = static_cast<BBChunkData *>(tls->userdata_chunk);
857 chunk->max_abs_val = max_dd(max_abs, chunk->max_abs_val);
858 }
859
860 struct BBPadData {
861 Array<BoundingBox> *face_bounding_box;
862 double pad;
863
BBPadDatablender::meshintersect::BBPadData864 BBPadData(Array<BoundingBox> *fbb, double pad) : face_bounding_box(fbb), pad(pad){};
865 };
866
pad_face_bb_range_func(void * __restrict userdata,const int iter,const TaskParallelTLS * __restrict UNUSED (tls))867 static void pad_face_bb_range_func(void *__restrict userdata,
868 const int iter,
869 const TaskParallelTLS *__restrict UNUSED(tls))
870 {
871 BBPadData *pad_data = static_cast<BBPadData *>(userdata);
872 (*pad_data->face_bounding_box)[iter].expand(pad_data->pad);
873 }
874
calc_face_bb_reduce(const void * __restrict UNUSED (userdata),void * __restrict chunk_join,void * __restrict chunk)875 static void calc_face_bb_reduce(const void *__restrict UNUSED(userdata),
876 void *__restrict chunk_join,
877 void *__restrict chunk)
878 {
879 BBChunkData *bbchunk_join = static_cast<BBChunkData *>(chunk_join);
880 BBChunkData *bbchunk = static_cast<BBChunkData *>(chunk);
881 bbchunk_join->max_abs_val = max_dd(bbchunk_join->max_abs_val, bbchunk->max_abs_val);
882 }
883
884 /**
885 * We will expand the bounding boxes by an epsilon on all sides so that
886 * the "less than" tests in isect_aabb_aabb_v3 are sufficient to detect
887 * touching or overlap.
888 */
calc_face_bounding_boxes(const IMesh & m)889 static Array<BoundingBox> calc_face_bounding_boxes(const IMesh &m)
890 {
891 int n = m.face_size();
892 Array<BoundingBox> ans(n);
893 TaskParallelSettings settings;
894 BBCalcData data(m, &ans);
895 BBChunkData chunk_data;
896 BLI_parallel_range_settings_defaults(&settings);
897 settings.userdata_chunk = &chunk_data;
898 settings.userdata_chunk_size = sizeof(chunk_data);
899 settings.func_reduce = calc_face_bb_reduce;
900 settings.min_iter_per_thread = 1000;
901 settings.use_threading = intersect_use_threading;
902 BLI_task_parallel_range(0, n, &data, calc_face_bb_range_func, &settings);
903 double max_abs_val = chunk_data.max_abs_val;
904 constexpr float pad_factor = 10.0f;
905 float pad = max_abs_val == 0.0f ? FLT_EPSILON : 2 * FLT_EPSILON * max_abs_val;
906 pad *= pad_factor; /* For extra safety. */
907 TaskParallelSettings pad_settings;
908 BLI_parallel_range_settings_defaults(&pad_settings);
909 settings.min_iter_per_thread = 1000;
910 settings.use_threading = intersect_use_threading;
911 BBPadData pad_data(&ans, pad);
912 BLI_task_parallel_range(0, n, &pad_data, pad_face_bb_range_func, &pad_settings);
913 return ans;
914 }
915
916 /**
917 * A cluster of co-planar triangles, by index.
918 * A pair of triangles T0 and T1 is said to "non-trivially co-planar-intersect"
919 * if they are co-planar, intersect, and their intersection is not just existing
920 * elements (verts, edges) of both triangles.
921 * A co-planar cluster is said to be "nontrivial" if it has more than one triangle
922 * and every triangle in it non-trivially co-planar-intersects with at least one other
923 * triangle in the cluster.
924 */
925 class CoplanarCluster {
926 Vector<int> tris_;
927 BoundingBox bb_;
928
929 public:
930 CoplanarCluster() = default;
CoplanarCluster(int t,const BoundingBox & bb)931 CoplanarCluster(int t, const BoundingBox &bb)
932 {
933 this->add_tri(t, bb);
934 }
CoplanarCluster(const CoplanarCluster & other)935 CoplanarCluster(const CoplanarCluster &other) : tris_(other.tris_), bb_(other.bb_)
936 {
937 }
CoplanarCluster(CoplanarCluster && other)938 CoplanarCluster(CoplanarCluster &&other) noexcept
939 : tris_(std::move(other.tris_)), bb_(std::move(other.bb_))
940 {
941 }
942 ~CoplanarCluster() = default;
operator =(const CoplanarCluster & other)943 CoplanarCluster &operator=(const CoplanarCluster &other)
944 {
945 if (this != &other) {
946 tris_ = other.tris_;
947 bb_ = other.bb_;
948 }
949 return *this;
950 }
operator =(CoplanarCluster && other)951 CoplanarCluster &operator=(CoplanarCluster &&other) noexcept
952 {
953 tris_ = std::move(other.tris_);
954 bb_ = std::move(other.bb_);
955 return *this;
956 }
957
958 /* Assume that caller knows this will not be a duplicate. */
add_tri(int t,const BoundingBox & bb)959 void add_tri(int t, const BoundingBox &bb)
960 {
961 tris_.append(t);
962 bb_.combine(bb);
963 }
tot_tri() const964 int tot_tri() const
965 {
966 return tris_.size();
967 }
tri(int index) const968 int tri(int index) const
969 {
970 return tris_[index];
971 }
begin() const972 const int *begin() const
973 {
974 return tris_.begin();
975 }
end() const976 const int *end() const
977 {
978 return tris_.end();
979 }
980
bounding_box() const981 const BoundingBox &bounding_box() const
982 {
983 return bb_;
984 }
985 };
986
987 /**
988 * Maintains indexed set of #CoplanarCluster, with the added ability
989 * to efficiently find the cluster index of any given triangle
990 * (the max triangle index needs to be given in the initializer).
991 * The #tri_cluster(t) function returns -1 if t is not part of any cluster.
992 */
993 class CoplanarClusterInfo {
994 Vector<CoplanarCluster> clusters_;
995 Array<int> tri_cluster_;
996
997 public:
998 CoplanarClusterInfo() = default;
CoplanarClusterInfo(int numtri)999 explicit CoplanarClusterInfo(int numtri) : tri_cluster_(Array<int>(numtri))
1000 {
1001 tri_cluster_.fill(-1);
1002 }
1003
tri_cluster(int t) const1004 int tri_cluster(int t) const
1005 {
1006 BLI_assert(t < tri_cluster_.size());
1007 return tri_cluster_[t];
1008 }
1009
add_cluster(CoplanarCluster cl)1010 int add_cluster(CoplanarCluster cl)
1011 {
1012 int c_index = clusters_.append_and_get_index(cl);
1013 for (int t : cl) {
1014 BLI_assert(t < tri_cluster_.size());
1015 tri_cluster_[t] = c_index;
1016 }
1017 return c_index;
1018 }
1019
tot_cluster() const1020 int tot_cluster() const
1021 {
1022 return clusters_.size();
1023 }
1024
begin()1025 const CoplanarCluster *begin()
1026 {
1027 return clusters_.begin();
1028 }
1029
end()1030 const CoplanarCluster *end()
1031 {
1032 return clusters_.end();
1033 }
1034
index_range() const1035 IndexRange index_range() const
1036 {
1037 return clusters_.index_range();
1038 }
1039
cluster(int index) const1040 const CoplanarCluster &cluster(int index) const
1041 {
1042 BLI_assert(index < clusters_.size());
1043 return clusters_[index];
1044 }
1045 };
1046
1047 static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl);
1048
1049 static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo);
1050
1051 enum ITT_value_kind { INONE, IPOINT, ISEGMENT, ICOPLANAR };
1052
1053 struct ITT_value {
1054 enum ITT_value_kind kind;
1055 mpq3 p1; /* Only relevant for IPOINT and ISEGMENT kind. */
1056 mpq3 p2; /* Only relevant for ISEGMENT kind. */
1057 int t_source; /* Index of the source triangle that intersected the target one. */
1058
ITT_valueblender::meshintersect::ITT_value1059 ITT_value() : kind(INONE), t_source(-1)
1060 {
1061 }
ITT_valueblender::meshintersect::ITT_value1062 ITT_value(ITT_value_kind k) : kind(k), t_source(-1)
1063 {
1064 }
ITT_valueblender::meshintersect::ITT_value1065 ITT_value(ITT_value_kind k, int tsrc) : kind(k), t_source(tsrc)
1066 {
1067 }
ITT_valueblender::meshintersect::ITT_value1068 ITT_value(ITT_value_kind k, const mpq3 &p1) : kind(k), p1(p1), t_source(-1)
1069 {
1070 }
ITT_valueblender::meshintersect::ITT_value1071 ITT_value(ITT_value_kind k, const mpq3 &p1, const mpq3 &p2)
1072 : kind(k), p1(p1), p2(p2), t_source(-1)
1073 {
1074 }
ITT_valueblender::meshintersect::ITT_value1075 ITT_value(const ITT_value &other)
1076 : kind(other.kind), p1(other.p1), p2(other.p2), t_source(other.t_source)
1077 {
1078 }
ITT_valueblender::meshintersect::ITT_value1079 ITT_value(ITT_value &&other) noexcept
1080 : kind(other.kind),
1081 p1(std::move(other.p1)),
1082 p2(std::move(other.p2)),
1083 t_source(other.t_source)
1084 {
1085 }
~ITT_valueblender::meshintersect::ITT_value1086 ~ITT_value()
1087 {
1088 }
operator =blender::meshintersect::ITT_value1089 ITT_value &operator=(const ITT_value &other)
1090 {
1091 if (this != &other) {
1092 kind = other.kind;
1093 p1 = other.p1;
1094 p2 = other.p2;
1095 t_source = other.t_source;
1096 }
1097 return *this;
1098 }
operator =blender::meshintersect::ITT_value1099 ITT_value &operator=(ITT_value &&other) noexcept
1100 {
1101 kind = other.kind;
1102 p1 = std::move(other.p1);
1103 p2 = std::move(other.p2);
1104 t_source = other.t_source;
1105 return *this;
1106 }
1107 };
1108
1109 static std::ostream &operator<<(std::ostream &os, const ITT_value &itt);
1110
1111 /**
1112 * Project a 3d vert to a 2d one by eliding proj_axis. This does not create
1113 * degeneracies as long as the projection axis is one where the corresponding
1114 * component of the originating plane normal is non-zero.
1115 */
project_3d_to_2d(const mpq3 & p3d,int proj_axis)1116 static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
1117 {
1118 mpq2 p2d;
1119 switch (proj_axis) {
1120 case (0): {
1121 p2d[0] = p3d[1];
1122 p2d[1] = p3d[2];
1123 break;
1124 }
1125 case (1): {
1126 p2d[0] = p3d[0];
1127 p2d[1] = p3d[2];
1128 break;
1129 }
1130 case (2): {
1131 p2d[0] = p3d[0];
1132 p2d[1] = p3d[1];
1133 break;
1134 }
1135 default:
1136 BLI_assert(false);
1137 }
1138 return p2d;
1139 }
1140
1141 /**
1142 * The sup and index functions are defined in the paper:
1143 * EXACT GEOMETRIC COMPUTATION USING CASCADING, by
1144 * Burnikel, Funke, and Seel. They are used to find absolute
1145 * bounds on the error due to doing a calculation in double
1146 * instead of exactly. For calculations involving only +, -, and *,
1147 * the supremum is the same function except using absolute values
1148 * on inputs and using + instead of -.
1149 * The index function follows these rules:
1150 * index(x op y) = 1 + max(index(x), index(y)) for op + or -
1151 * index(x * y) = 1 + index(x) + index(y)
1152 * index(x) = 0 if input x can be represented exactly as a double
1153 * index(x) = 1 otherwise.
1154 *
1155 * With these rules in place, we know an absolute error bound:
1156 *
1157 * |E_exact - E| <= supremum(E) * index(E) * DBL_EPSILON
1158 *
1159 * where E_exact is what would have been the exact value of the
1160 * expression and E is the one calculated with doubles.
1161 *
1162 * So the sign of E is the same as the sign of E_exact if
1163 * |E| > supremum(E) * index(E) * DBL_EPSILON
1164 *
1165 * Note: a possible speedup would be to have a simple function
1166 * that calculates the error bound if one knows that all values
1167 * are less than some global maximum - most of the function would
1168 * be calculated ahead of time. The global max could be passed
1169 * from above.
1170 */
supremum_dot_cross(const double3 & a,const double3 & b)1171 static double supremum_dot_cross(const double3 &a, const double3 &b)
1172 {
1173 double3 abs_a = double3::abs(a);
1174 double3 abs_b = double3::abs(b);
1175 double3 c;
1176 /* This is dot(cross(a, b), cross(a,b)) but using absolute values for a and b
1177 * and always using + when operation is + or -. */
1178 c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1];
1179 c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2];
1180 c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0];
1181 return double3::dot(c, c);
1182 }
1183
1184 /* The index of dot when inputs are plane_coords with index 1 is much higher.
1185 * Plane coords have index 6.
1186 */
1187 constexpr int index_dot_plane_coords = 15;
1188
1189 /**
1190 * Used with supremum to get error bound. See Burnikel et al paper.
1191 * index_plane_coord is the index of a plane coordinate calculated
1192 * for a triangle using the usual formula, assuming the input
1193 * coordinates have index 1.
1194 * index_cross is the index of each coordinate of the cross product.
1195 * It is actually 2 + 2 * (max index of input coords).
1196 * index_dot_cross is the index of the dot product of two cross products.
1197 * It is actually 7 + 4 * (max index of input coords)
1198 */
1199 constexpr int index_dot_cross = 11;
1200
1201 /**
1202 * Return the approximate side of point p on a plane with normal plane_no and point plane_p.
1203 * The answer will be 1 if p is definitely above the plane, -1 if it is definitely below.
1204 * If the answer is 0, we are unsure about which side of the plane (or if it is on the plane).
1205 * In exact arithmetic, the answer is just `sgn(dot(p - plane_p, plane_no))`.
1206 *
1207 * The plane_no input is constructed, so has a higher index.
1208 */
1209 constexpr int index_plane_side = 3 + 2 * index_dot_plane_coords;
1210
filter_plane_side(const double3 & p,const double3 & plane_p,const double3 & plane_no,const double3 & abs_p,const double3 & abs_plane_p,const double3 & abs_plane_no)1211 static int filter_plane_side(const double3 &p,
1212 const double3 &plane_p,
1213 const double3 &plane_no,
1214 const double3 &abs_p,
1215 const double3 &abs_plane_p,
1216 const double3 &abs_plane_no)
1217 {
1218 double d = double3::dot(p - plane_p, plane_no);
1219 if (d == 0.0) {
1220 return 0;
1221 }
1222 double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no);
1223 double err_bound = supremum * index_plane_side * DBL_EPSILON;
1224 if (fabs(d) > err_bound) {
1225 return d > 0 ? 1 : -1;
1226 }
1227 return 0;
1228 }
1229
1230 /*
1231 * interesect_tri_tri and helper functions.
1232 * This code uses the algorithm of Guigue and Devillers, as described
1233 * in "Faster Triangle-Triangle Intersection Tests".
1234 * Adapted from github code by Eric Haines:
1235 * github.com/erich666/jgt-code/tree/master/Volume_08/Number_1/Guigue2003
1236 */
1237
1238 /**
1239 * Return the point on ab where the plane with normal n containing point c intersects it.
1240 * Assumes ab is not perpendicular to n.
1241 * This works because the ratio of the projections of ab and ac onto n is the same as
1242 * the ratio along the line ab of the intersection point to the whole of ab.
1243 */
tti_interp(const mpq3 & a,const mpq3 & b,const mpq3 & c,const mpq3 & n)1244 static inline mpq3 tti_interp(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &n)
1245 {
1246 mpq3 ab = a - b;
1247 mpq_class den = mpq3::dot(ab, n);
1248 BLI_assert(den != 0);
1249 mpq_class alpha = mpq3::dot(a - c, n) / den;
1250 return a - alpha * ab;
1251 }
1252
1253 /**
1254 * Return +1, 0, -1 as a + ad is above, on, or below the oriented plane containing a, b, c in CCW
1255 * order. This is the same as -oriented(a, b, c, a + ad), but uses fewer arithmetic operations.
1256 * TODO: change arguments to `const Vert *` and use floating filters.
1257 */
tti_above(const mpq3 & a,const mpq3 & b,const mpq3 & c,const mpq3 & ad)1258 static inline int tti_above(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &ad)
1259 {
1260 mpq3 n = mpq3::cross(b - a, c - a);
1261 return sgn(mpq3::dot(ad, n));
1262 }
1263
1264 /**
1265 * Given that triangles (p1, q1, r1) and (p2, q2, r2) are in canonical order,
1266 * use the classification chart in the Guigue and Devillers paper to find out
1267 * how the intervals [i,j] and [k,l] overlap, where [i,j] is where p1r1 and p1q1
1268 * intersect the plane-plane intersection line, L, and [k,l] is where p2q2 and p2r2
1269 * intersect L. By the canonicalization, those segments intersect L exactly once.
1270 * Canonicalization has made it so that for p1, q1, r1, either:
1271 * (a)) p1 is off the second triangle's plane and both q1 and r1 are either
1272 * on the plane or on the other side of it from p1; or
1273 * (b) p1 is on the plane both q1 and r1 are on the same side
1274 * of the plane and at least one of q1 and r1 are off the plane.
1275 * Similarly for p2, q2, r2 with respect to the first triangle's plane.
1276 */
itt_canon2(const mpq3 & p1,const mpq3 & q1,const mpq3 & r1,const mpq3 & p2,const mpq3 & q2,const mpq3 & r2,const mpq3 & n1,const mpq3 & n2)1277 static ITT_value itt_canon2(const mpq3 &p1,
1278 const mpq3 &q1,
1279 const mpq3 &r1,
1280 const mpq3 &p2,
1281 const mpq3 &q2,
1282 const mpq3 &r2,
1283 const mpq3 &n1,
1284 const mpq3 &n2)
1285 {
1286 constexpr int dbg_level = 0;
1287 if (dbg_level > 0) {
1288 std::cout << "\ntri_tri_intersect_canon:\n";
1289 std::cout << "p1=" << p1 << " q1=" << q1 << " r1=" << r1 << "\n";
1290 std::cout << "p2=" << p2 << " q2=" << q2 << " r2=" << r2 << "\n";
1291 std::cout << "n1=" << n1 << " n2=" << n2 << "\n";
1292 std::cout << "approximate values:\n";
1293 std::cout << "p1=(" << p1[0].get_d() << "," << p1[1].get_d() << "," << p1[2].get_d() << ")\n";
1294 std::cout << "q1=(" << q1[0].get_d() << "," << q1[1].get_d() << "," << q1[2].get_d() << ")\n";
1295 std::cout << "r1=(" << r1[0].get_d() << "," << r1[1].get_d() << "," << r1[2].get_d() << ")\n";
1296 std::cout << "p2=(" << p2[0].get_d() << "," << p2[1].get_d() << "," << p2[2].get_d() << ")\n";
1297 std::cout << "q2=(" << q2[0].get_d() << "," << q2[1].get_d() << "," << q2[2].get_d() << ")\n";
1298 std::cout << "r2=(" << r2[0].get_d() << "," << r2[1].get_d() << "," << r2[2].get_d() << ")\n";
1299 std::cout << "n1=(" << n1[0].get_d() << "," << n1[1].get_d() << "," << n1[2].get_d() << ")\n";
1300 std::cout << "n2=(" << n2[0].get_d() << "," << n2[1].get_d() << "," << n2[2].get_d() << ")\n";
1301 }
1302 mpq3 p1p2 = p2 - p1;
1303 mpq3 intersect_1;
1304 mpq3 intersect_2;
1305 bool no_overlap = false;
1306 /* Top test in classification tree. */
1307 if (tti_above(p1, q1, r2, p1p2) > 0) {
1308 /* Middle right test in classification tree. */
1309 if (tti_above(p1, r1, r2, p1p2) <= 0) {
1310 /* Bottom right test in classification tree. */
1311 if (tti_above(p1, r1, q2, p1p2) > 0) {
1312 /* Overlap is [k [i l] j]. */
1313 if (dbg_level > 0) {
1314 std::cout << "overlap [k [i l] j]\n";
1315 }
1316 /* i is intersect with p1r1. l is intersect with p2r2. */
1317 intersect_1 = tti_interp(p1, r1, p2, n2);
1318 intersect_2 = tti_interp(p2, r2, p1, n1);
1319 }
1320 else {
1321 /* Overlap is [i [k l] j]. */
1322 if (dbg_level > 0) {
1323 std::cout << "overlap [i [k l] j]\n";
1324 }
1325 /* k is intersect with p2q2. l is intersect is p2r2. */
1326 intersect_1 = tti_interp(p2, q2, p1, n1);
1327 intersect_2 = tti_interp(p2, r2, p1, n1);
1328 }
1329 }
1330 else {
1331 /* No overlap: [k l] [i j]. */
1332 if (dbg_level > 0) {
1333 std::cout << "no overlap: [k l] [i j]\n";
1334 }
1335 no_overlap = true;
1336 }
1337 }
1338 else {
1339 /* Middle left test in classification tree. */
1340 if (tti_above(p1, q1, q2, p1p2) < 0) {
1341 /* No overlap: [i j] [k l]. */
1342 if (dbg_level > 0) {
1343 std::cout << "no overlap: [i j] [k l]\n";
1344 }
1345 no_overlap = true;
1346 }
1347 else {
1348 /* Bottom left test in classification tree. */
1349 if (tti_above(p1, r1, q2, p1p2) >= 0) {
1350 /* Overlap is [k [i j] l]. */
1351 if (dbg_level > 0) {
1352 std::cout << "overlap [k [i j] l]\n";
1353 }
1354 /* i is intersect with p1r1. j is intersect with p1q1. */
1355 intersect_1 = tti_interp(p1, r1, p2, n2);
1356 intersect_2 = tti_interp(p1, q1, p2, n2);
1357 }
1358 else {
1359 /* Overlap is [i [k j] l]. */
1360 if (dbg_level > 0) {
1361 std::cout << "overlap [i [k j] l]\n";
1362 }
1363 /* k is intersect with p2q2. j is intersect with p1q1. */
1364 intersect_1 = tti_interp(p2, q2, p1, n1);
1365 intersect_2 = tti_interp(p1, q1, p2, n2);
1366 }
1367 }
1368 }
1369 if (no_overlap) {
1370 return ITT_value(INONE);
1371 }
1372 if (intersect_1 == intersect_2) {
1373 if (dbg_level > 0) {
1374 std::cout << "single intersect: " << intersect_1 << "\n";
1375 }
1376 return ITT_value(IPOINT, intersect_1);
1377 }
1378 if (dbg_level > 0) {
1379 std::cout << "intersect segment: " << intersect_1 << ", " << intersect_2 << "\n";
1380 }
1381 return ITT_value(ISEGMENT, intersect_1, intersect_2);
1382 }
1383
1384 /* Helper function for intersect_tri_tri. Args have been canonicalized for triangle 1. */
1385
itt_canon1(const mpq3 & p1,const mpq3 & q1,const mpq3 & r1,const mpq3 & p2,const mpq3 & q2,const mpq3 & r2,const mpq3 & n1,const mpq3 & n2,int sp2,int sq2,int sr2)1386 static ITT_value itt_canon1(const mpq3 &p1,
1387 const mpq3 &q1,
1388 const mpq3 &r1,
1389 const mpq3 &p2,
1390 const mpq3 &q2,
1391 const mpq3 &r2,
1392 const mpq3 &n1,
1393 const mpq3 &n2,
1394 int sp2,
1395 int sq2,
1396 int sr2)
1397 {
1398 constexpr int dbg_level = 0;
1399 if (sp2 > 0) {
1400 if (sq2 > 0) {
1401 return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2);
1402 }
1403 if (sr2 > 0) {
1404 return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2);
1405 }
1406 return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2);
1407 }
1408 if (sp2 < 0) {
1409 if (sq2 < 0) {
1410 return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2);
1411 }
1412 if (sr2 < 0) {
1413 return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2);
1414 }
1415 return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2);
1416 }
1417 if (sq2 < 0) {
1418 if (sr2 >= 0) {
1419 return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2);
1420 }
1421 return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2);
1422 }
1423 if (sq2 > 0) {
1424 if (sr2 > 0) {
1425 return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2);
1426 }
1427 return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2);
1428 }
1429 if (sr2 > 0) {
1430 return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2);
1431 }
1432 if (sr2 < 0) {
1433 return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2);
1434 }
1435 if (dbg_level > 0) {
1436 std::cout << "triangles are co-planar\n";
1437 }
1438 return ITT_value(ICOPLANAR);
1439 }
1440
intersect_tri_tri(const IMesh & tm,int t1,int t2)1441 static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
1442 {
1443 constexpr int dbg_level = 0;
1444 # ifdef PERFDEBUG
1445 incperfcount(1); /* Intersect_tri_tri calls. */
1446 # endif
1447 const Face &tri1 = *tm.face(t1);
1448 const Face &tri2 = *tm.face(t2);
1449 BLI_assert(tri1.plane_populated() && tri2.plane_populated());
1450 const Vert *vp1 = tri1[0];
1451 const Vert *vq1 = tri1[1];
1452 const Vert *vr1 = tri1[2];
1453 const Vert *vp2 = tri2[0];
1454 const Vert *vq2 = tri2[1];
1455 const Vert *vr2 = tri2[2];
1456 if (dbg_level > 0) {
1457 std::cout << "\nINTERSECT_TRI_TRI t1=" << t1 << ", t2=" << t2 << "\n";
1458 std::cout << " p1 = " << vp1 << "\n";
1459 std::cout << " q1 = " << vq1 << "\n";
1460 std::cout << " r1 = " << vr1 << "\n";
1461 std::cout << " p2 = " << vp2 << "\n";
1462 std::cout << " q2 = " << vq2 << "\n";
1463 std::cout << " r2 = " << vr2 << "\n";
1464 }
1465
1466 /* Get signs of t1's vertices' distances to plane of t2 and vice versa. */
1467
1468 /* Try first getting signs with double arithmetic, with error bounds.
1469 * If the signs calculated in this section are not 0, they are the same
1470 * as what they would be using exact arithmetic. */
1471 const double3 &d_p1 = vp1->co;
1472 const double3 &d_q1 = vq1->co;
1473 const double3 &d_r1 = vr1->co;
1474 const double3 &d_p2 = vp2->co;
1475 const double3 &d_q2 = vq2->co;
1476 const double3 &d_r2 = vr2->co;
1477 const double3 &d_n2 = tri2.plane->norm;
1478
1479 const double3 &abs_d_p1 = double3::abs(d_p1);
1480 const double3 &abs_d_q1 = double3::abs(d_q1);
1481 const double3 &abs_d_r1 = double3::abs(d_r1);
1482 const double3 &abs_d_r2 = double3::abs(d_r2);
1483 const double3 &abs_d_n2 = double3::abs(d_n2);
1484
1485 int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2);
1486 int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2);
1487 int sr1 = filter_plane_side(d_r1, d_r2, d_n2, abs_d_r1, abs_d_r2, abs_d_n2);
1488 if ((sp1 > 0 && sq1 > 0 && sr1 > 0) || (sp1 < 0 && sq1 < 0 && sr1 < 0)) {
1489 # ifdef PERFDEBUG
1490 incperfcount(2); /* Tri tri intersects decided by filter plane tests. */
1491 # endif
1492 if (dbg_level > 0) {
1493 std::cout << "no intersection, all t1's verts above or below t2\n";
1494 }
1495 return ITT_value(INONE);
1496 }
1497
1498 const double3 &d_n1 = tri1.plane->norm;
1499 const double3 &abs_d_p2 = double3::abs(d_p2);
1500 const double3 &abs_d_q2 = double3::abs(d_q2);
1501 const double3 &abs_d_n1 = double3::abs(d_n1);
1502
1503 int sp2 = filter_plane_side(d_p2, d_r1, d_n1, abs_d_p2, abs_d_r1, abs_d_n1);
1504 int sq2 = filter_plane_side(d_q2, d_r1, d_n1, abs_d_q2, abs_d_r1, abs_d_n1);
1505 int sr2 = filter_plane_side(d_r2, d_r1, d_n1, abs_d_r2, abs_d_r1, abs_d_n1);
1506 if ((sp2 > 0 && sq2 > 0 && sr2 > 0) || (sp2 < 0 && sq2 < 0 && sr2 < 0)) {
1507 # ifdef PERFDEBUG
1508 incperfcount(2); /* Tri tri intersects decided by filter plane tests. */
1509 # endif
1510 if (dbg_level > 0) {
1511 std::cout << "no intersection, all t2's verts above or below t1\n";
1512 }
1513 return ITT_value(INONE);
1514 }
1515
1516 const mpq3 &p1 = vp1->co_exact;
1517 const mpq3 &q1 = vq1->co_exact;
1518 const mpq3 &r1 = vr1->co_exact;
1519 const mpq3 &p2 = vp2->co_exact;
1520 const mpq3 &q2 = vq2->co_exact;
1521 const mpq3 &r2 = vr2->co_exact;
1522
1523 const mpq3 &n2 = tri2.plane->norm_exact;
1524 if (sp1 == 0) {
1525 sp1 = sgn(mpq3::dot(p1 - r2, n2));
1526 }
1527 if (sq1 == 0) {
1528 sq1 = sgn(mpq3::dot(q1 - r2, n2));
1529 }
1530 if (sr1 == 0) {
1531 sr1 = sgn(mpq3::dot(r1 - r2, n2));
1532 }
1533
1534 if (dbg_level > 1) {
1535 std::cout << " sp1=" << sp1 << " sq1=" << sq1 << " sr1=" << sr1 << "\n";
1536 }
1537
1538 if ((sp1 * sq1 > 0) && (sp1 * sr1 > 0)) {
1539 if (dbg_level > 0) {
1540 std::cout << "no intersection, all t1's verts above or below t2 (exact)\n";
1541 }
1542 # ifdef PERFDEBUG
1543 incperfcount(3); /* Tri tri intersects decided by exact plane tests. */
1544 # endif
1545 return ITT_value(INONE);
1546 }
1547
1548 /* Repeat for signs of t2's vertices with respect to plane of t1. */
1549 const mpq3 &n1 = tri1.plane->norm_exact;
1550 if (sp2 == 0) {
1551 sp2 = sgn(mpq3::dot(p2 - r1, n1));
1552 }
1553 if (sq2 == 0) {
1554 sq2 = sgn(mpq3::dot(q2 - r1, n1));
1555 }
1556 if (sr2 == 0) {
1557 sr2 = sgn(mpq3::dot(r2 - r1, n1));
1558 }
1559
1560 if (dbg_level > 1) {
1561 std::cout << " sp2=" << sp2 << " sq2=" << sq2 << " sr2=" << sr2 << "\n";
1562 }
1563
1564 if ((sp2 * sq2 > 0) && (sp2 * sr2 > 0)) {
1565 if (dbg_level > 0) {
1566 std::cout << "no intersection, all t2's verts above or below t1 (exact)\n";
1567 }
1568 # ifdef PERFDEBUG
1569 incperfcount(3); /* Tri tri intersects decided by exact plane tests. */
1570 # endif
1571 return ITT_value(INONE);
1572 }
1573
1574 /* Do rest of the work with vertices in a canonical order, where p1 is on
1575 * positive side of plane and q1, r1 are not, or p1 is on the plane and
1576 * q1 and r1 are off the plane on the same side. */
1577 ITT_value ans;
1578 if (sp1 > 0) {
1579 if (sq1 > 0) {
1580 ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1581 }
1582 else if (sr1 > 0) {
1583 ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1584 }
1585 else {
1586 ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1587 }
1588 }
1589 else if (sp1 < 0) {
1590 if (sq1 < 0) {
1591 ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1592 }
1593 else if (sr1 < 0) {
1594 ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1595 }
1596 else {
1597 ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1598 }
1599 }
1600 else {
1601 if (sq1 < 0) {
1602 if (sr1 >= 0) {
1603 ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1604 }
1605 else {
1606 ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1607 }
1608 }
1609 else if (sq1 > 0) {
1610 if (sr1 > 0) {
1611 ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1612 }
1613 else {
1614 ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1615 }
1616 }
1617 else {
1618 if (sr1 > 0) {
1619 ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1620 }
1621 else if (sr1 < 0) {
1622 ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1623 }
1624 else {
1625 if (dbg_level > 0) {
1626 std::cout << "triangles are co-planar\n";
1627 }
1628 ans = ITT_value(ICOPLANAR);
1629 }
1630 }
1631 }
1632 if (ans.kind == ICOPLANAR) {
1633 ans.t_source = t2;
1634 }
1635
1636 # ifdef PERFDEBUG
1637 if (ans.kind != INONE) {
1638 incperfcount(4);
1639 }
1640 # endif
1641 return ans;
1642 }
1643
1644 struct CDT_data {
1645 const Plane *t_plane;
1646 Vector<mpq2> vert;
1647 Vector<std::pair<int, int>> edge;
1648 Vector<Vector<int>> face;
1649 /** Parallels face, gives id from input #IMesh of input face. */
1650 Vector<int> input_face;
1651 /** Parallels face, says if input face orientation is opposite. */
1652 Vector<bool> is_reversed;
1653 /** Result of running CDT on input with (vert, edge, face). */
1654 CDT_result<mpq_class> cdt_out;
1655 int proj_axis;
1656 };
1657
1658 /**
1659 * We could de-duplicate verts here, but CDT routine will do that anyway.
1660 */
prepare_need_vert(CDT_data & cd,const mpq3 & p3d)1661 static int prepare_need_vert(CDT_data &cd, const mpq3 &p3d)
1662 {
1663 mpq2 p2d = project_3d_to_2d(p3d, cd.proj_axis);
1664 int v = cd.vert.append_and_get_index(p2d);
1665 return v;
1666 }
1667
1668 /**
1669 * To un-project a 2d vert that was projected along cd.proj_axis, we copy the coordinates
1670 * from the two axes not involved in the projection, and use the plane equation of the
1671 * originating 3d plane, cd.t_plane, to derive the coordinate of the projected axis.
1672 * The plane equation says a point p is on the plane if dot(p, plane.n()) + plane.d() == 0.
1673 * Assume that the projection axis is such that plane.n()[proj_axis] != 0.
1674 */
unproject_cdt_vert(const CDT_data & cd,const mpq2 & p2d)1675 static mpq3 unproject_cdt_vert(const CDT_data &cd, const mpq2 &p2d)
1676 {
1677 mpq3 p3d;
1678 BLI_assert(cd.t_plane->exact_populated());
1679 BLI_assert(cd.t_plane->norm_exact[cd.proj_axis] != 0);
1680 const mpq3 &n = cd.t_plane->norm_exact;
1681 const mpq_class &d = cd.t_plane->d_exact;
1682 switch (cd.proj_axis) {
1683 case (0): {
1684 mpq_class num = n[1] * p2d[0] + n[2] * p2d[1] + d;
1685 num = -num;
1686 p3d[0] = num / n[0];
1687 p3d[1] = p2d[0];
1688 p3d[2] = p2d[1];
1689 break;
1690 }
1691 case (1): {
1692 p3d[0] = p2d[0];
1693 mpq_class num = n[0] * p2d[0] + n[2] * p2d[1] + d;
1694 num = -num;
1695 p3d[1] = num / n[1];
1696 p3d[2] = p2d[1];
1697 break;
1698 }
1699 case (2): {
1700 p3d[0] = p2d[0];
1701 p3d[1] = p2d[1];
1702 mpq_class num = n[0] * p2d[0] + n[1] * p2d[1] + d;
1703 num = -num;
1704 p3d[2] = num / n[2];
1705 break;
1706 }
1707 default:
1708 BLI_assert(false);
1709 }
1710 return p3d;
1711 }
1712
prepare_need_edge(CDT_data & cd,const mpq3 & p1,const mpq3 & p2)1713 static void prepare_need_edge(CDT_data &cd, const mpq3 &p1, const mpq3 &p2)
1714 {
1715 int v1 = prepare_need_vert(cd, p1);
1716 int v2 = prepare_need_vert(cd, p2);
1717 cd.edge.append(std::pair<int, int>(v1, v2));
1718 }
1719
prepare_need_tri(CDT_data & cd,const IMesh & tm,int t)1720 static void prepare_need_tri(CDT_data &cd, const IMesh &tm, int t)
1721 {
1722 const Face &tri = *tm.face(t);
1723 int v0 = prepare_need_vert(cd, tri[0]->co_exact);
1724 int v1 = prepare_need_vert(cd, tri[1]->co_exact);
1725 int v2 = prepare_need_vert(cd, tri[2]->co_exact);
1726 bool rev;
1727 /* How to get CCW orientation of projected triangle? Note that when look down y axis
1728 * as opposed to x or z, the orientation of the other two axes is not right-and-up. */
1729 BLI_assert(cd.t_plane->exact_populated());
1730 if (tri.plane->norm_exact[cd.proj_axis] >= 0) {
1731 rev = cd.proj_axis == 1;
1732 }
1733 else {
1734 rev = cd.proj_axis != 1;
1735 }
1736 int cd_t = cd.face.append_and_get_index(Vector<int>());
1737 cd.face[cd_t].append(v0);
1738 if (rev) {
1739 cd.face[cd_t].append(v2);
1740 cd.face[cd_t].append(v1);
1741 }
1742 else {
1743 cd.face[cd_t].append(v1);
1744 cd.face[cd_t].append(v2);
1745 }
1746 cd.input_face.append(t);
1747 cd.is_reversed.append(rev);
1748 }
1749
prepare_cdt_input(const IMesh & tm,int t,const Vector<ITT_value> itts)1750 static CDT_data prepare_cdt_input(const IMesh &tm, int t, const Vector<ITT_value> itts)
1751 {
1752 CDT_data ans;
1753 BLI_assert(tm.face(t)->plane_populated());
1754 ans.t_plane = tm.face(t)->plane;
1755 BLI_assert(ans.t_plane->exact_populated());
1756 ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
1757 prepare_need_tri(ans, tm, t);
1758 for (const ITT_value &itt : itts) {
1759 switch (itt.kind) {
1760 case INONE:
1761 break;
1762 case IPOINT: {
1763 prepare_need_vert(ans, itt.p1);
1764 break;
1765 }
1766 case ISEGMENT: {
1767 prepare_need_edge(ans, itt.p1, itt.p2);
1768 break;
1769 }
1770 case ICOPLANAR: {
1771 prepare_need_tri(ans, tm, itt.t_source);
1772 break;
1773 }
1774 }
1775 }
1776 return ans;
1777 }
1778
prepare_cdt_input_for_cluster(const IMesh & tm,const CoplanarClusterInfo & clinfo,int c,const Vector<ITT_value> itts)1779 static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm,
1780 const CoplanarClusterInfo &clinfo,
1781 int c,
1782 const Vector<ITT_value> itts)
1783 {
1784 CDT_data ans;
1785 BLI_assert(c < clinfo.tot_cluster());
1786 const CoplanarCluster &cl = clinfo.cluster(c);
1787 BLI_assert(cl.tot_tri() > 0);
1788 int t0 = cl.tri(0);
1789 BLI_assert(tm.face(t0)->plane_populated());
1790 ans.t_plane = tm.face(t0)->plane;
1791 BLI_assert(ans.t_plane->exact_populated());
1792 ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
1793 for (const int t : cl) {
1794 prepare_need_tri(ans, tm, t);
1795 }
1796 for (const ITT_value &itt : itts) {
1797 switch (itt.kind) {
1798 case IPOINT: {
1799 prepare_need_vert(ans, itt.p1);
1800 break;
1801 }
1802 case ISEGMENT: {
1803 prepare_need_edge(ans, itt.p1, itt.p2);
1804 break;
1805 }
1806 default:
1807 break;
1808 }
1809 }
1810 return ans;
1811 }
1812
1813 /**
1814 * Fills in cd.cdt_out with result of doing the cdt calculation on (vert, edge, face).
1815 */
do_cdt(CDT_data & cd)1816 static void do_cdt(CDT_data &cd)
1817 {
1818 constexpr int dbg_level = 0;
1819 CDT_input<mpq_class> cdt_in;
1820 cdt_in.vert = Span<mpq2>(cd.vert);
1821 cdt_in.edge = Span<std::pair<int, int>>(cd.edge);
1822 cdt_in.face = Span<Vector<int>>(cd.face);
1823 if (dbg_level > 0) {
1824 std::cout << "CDT input\nVerts:\n";
1825 for (int i : cdt_in.vert.index_range()) {
1826 std::cout << "v" << i << ": " << cdt_in.vert[i] << "=(" << cdt_in.vert[i][0].get_d() << ","
1827 << cdt_in.vert[i][1].get_d() << ")\n";
1828 }
1829 std::cout << "Edges:\n";
1830 for (int i : cdt_in.edge.index_range()) {
1831 std::cout << "e" << i << ": (" << cdt_in.edge[i].first << ", " << cdt_in.edge[i].second
1832 << ")\n";
1833 }
1834 std::cout << "Tris\n";
1835 for (int f : cdt_in.face.index_range()) {
1836 std::cout << "f" << f << ": ";
1837 for (int j : cdt_in.face[f].index_range()) {
1838 std::cout << cdt_in.face[f][j] << " ";
1839 }
1840 std::cout << "\n";
1841 }
1842 }
1843 cdt_in.epsilon = 0; /* TODO: needs attention for non-exact T. */
1844 cd.cdt_out = blender::meshintersect::delaunay_2d_calc(cdt_in, CDT_INSIDE);
1845 if (dbg_level > 0) {
1846 std::cout << "\nCDT result\nVerts:\n";
1847 for (int i : cd.cdt_out.vert.index_range()) {
1848 std::cout << "v" << i << ": " << cd.cdt_out.vert[i] << "=(" << cd.cdt_out.vert[i][0].get_d()
1849 << "," << cd.cdt_out.vert[i][1].get_d() << "\n";
1850 }
1851 std::cout << "Tris\n";
1852 for (int f : cd.cdt_out.face.index_range()) {
1853 std::cout << "f" << f << ": ";
1854 for (int j : cd.cdt_out.face[f].index_range()) {
1855 std::cout << cd.cdt_out.face[f][j] << " ";
1856 }
1857 std::cout << "orig: ";
1858 for (int j : cd.cdt_out.face_orig[f].index_range()) {
1859 std::cout << cd.cdt_out.face_orig[f][j] << " ";
1860 }
1861 std::cout << "\n";
1862 }
1863 std::cout << "Edges\n";
1864 for (int e : cd.cdt_out.edge.index_range()) {
1865 std::cout << "e" << e << ": (" << cd.cdt_out.edge[e].first << ", "
1866 << cd.cdt_out.edge[e].second << ") ";
1867 std::cout << "orig: ";
1868 for (int j : cd.cdt_out.edge_orig[e].index_range()) {
1869 std::cout << cd.cdt_out.edge_orig[e][j] << " ";
1870 }
1871 std::cout << "\n";
1872 }
1873 }
1874 }
1875
get_cdt_edge_orig(int i0,int i1,const CDT_data & cd,const IMesh & in_tm,bool * r_is_intersect)1876 static int get_cdt_edge_orig(
1877 int i0, int i1, const CDT_data &cd, const IMesh &in_tm, bool *r_is_intersect)
1878 {
1879 int foff = cd.cdt_out.face_edge_offset;
1880 *r_is_intersect = false;
1881 for (int e : cd.cdt_out.edge.index_range()) {
1882 std::pair<int, int> edge = cd.cdt_out.edge[e];
1883 if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) {
1884 /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */
1885 /* TODO: if edge has origs from more than on part of the nary input,
1886 * then want to set *r_is_intersect to true. */
1887 for (int orig_index : cd.cdt_out.edge_orig[e]) {
1888 /* orig_index encodes the triangle and pos within the triangle of the input edge. */
1889 if (orig_index >= foff) {
1890 int in_face_index = (orig_index / foff) - 1;
1891 int pos = orig_index % foff;
1892 /* We need to retrieve the edge orig field from the Face used to populate the
1893 * in_face_index'th face of the CDT, at the pos'th position of the face. */
1894 int in_tm_face_index = cd.input_face[in_face_index];
1895 BLI_assert(in_tm_face_index < in_tm.face_size());
1896 const Face *facep = in_tm.face(in_tm_face_index);
1897 BLI_assert(pos < facep->size());
1898 bool is_rev = cd.is_reversed[in_face_index];
1899 int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
1900 if (eorig != NO_INDEX) {
1901 return eorig;
1902 }
1903 }
1904 else {
1905 /* This edge came from an edge input to the CDT problem,
1906 * so it is an intersect edge. */
1907 *r_is_intersect = true;
1908 /* TODO: maybe there is an orig index:
1909 * This happens if an input edge was formed by an input face having
1910 * an edge that is co-planar with the cluster, while the face as a whole is not. */
1911 return NO_INDEX;
1912 }
1913 }
1914 return NO_INDEX;
1915 }
1916 }
1917 return NO_INDEX;
1918 }
1919
1920 /**
1921 * Using the result of CDT in cd.cdt_out, extract an #IMesh representing the subdivision
1922 * of input triangle t, which should be an element of cd.input_face.
1923 */
extract_subdivided_tri(const CDT_data & cd,const IMesh & in_tm,int t,IMeshArena * arena)1924 static IMesh extract_subdivided_tri(const CDT_data &cd,
1925 const IMesh &in_tm,
1926 int t,
1927 IMeshArena *arena)
1928 {
1929 const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
1930 int t_in_cdt = -1;
1931 for (int i = 0; i < cd.input_face.size(); ++i) {
1932 if (cd.input_face[i] == t) {
1933 t_in_cdt = i;
1934 }
1935 }
1936 if (t_in_cdt == -1) {
1937 std::cout << "Could not find " << t << " in cdt input tris\n";
1938 BLI_assert(false);
1939 return IMesh();
1940 }
1941 int t_orig = in_tm.face(t)->orig;
1942 constexpr int inline_buf_size = 20;
1943 Vector<Face *, inline_buf_size> faces;
1944 for (int f : cdt_out.face.index_range()) {
1945 if (cdt_out.face_orig[f].contains(t_in_cdt)) {
1946 BLI_assert(cdt_out.face[f].size() == 3);
1947 int i0 = cdt_out.face[f][0];
1948 int i1 = cdt_out.face[f][1];
1949 int i2 = cdt_out.face[f][2];
1950 mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]);
1951 mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]);
1952 mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]);
1953 /* No need to provide an original index: if coord matches
1954 * an original one, then it will already be in the arena
1955 * with the correct orig field. */
1956 const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX);
1957 const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX);
1958 const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX);
1959 Face *facep;
1960 bool is_isect0;
1961 bool is_isect1;
1962 bool is_isect2;
1963 if (cd.is_reversed[t_in_cdt]) {
1964 int oe0 = get_cdt_edge_orig(i0, i2, cd, in_tm, &is_isect0);
1965 int oe1 = get_cdt_edge_orig(i2, i1, cd, in_tm, &is_isect1);
1966 int oe2 = get_cdt_edge_orig(i1, i0, cd, in_tm, &is_isect2);
1967 facep = arena->add_face(
1968 {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
1969 }
1970 else {
1971 int oe0 = get_cdt_edge_orig(i0, i1, cd, in_tm, &is_isect0);
1972 int oe1 = get_cdt_edge_orig(i1, i2, cd, in_tm, &is_isect1);
1973 int oe2 = get_cdt_edge_orig(i2, i0, cd, in_tm, &is_isect2);
1974 facep = arena->add_face(
1975 {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
1976 }
1977 facep->populate_plane(false);
1978 faces.append(facep);
1979 }
1980 }
1981 return IMesh(faces);
1982 }
1983
extract_single_tri(const IMesh & tm,int t)1984 static IMesh extract_single_tri(const IMesh &tm, int t)
1985 {
1986 Face *f = tm.face(t);
1987 return IMesh({f});
1988 }
1989
bvhtreeverlap_cmp(const BVHTreeOverlap & a,const BVHTreeOverlap & b)1990 static bool bvhtreeverlap_cmp(const BVHTreeOverlap &a, const BVHTreeOverlap &b)
1991 {
1992 if (a.indexA < b.indexA) {
1993 return true;
1994 }
1995 if ((a.indexA == b.indexA) & (a.indexB < b.indexB)) {
1996 return true;
1997 }
1998 return false;
1999 }
2000 class TriOverlaps {
2001 BVHTree *tree_{nullptr};
2002 BVHTree *tree_b_{nullptr};
2003 BVHTreeOverlap *overlap_{nullptr};
2004 Array<int> first_overlap_;
2005 uint overlap_tot_{0};
2006
2007 struct CBData {
2008 const IMesh &tm;
2009 std::function<int(int)> shape_fn;
2010 int nshapes;
2011 bool use_self;
2012 };
2013
2014 public:
TriOverlaps(const IMesh & tm,const Array<BoundingBox> & tri_bb,int nshapes,std::function<int (int)> shape_fn,bool use_self)2015 TriOverlaps(const IMesh &tm,
2016 const Array<BoundingBox> &tri_bb,
2017 int nshapes,
2018 std::function<int(int)> shape_fn,
2019 bool use_self)
2020 {
2021 constexpr int dbg_level = 0;
2022 if (dbg_level > 0) {
2023 std::cout << "TriOverlaps construction\n";
2024 }
2025 /* Tree type is 8 => octtree; axis = 6 => using XYZ axes only. */
2026 tree_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6);
2027 /* In the common case of a binary boolean and no self intersection in
2028 * each shape, we will use two trees and simple bounding box overlap. */
2029 bool two_trees_no_self = nshapes == 2 && !use_self;
2030 if (two_trees_no_self) {
2031 tree_b_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6);
2032 }
2033 float bbpts[6];
2034 for (int t : tm.face_index_range()) {
2035 const BoundingBox &bb = tri_bb[t];
2036 copy_v3_v3(bbpts, bb.min);
2037 copy_v3_v3(bbpts + 3, bb.max);
2038 int shape = shape_fn(tm.face(t)->orig);
2039 if (two_trees_no_self) {
2040 if (shape == 0) {
2041 BLI_bvhtree_insert(tree_, t, bbpts, 2);
2042 }
2043 else if (shape == 1) {
2044 BLI_bvhtree_insert(tree_b_, t, bbpts, 2);
2045 }
2046 }
2047 else {
2048 if (shape != -1) {
2049 BLI_bvhtree_insert(tree_, t, bbpts, 2);
2050 }
2051 }
2052 }
2053 BLI_bvhtree_balance(tree_);
2054 if (two_trees_no_self) {
2055 BLI_bvhtree_balance(tree_b_);
2056 /* Don't expect a lot of trivial intersects in this case. */
2057 overlap_ = BLI_bvhtree_overlap(tree_, tree_b_, &overlap_tot_, NULL, NULL);
2058 }
2059 else {
2060 CBData cbdata{tm, shape_fn, nshapes, use_self};
2061 if (nshapes == 1) {
2062 overlap_ = BLI_bvhtree_overlap(tree_, tree_, &overlap_tot_, NULL, NULL);
2063 }
2064 else {
2065 overlap_ = BLI_bvhtree_overlap(
2066 tree_, tree_, &overlap_tot_, only_different_shapes, &cbdata);
2067 }
2068 }
2069 /* The rest of the code is simpler and easier to parallelize if, in the two-trees case,
2070 * we repeat the overlaps with indexA and indexB reversed. It is important that
2071 * in the repeated part, sorting will then bring things with indexB together. */
2072 if (two_trees_no_self) {
2073 overlap_ = static_cast<BVHTreeOverlap *>(
2074 MEM_reallocN(overlap_, 2 * overlap_tot_ * sizeof(overlap_[0])));
2075 for (uint i = 0; i < overlap_tot_; ++i) {
2076 overlap_[overlap_tot_ + i].indexA = overlap_[i].indexB;
2077 overlap_[overlap_tot_ + i].indexB = overlap_[i].indexA;
2078 }
2079 overlap_tot_ += overlap_tot_;
2080 }
2081 /* Sort the overlaps to bring all the intersects with a given indexA together. */
2082 std::sort(overlap_, overlap_ + overlap_tot_, bvhtreeverlap_cmp);
2083 if (dbg_level > 0) {
2084 std::cout << overlap_tot_ << " overlaps found:\n";
2085 for (BVHTreeOverlap ov : overlap()) {
2086 std::cout << "A: " << ov.indexA << ", B: " << ov.indexB << "\n";
2087 }
2088 }
2089 first_overlap_ = Array<int>(tm.face_size(), -1);
2090 for (int i = 0; i < static_cast<int>(overlap_tot_); ++i) {
2091 int t = overlap_[i].indexA;
2092 if (first_overlap_[t] == -1) {
2093 first_overlap_[t] = i;
2094 }
2095 }
2096 }
2097
~TriOverlaps()2098 ~TriOverlaps()
2099 {
2100 if (tree_) {
2101 BLI_bvhtree_free(tree_);
2102 }
2103 if (tree_b_) {
2104 BLI_bvhtree_free(tree_b_);
2105 }
2106 if (overlap_) {
2107 MEM_freeN(overlap_);
2108 }
2109 }
2110
overlap() const2111 Span<BVHTreeOverlap> overlap() const
2112 {
2113 return Span<BVHTreeOverlap>(overlap_, overlap_tot_);
2114 }
2115
first_overlap_index(int t) const2116 int first_overlap_index(int t) const
2117 {
2118 return first_overlap_[t];
2119 }
2120
2121 private:
only_different_shapes(void * userdata,int index_a,int index_b,int UNUSED (thread))2122 static bool only_different_shapes(void *userdata, int index_a, int index_b, int UNUSED(thread))
2123 {
2124 CBData *cbdata = static_cast<CBData *>(userdata);
2125 return cbdata->tm.face(index_a)->orig != cbdata->tm.face(index_b)->orig;
2126 }
2127 };
2128
2129 /**
2130 * Data needed for parallelization of #calc_overlap_itts.
2131 */
2132 struct OverlapIttsData {
2133 Vector<std::pair<int, int>> intersect_pairs;
2134 Map<std::pair<int, int>, ITT_value> &itt_map;
2135 const IMesh &tm;
2136 IMeshArena *arena;
2137
OverlapIttsDatablender::meshintersect::OverlapIttsData2138 OverlapIttsData(Map<std::pair<int, int>, ITT_value> &itt_map, const IMesh &tm, IMeshArena *arena)
2139 : itt_map(itt_map), tm(tm), arena(arena)
2140 {
2141 }
2142 };
2143
2144 /**
2145 * Return a std::pair containing a and b in canonical order:
2146 * With a <= b.
2147 */
canon_int_pair(int a,int b)2148 static std::pair<int, int> canon_int_pair(int a, int b)
2149 {
2150 if (a > b) {
2151 std::swap(a, b);
2152 }
2153 return std::pair<int, int>(a, b);
2154 }
2155
calc_overlap_itts_range_func(void * __restrict userdata,const int iter,const TaskParallelTLS * __restrict UNUSED (tls))2156 static void calc_overlap_itts_range_func(void *__restrict userdata,
2157 const int iter,
2158 const TaskParallelTLS *__restrict UNUSED(tls))
2159 {
2160 constexpr int dbg_level = 0;
2161 OverlapIttsData *data = static_cast<OverlapIttsData *>(userdata);
2162 std::pair<int, int> tri_pair = data->intersect_pairs[iter];
2163 int a = tri_pair.first;
2164 int b = tri_pair.second;
2165 if (dbg_level > 0) {
2166 std::cout << "calc_overlap_itts_range_func a=" << a << ", b=" << b << "\n";
2167 }
2168 ITT_value itt = intersect_tri_tri(data->tm, a, b);
2169 if (dbg_level > 0) {
2170 std::cout << "result of intersecting " << a << " and " << b << " = " << itt << "\n";
2171 }
2172 BLI_assert(data->itt_map.contains(tri_pair));
2173 data->itt_map.add_overwrite(tri_pair, itt);
2174 }
2175
2176 /**
2177 * Fill in itt_map with the vector of ITT_values that result from intersecting the triangles in ov.
2178 * Use a canonical order for triangles: (a,b) where a < b.
2179 */
calc_overlap_itts(Map<std::pair<int,int>,ITT_value> & itt_map,const IMesh & tm,const TriOverlaps & ov,IMeshArena * arena)2180 static void calc_overlap_itts(Map<std::pair<int, int>, ITT_value> &itt_map,
2181 const IMesh &tm,
2182 const TriOverlaps &ov,
2183 IMeshArena *arena)
2184 {
2185 OverlapIttsData data(itt_map, tm, arena);
2186 /* Put dummy values in `itt_map` initially,
2187 * so map entries will exist when doing the range function.
2188 * This means we won't have to protect the `itt_map.add_overwrite` function with a lock. */
2189 for (const BVHTreeOverlap &olap : ov.overlap()) {
2190 std::pair<int, int> key = canon_int_pair(olap.indexA, olap.indexB);
2191 if (!itt_map.contains(key)) {
2192 itt_map.add_new(key, ITT_value());
2193 data.intersect_pairs.append(key);
2194 }
2195 }
2196 int tot_intersect_pairs = data.intersect_pairs.size();
2197 TaskParallelSettings settings;
2198 BLI_parallel_range_settings_defaults(&settings);
2199 settings.min_iter_per_thread = 1000;
2200 settings.use_threading = intersect_use_threading;
2201 BLI_task_parallel_range(0, tot_intersect_pairs, &data, calc_overlap_itts_range_func, &settings);
2202 }
2203
2204 /**
2205 * Data needed for parallelization of calc_subdivided_tris.
2206 */
2207 struct OverlapTriRange {
2208 int tri_index;
2209 int overlap_start;
2210 int len;
2211 };
2212 struct SubdivideTrisData {
2213 Array<IMesh> &r_tri_subdivided;
2214 const IMesh &tm;
2215 const Map<std::pair<int, int>, ITT_value> &itt_map;
2216 Span<BVHTreeOverlap> overlap;
2217 IMeshArena *arena;
2218
2219 /* This vector gives, for each triangle in tm that has an intersection
2220 * we want to calculate: what the index of that triangle in tm is,
2221 * where it starts in the ov structure as indexA, and how many
2222 * overlap pairs have that same indexA (they will be continuous). */
2223 Vector<OverlapTriRange> overlap_tri_range;
2224
SubdivideTrisDatablender::meshintersect::SubdivideTrisData2225 SubdivideTrisData(Array<IMesh> &r_tri_subdivided,
2226 const IMesh &tm,
2227 const Map<std::pair<int, int>, ITT_value> &itt_map,
2228 Span<BVHTreeOverlap> overlap,
2229 IMeshArena *arena)
2230 : r_tri_subdivided(r_tri_subdivided),
2231 tm(tm),
2232 itt_map(itt_map),
2233 overlap(overlap),
2234 arena(arena),
2235 overlap_tri_range{}
2236 {
2237 }
2238 };
2239
calc_subdivided_tri_range_func(void * __restrict userdata,const int iter,const TaskParallelTLS * __restrict UNUSED (tls))2240 static void calc_subdivided_tri_range_func(void *__restrict userdata,
2241 const int iter,
2242 const TaskParallelTLS *__restrict UNUSED(tls))
2243 {
2244 constexpr int dbg_level = 0;
2245 SubdivideTrisData *data = static_cast<SubdivideTrisData *>(userdata);
2246 OverlapTriRange &otr = data->overlap_tri_range[iter];
2247 int t = otr.tri_index;
2248 if (dbg_level > 0) {
2249 std::cout << "calc_subdivided_tri_range_func\nt=" << t << " start=" << otr.overlap_start
2250 << " len=" << otr.len << "\n";
2251 }
2252 constexpr int inline_capacity = 100;
2253 Vector<ITT_value, inline_capacity> itts(otr.len);
2254 for (int j = otr.overlap_start; j < otr.overlap_start + otr.len; ++j) {
2255 int t_other = data->overlap[j].indexB;
2256 std::pair<int, int> key = canon_int_pair(t, t_other);
2257 ITT_value itt;
2258 if (data->itt_map.contains(key)) {
2259 itt = data->itt_map.lookup(key);
2260 }
2261 if (itt.kind != INONE) {
2262 itts.append(itt);
2263 }
2264 if (dbg_level > 0) {
2265 std::cout << " tri t" << t_other << "; result = " << itt << "\n";
2266 }
2267 }
2268 if (itts.size() > 0) {
2269 CDT_data cd_data = prepare_cdt_input(data->tm, t, itts);
2270 do_cdt(cd_data);
2271 data->r_tri_subdivided[t] = extract_subdivided_tri(cd_data, data->tm, t, data->arena);
2272 if (dbg_level > 0) {
2273 std::cout << "subdivide output\n" << data->r_tri_subdivided[t];
2274 }
2275 }
2276 }
2277
2278 /**
2279 * For each triangle in tm, fill in the corresponding slot in
2280 * r_tri_subdivided with the result of intersecting it with
2281 * all the other triangles in the mesh, if it intersects any others.
2282 * But don't do this for triangles that are part of a cluster.
2283 * Also, do nothing here if the answer is just the triangle itself.
2284 */
calc_subdivided_tris(Array<IMesh> & r_tri_subdivided,const IMesh & tm,const Map<std::pair<int,int>,ITT_value> & itt_map,const CoplanarClusterInfo & clinfo,const TriOverlaps & ov,IMeshArena * arena)2285 static void calc_subdivided_tris(Array<IMesh> &r_tri_subdivided,
2286 const IMesh &tm,
2287 const Map<std::pair<int, int>, ITT_value> &itt_map,
2288 const CoplanarClusterInfo &clinfo,
2289 const TriOverlaps &ov,
2290 IMeshArena *arena)
2291 {
2292 const int dbg_level = 0;
2293 if (dbg_level > 0) {
2294 std::cout << "\nCALC_SUBDIVIDED_TRIS\n\n";
2295 }
2296 Span<BVHTreeOverlap> overlap = ov.overlap();
2297 SubdivideTrisData data(r_tri_subdivided, tm, itt_map, overlap, arena);
2298 int overlap_tot = overlap.size();
2299 data.overlap_tri_range = Vector<OverlapTriRange>();
2300 data.overlap_tri_range.reserve(overlap_tot);
2301 int overlap_index = 0;
2302 while (overlap_index < overlap_tot) {
2303 int t = overlap[overlap_index].indexA;
2304 int i = overlap_index;
2305 while (i + 1 < overlap_tot && overlap[i + 1].indexA == t) {
2306 ++i;
2307 }
2308 /* Now overlap[overlap_index] to overlap[i] have indexA == t.
2309 * We only record ranges for triangles that are not in clusters,
2310 * because the ones in clusters are handled separately. */
2311 if (clinfo.tri_cluster(t) == NO_INDEX) {
2312 int len = i - overlap_index + 1;
2313 if (!(len == 1 && overlap[overlap_index].indexB == t)) {
2314 OverlapTriRange range = {t, overlap_index, len};
2315 data.overlap_tri_range.append(range);
2316 # ifdef PERFDEBUG
2317 bumpperfcount(0, len); /* Non-cluster overlaps. */
2318 # endif
2319 }
2320 }
2321 overlap_index = i + 1;
2322 }
2323 int overlap_tri_range_tot = data.overlap_tri_range.size();
2324 TaskParallelSettings settings;
2325 BLI_parallel_range_settings_defaults(&settings);
2326 settings.min_iter_per_thread = 50;
2327 settings.use_threading = intersect_use_threading;
2328 BLI_task_parallel_range(
2329 0, overlap_tri_range_tot, &data, calc_subdivided_tri_range_func, &settings);
2330 }
2331
calc_cluster_subdivided(const CoplanarClusterInfo & clinfo,int c,const IMesh & tm,const TriOverlaps & ov,const Map<std::pair<int,int>,ITT_value> & itt_map,IMeshArena * UNUSED (arena))2332 static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo,
2333 int c,
2334 const IMesh &tm,
2335 const TriOverlaps &ov,
2336 const Map<std::pair<int, int>, ITT_value> &itt_map,
2337 IMeshArena *UNUSED(arena))
2338 {
2339 constexpr int dbg_level = 0;
2340 BLI_assert(c < clinfo.tot_cluster());
2341 const CoplanarCluster &cl = clinfo.cluster(c);
2342 /* Make a CDT input with triangles from C and intersects from other triangles in tm. */
2343 if (dbg_level > 0) {
2344 std::cout << "CALC_CLUSTER_SUBDIVIDED for cluster " << c << " = " << cl << "\n";
2345 }
2346 /* Get vector itts of all intersections of a triangle of cl with any triangle of tm not
2347 * in cl and not co-planar with it (for that latter, if there were an intersection,
2348 * it should already be in cluster cl). */
2349 Vector<ITT_value> itts;
2350 Span<BVHTreeOverlap> ovspan = ov.overlap();
2351 for (int t : cl) {
2352 if (dbg_level > 0) {
2353 std::cout << "find intersects with triangle " << t << " of cluster\n";
2354 }
2355 int first_i = ov.first_overlap_index(t);
2356 if (first_i == -1) {
2357 continue;
2358 }
2359 for (int i = first_i; i < ovspan.size() && ovspan[i].indexA == t; ++i) {
2360 int t_other = ovspan[i].indexB;
2361 if (clinfo.tri_cluster(t_other) != c) {
2362 if (dbg_level > 0) {
2363 std::cout << "use intersect(" << t << "," << t_other << "\n";
2364 }
2365 std::pair<int, int> key = canon_int_pair(t, t_other);
2366 if (itt_map.contains(key)) {
2367 ITT_value itt = itt_map.lookup(key);
2368 if (itt.kind != INONE && itt.kind != ICOPLANAR) {
2369 itts.append(itt);
2370 if (dbg_level > 0) {
2371 std::cout << " itt = " << itt << "\n";
2372 }
2373 }
2374 }
2375 }
2376 }
2377 }
2378 /* Use CDT to subdivide the cluster triangles and the points and segs in itts. */
2379 CDT_data cd_data = prepare_cdt_input_for_cluster(tm, clinfo, c, itts);
2380 do_cdt(cd_data);
2381 return cd_data;
2382 }
2383
union_tri_subdivides(const blender::Array<IMesh> & tri_subdivided)2384 static IMesh union_tri_subdivides(const blender::Array<IMesh> &tri_subdivided)
2385 {
2386 int tot_tri = 0;
2387 for (const IMesh &m : tri_subdivided) {
2388 tot_tri += m.face_size();
2389 }
2390 Array<Face *> faces(tot_tri);
2391 int face_index = 0;
2392 for (const IMesh &m : tri_subdivided) {
2393 for (Face *f : m.faces()) {
2394 faces[face_index++] = f;
2395 }
2396 }
2397 return IMesh(faces);
2398 }
2399
find_clusters(const IMesh & tm,const Array<BoundingBox> & tri_bb,const Map<std::pair<int,int>,ITT_value> & itt_map)2400 static CoplanarClusterInfo find_clusters(const IMesh &tm,
2401 const Array<BoundingBox> &tri_bb,
2402 const Map<std::pair<int, int>, ITT_value> &itt_map)
2403 {
2404 constexpr int dbg_level = 0;
2405 if (dbg_level > 0) {
2406 std::cout << "FIND_CLUSTERS\n";
2407 }
2408 CoplanarClusterInfo ans(tm.face_size());
2409 /* Use a VectorSet to get stable order from run to run. */
2410 VectorSet<int> maybe_coplanar_tris;
2411 maybe_coplanar_tris.reserve(2 * itt_map.size());
2412 for (auto item : itt_map.items()) {
2413 if (item.value.kind == ICOPLANAR) {
2414 int t1 = item.key.first;
2415 int t2 = item.key.second;
2416 maybe_coplanar_tris.add_multiple({t1, t2});
2417 }
2418 }
2419 if (dbg_level > 0) {
2420 std::cout << "found " << maybe_coplanar_tris.size() << " possible coplanar tris\n";
2421 }
2422 if (maybe_coplanar_tris.size() == 0) {
2423 if (dbg_level > 0) {
2424 std::cout << "No possible coplanar tris, so no clusters\n";
2425 }
2426 return ans;
2427 }
2428 /* There can be more than one #CoplanarCluster per plane. Accumulate them in
2429 * a Vector. We will have to merge some elements of the Vector as we discover
2430 * triangles that form intersection bridges between two or more clusters. */
2431 Map<Plane, Vector<CoplanarCluster>> plane_cls;
2432 plane_cls.reserve(maybe_coplanar_tris.size());
2433 for (int t : maybe_coplanar_tris) {
2434 /* Use a canonical version of the plane for map index.
2435 * We can't just store the canonical version in the face
2436 * since canonicalizing loses the orientation of the normal. */
2437 Plane tplane = *tm.face(t)->plane;
2438 BLI_assert(tplane.exact_populated());
2439 tplane.make_canonical();
2440 if (dbg_level > 0) {
2441 std::cout << "plane for tri " << t << " = " << &tplane << "\n";
2442 }
2443 /* Assume all planes are in canonical from (see canon_plane()). */
2444 if (plane_cls.contains(tplane)) {
2445 Vector<CoplanarCluster> &curcls = plane_cls.lookup(tplane);
2446 if (dbg_level > 0) {
2447 std::cout << "already has " << curcls.size() << " clusters\n";
2448 }
2449 /* Partition `curcls` into those that intersect t non-trivially, and those that don't. */
2450 Vector<CoplanarCluster *> int_cls;
2451 Vector<CoplanarCluster *> no_int_cls;
2452 for (CoplanarCluster &cl : curcls) {
2453 if (dbg_level > 1) {
2454 std::cout << "consider intersecting with cluster " << cl << "\n";
2455 }
2456 if (bbs_might_intersect(tri_bb[t], cl.bounding_box())) {
2457 if (dbg_level > 1) {
2458 std::cout << "append to int_cls\n";
2459 }
2460 int_cls.append(&cl);
2461 }
2462 else {
2463 if (dbg_level > 1) {
2464 std::cout << "append to no_int_cls\n";
2465 }
2466 no_int_cls.append(&cl);
2467 }
2468 }
2469 if (int_cls.size() == 0) {
2470 /* t doesn't intersect any existing cluster in its plane, so make one just for it. */
2471 if (dbg_level > 1) {
2472 std::cout << "no intersecting clusters for t, make a new one\n";
2473 }
2474 curcls.append(CoplanarCluster(t, tri_bb[t]));
2475 }
2476 else if (int_cls.size() == 1) {
2477 /* t intersects exactly one existing cluster, so can add t to that cluster. */
2478 if (dbg_level > 1) {
2479 std::cout << "exactly one existing cluster, " << int_cls[0] << ", adding to it\n";
2480 }
2481 int_cls[0]->add_tri(t, tri_bb[t]);
2482 }
2483 else {
2484 /* t intersections 2 or more existing clusters: need to merge them and replace all the
2485 * originals with the merged one in `curcls`. */
2486 if (dbg_level > 1) {
2487 std::cout << "merging\n";
2488 }
2489 CoplanarCluster mergecl;
2490 mergecl.add_tri(t, tri_bb[t]);
2491 for (CoplanarCluster *cl : int_cls) {
2492 for (int t : *cl) {
2493 mergecl.add_tri(t, tri_bb[t]);
2494 }
2495 }
2496 Vector<CoplanarCluster> newvec;
2497 newvec.append(mergecl);
2498 for (CoplanarCluster *cl_no_int : no_int_cls) {
2499 newvec.append(*cl_no_int);
2500 }
2501 plane_cls.add_overwrite(tplane, newvec);
2502 }
2503 }
2504 else {
2505 if (dbg_level > 0) {
2506 std::cout << "first cluster for its plane\n";
2507 }
2508 plane_cls.add_new(tplane, Vector<CoplanarCluster>{CoplanarCluster(t, tri_bb[t])});
2509 }
2510 }
2511 /* Does this give deterministic order for cluster ids? I think so, since
2512 * hash for planes is on their values, not their addresses. */
2513 for (auto item : plane_cls.items()) {
2514 for (const CoplanarCluster &cl : item.value) {
2515 if (cl.tot_tri() > 1) {
2516 ans.add_cluster(cl);
2517 }
2518 }
2519 }
2520
2521 return ans;
2522 }
2523
face_is_degenerate(const Face * f)2524 static bool face_is_degenerate(const Face *f)
2525 {
2526 const Face &face = *f;
2527 const Vert *v0 = face[0];
2528 const Vert *v1 = face[1];
2529 const Vert *v2 = face[2];
2530 if (v0 == v1 || v0 == v2 || v1 == v2) {
2531 return true;
2532 }
2533 double3 da = v2->co - v0->co;
2534 double3 db = v2->co - v1->co;
2535 double3 dab = double3::cross_high_precision(da, db);
2536 double dab_length_squared = dab.length_squared();
2537 double err_bound = supremum_dot_cross(dab, dab) * index_dot_cross * DBL_EPSILON;
2538 if (dab_length_squared > err_bound) {
2539 return false;
2540 }
2541 mpq3 a = v2->co_exact - v0->co_exact;
2542 mpq3 b = v2->co_exact - v1->co_exact;
2543 mpq3 ab = mpq3::cross(a, b);
2544 if (ab.x == 0 && ab.y == 0 && ab.z == 0) {
2545 return true;
2546 }
2547
2548 return false;
2549 }
2550
2551 /* Data and functions to test triangle degeneracy in parallel. */
2552 struct DegenData {
2553 const IMesh &tm;
2554 };
2555
2556 struct DegenChunkData {
2557 bool has_degenerate_tri = false;
2558 };
2559
degenerate_range_func(void * __restrict userdata,const int iter,const TaskParallelTLS * __restrict tls)2560 static void degenerate_range_func(void *__restrict userdata,
2561 const int iter,
2562 const TaskParallelTLS *__restrict tls)
2563 {
2564 DegenData *data = static_cast<DegenData *>(userdata);
2565 DegenChunkData *chunk_data = static_cast<DegenChunkData *>(tls->userdata_chunk);
2566 const Face *f = data->tm.face(iter);
2567 bool is_degenerate = face_is_degenerate(f);
2568 chunk_data->has_degenerate_tri |= is_degenerate;
2569 }
2570
degenerate_reduce(const void * __restrict UNUSED (userdata),void * __restrict chunk_join,void * __restrict chunk)2571 static void degenerate_reduce(const void *__restrict UNUSED(userdata),
2572 void *__restrict chunk_join,
2573 void *__restrict chunk)
2574 {
2575 DegenChunkData *degen_chunk_join = static_cast<DegenChunkData *>(chunk_join);
2576 DegenChunkData *degen_chunk = static_cast<DegenChunkData *>(chunk);
2577 degen_chunk_join->has_degenerate_tri |= degen_chunk->has_degenerate_tri;
2578 }
2579
2580 /* Does triangle #IMesh tm have any triangles with zero area? */
has_degenerate_tris(const IMesh & tm)2581 static bool has_degenerate_tris(const IMesh &tm)
2582 {
2583 DegenData degen_data = {tm};
2584 DegenChunkData degen_chunk_data;
2585 TaskParallelSettings settings;
2586 BLI_parallel_range_settings_defaults(&settings);
2587 settings.userdata_chunk = °en_chunk_data;
2588 settings.userdata_chunk_size = sizeof(degen_chunk_data);
2589 settings.func_reduce = degenerate_reduce;
2590 settings.min_iter_per_thread = 1000;
2591 settings.use_threading = intersect_use_threading;
2592 BLI_task_parallel_range(0, tm.face_size(), °en_data, degenerate_range_func, &settings);
2593 return degen_chunk_data.has_degenerate_tri;
2594 }
2595
remove_degenerate_tris(const IMesh & tm_in)2596 static IMesh remove_degenerate_tris(const IMesh &tm_in)
2597 {
2598 IMesh ans;
2599 Vector<Face *> new_faces;
2600 new_faces.reserve(tm_in.face_size());
2601 for (Face *f : tm_in.faces()) {
2602 if (!face_is_degenerate(f)) {
2603 new_faces.append(f);
2604 }
2605 }
2606 ans.set_faces(new_faces);
2607 return ans;
2608 }
2609
2610 /* This is the main routine for calculating the self_intersection of a triangle mesh. */
trimesh_self_intersect(const IMesh & tm_in,IMeshArena * arena)2611 IMesh trimesh_self_intersect(const IMesh &tm_in, IMeshArena *arena)
2612 {
2613 return trimesh_nary_intersect(
2614 tm_in, 1, [](int UNUSED(t)) { return 0; }, true, arena);
2615 }
2616
trimesh_nary_intersect(const IMesh & tm_in,int nshapes,std::function<int (int)> shape_fn,bool use_self,IMeshArena * arena)2617 IMesh trimesh_nary_intersect(const IMesh &tm_in,
2618 int nshapes,
2619 std::function<int(int)> shape_fn,
2620 bool use_self,
2621 IMeshArena *arena)
2622 {
2623 constexpr int dbg_level = 0;
2624 if (dbg_level > 0) {
2625 std::cout << "\nTRIMESH_NARY_INTERSECT nshapes=" << nshapes << " use_self=" << use_self
2626 << "\n";
2627 for (const Face *f : tm_in.faces()) {
2628 BLI_assert(f->is_tri());
2629 UNUSED_VARS_NDEBUG(f);
2630 }
2631 if (dbg_level > 1) {
2632 std::cout << "input mesh:\n" << tm_in;
2633 for (int t : tm_in.face_index_range()) {
2634 std::cout << "shape(" << t << ") = " << shape_fn(tm_in.face(t)->orig) << "\n";
2635 }
2636 write_obj_mesh(const_cast<IMesh &>(tm_in), "trimesh_input");
2637 }
2638 }
2639 # ifdef PERFDEBUG
2640 perfdata_init();
2641 double start_time = PIL_check_seconds_timer();
2642 std::cout << "trimesh_nary_intersect start\n";
2643 # endif
2644 /* Usually can use tm_in but if it has degenerate or illegal triangles,
2645 * then need to work on a copy of it without those triangles. */
2646 const IMesh *tm_clean = &tm_in;
2647 IMesh tm_cleaned;
2648 if (has_degenerate_tris(tm_in)) {
2649 if (dbg_level > 0) {
2650 std::cout << "cleaning degenerate triangles\n";
2651 }
2652 tm_cleaned = remove_degenerate_tris(tm_in);
2653 tm_clean = &tm_cleaned;
2654 if (dbg_level > 1) {
2655 std::cout << "cleaned input mesh:\n" << tm_cleaned;
2656 }
2657 }
2658 # ifdef PERFDEBUG
2659 double clean_time = PIL_check_seconds_timer();
2660 std::cout << "cleaned, time = " << clean_time - start_time << "\n";
2661 # endif
2662 Array<BoundingBox> tri_bb = calc_face_bounding_boxes(*tm_clean);
2663 # ifdef PERFDEBUG
2664 double bb_calc_time = PIL_check_seconds_timer();
2665 std::cout << "bbs calculated, time = " << bb_calc_time - clean_time << "\n";
2666 # endif
2667 TriOverlaps tri_ov(*tm_clean, tri_bb, nshapes, shape_fn, use_self);
2668 # ifdef PERFDEBUG
2669 double overlap_time = PIL_check_seconds_timer();
2670 std::cout << "intersect overlaps calculated, time = " << overlap_time - bb_calc_time << "\n";
2671 # endif
2672 for (int t : tm_clean->face_index_range()) {
2673 if (tri_ov.first_overlap_index(t) != -1) {
2674 tm_clean->face(t)->populate_plane(true);
2675 }
2676 }
2677 # ifdef PERFDEBUG
2678 double plane_populate = PIL_check_seconds_timer();
2679 std::cout << "planes populated, time = " << plane_populate - overlap_time << "\n";
2680 # endif
2681 /* itt_map((a,b)) will hold the intersection value resulting from intersecting
2682 * triangles with indices a and b, where a < b. */
2683 Map<std::pair<int, int>, ITT_value> itt_map;
2684 itt_map.reserve(tri_ov.overlap().size());
2685 calc_overlap_itts(itt_map, *tm_clean, tri_ov, arena);
2686 # ifdef PERFDEBUG
2687 double itt_time = PIL_check_seconds_timer();
2688 std::cout << "itts found, time = " << itt_time - plane_populate << "\n";
2689 # endif
2690 CoplanarClusterInfo clinfo = find_clusters(*tm_clean, tri_bb, itt_map);
2691 if (dbg_level > 1) {
2692 std::cout << clinfo;
2693 }
2694 # ifdef PERFDEBUG
2695 double find_cluster_time = PIL_check_seconds_timer();
2696 std::cout << "clusters found, time = " << find_cluster_time - itt_time << "\n";
2697 doperfmax(0, tm_in.face_size());
2698 doperfmax(1, clinfo.tot_cluster());
2699 doperfmax(2, tri_ov.overlap().size());
2700 # endif
2701 Array<IMesh> tri_subdivided(tm_clean->face_size());
2702 calc_subdivided_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena);
2703 # ifdef PERFDEBUG
2704 double subdivided_tris_time = PIL_check_seconds_timer();
2705 std::cout << "subdivided tris found, time = " << subdivided_tris_time - itt_time << "\n";
2706 # endif
2707 Array<CDT_data> cluster_subdivided(clinfo.tot_cluster());
2708 for (int c : clinfo.index_range()) {
2709 cluster_subdivided[c] = calc_cluster_subdivided(clinfo, c, *tm_clean, tri_ov, itt_map, arena);
2710 }
2711 # ifdef PERFDEBUG
2712 double cluster_subdivide_time = PIL_check_seconds_timer();
2713 std::cout << "subdivided clusters found, time = "
2714 << cluster_subdivide_time - subdivided_tris_time << "\n";
2715 # endif
2716 for (int t : tm_clean->face_index_range()) {
2717 int c = clinfo.tri_cluster(t);
2718 if (c != NO_INDEX) {
2719 BLI_assert(tri_subdivided[t].face_size() == 0);
2720 tri_subdivided[t] = extract_subdivided_tri(cluster_subdivided[c], *tm_clean, t, arena);
2721 }
2722 else if (tri_subdivided[t].face_size() == 0) {
2723 tri_subdivided[t] = extract_single_tri(*tm_clean, t);
2724 }
2725 }
2726 # ifdef PERFDEBUG
2727 double extract_time = PIL_check_seconds_timer();
2728 std::cout << "triangles extracted, time = " << extract_time - cluster_subdivide_time << "\n";
2729 # endif
2730 IMesh combined = union_tri_subdivides(tri_subdivided);
2731 if (dbg_level > 1) {
2732 std::cout << "TRIMESH_NARY_INTERSECT answer:\n";
2733 std::cout << combined;
2734 }
2735 # ifdef PERFDEBUG
2736 double end_time = PIL_check_seconds_timer();
2737 std::cout << "triangles combined, time = " << end_time - extract_time << "\n";
2738 std::cout << "trimesh_nary_intersect done, total time = " << end_time - start_time << "\n";
2739 dump_perfdata();
2740 # endif
2741 return combined;
2742 }
2743
operator <<(std::ostream & os,const CoplanarCluster & cl)2744 static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl)
2745 {
2746 os << "cl(";
2747 bool first = true;
2748 for (const int t : cl) {
2749 if (first) {
2750 first = false;
2751 }
2752 else {
2753 os << ",";
2754 }
2755 os << t;
2756 }
2757 os << ")";
2758 return os;
2759 }
2760
operator <<(std::ostream & os,const CoplanarClusterInfo & clinfo)2761 static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo)
2762 {
2763 os << "Coplanar Cluster Info:\n";
2764 for (int c : clinfo.index_range()) {
2765 os << c << ": " << clinfo.cluster(c) << "\n";
2766 }
2767 return os;
2768 }
2769
operator <<(std::ostream & os,const ITT_value & itt)2770 static std::ostream &operator<<(std::ostream &os, const ITT_value &itt)
2771 {
2772 switch (itt.kind) {
2773 case INONE:
2774 os << "none";
2775 break;
2776 case IPOINT:
2777 os << "point " << itt.p1;
2778 break;
2779 case ISEGMENT:
2780 os << "segment " << itt.p1 << " " << itt.p2;
2781 break;
2782 case ICOPLANAR:
2783 os << "co-planar t" << itt.t_source;
2784 break;
2785 }
2786 return os;
2787 }
2788
2789 /**
2790 * Writing the obj_mesh has the side effect of populating verts.
2791 */
write_obj_mesh(IMesh & m,const std::string & objname)2792 void write_obj_mesh(IMesh &m, const std::string &objname)
2793 {
2794 /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library.
2795 * This is just for developer debugging anyway,
2796 * and should never be called in production Blender. */
2797 # ifdef _WIN_32
2798 const char *objdir = BLI_getenv("HOME");
2799 # else
2800 const char *objdir = "/tmp/";
2801 # endif
2802 if (m.face_size() == 0) {
2803 return;
2804 }
2805
2806 std::string fname = std::string(objdir) + objname + std::string(".obj");
2807 std::ofstream f;
2808 f.open(fname);
2809 if (!f) {
2810 std::cout << "Could not open file " << fname << "\n";
2811 return;
2812 }
2813
2814 if (!m.has_verts()) {
2815 m.populate_vert();
2816 }
2817 for (const Vert *v : m.vertices()) {
2818 const double3 dv = v->co;
2819 f << "v " << dv[0] << " " << dv[1] << " " << dv[2] << "\n";
2820 }
2821 int i = 0;
2822 for (const Face *face : m.faces()) {
2823 /* OBJ files use 1-indexing for vertices. */
2824 f << "f ";
2825 for (const Vert *v : *face) {
2826 int i = m.lookup_vert(v);
2827 BLI_assert(i != NO_INDEX);
2828 /* OBJ files use 1-indexing for vertices. */
2829 f << i + 1 << " ";
2830 }
2831 f << "\n";
2832 ++i;
2833 }
2834 f.close();
2835 }
2836
2837 # ifdef PERFDEBUG
2838 struct PerfCounts {
2839 Vector<int> count;
2840 Vector<const char *> count_name;
2841 Vector<int> max;
2842 Vector<const char *> max_name;
2843 };
2844
2845 static PerfCounts *perfdata = nullptr;
2846
perfdata_init(void)2847 static void perfdata_init(void)
2848 {
2849 perfdata = new PerfCounts;
2850
2851 /* count 0. */
2852 perfdata->count.append(0);
2853 perfdata->count_name.append("Non-cluster overlaps");
2854
2855 /* count 1. */
2856 perfdata->count.append(0);
2857 perfdata->count_name.append("intersect_tri_tri calls");
2858
2859 /* count 2. */
2860 perfdata->count.append(0);
2861 perfdata->count_name.append("tri tri intersects decided by filter plane tests");
2862
2863 /* count 3. */
2864 perfdata->count.append(0);
2865 perfdata->count_name.append("tri tri intersects decided by exact plane tests");
2866
2867 /* count 4. */
2868 perfdata->count.append(0);
2869 perfdata->count_name.append("final non-NONE intersects");
2870
2871 /* max 0. */
2872 perfdata->max.append(0);
2873 perfdata->max_name.append("total faces");
2874
2875 /* max 1. */
2876 perfdata->max.append(0);
2877 perfdata->max_name.append("total clusters");
2878
2879 /* max 2. */
2880 perfdata->max.append(0);
2881 perfdata->max_name.append("total overlaps");
2882 }
2883
incperfcount(int countnum)2884 static void incperfcount(int countnum)
2885 {
2886 perfdata->count[countnum]++;
2887 }
2888
bumpperfcount(int countnum,int amt)2889 static void bumpperfcount(int countnum, int amt)
2890 {
2891 perfdata->count[countnum] += amt;
2892 }
2893
doperfmax(int maxnum,int val)2894 static void doperfmax(int maxnum, int val)
2895 {
2896 perfdata->max[maxnum] = max_ii(perfdata->max[maxnum], val);
2897 }
2898
dump_perfdata(void)2899 static void dump_perfdata(void)
2900 {
2901 std::cout << "\nPERFDATA\n";
2902 for (int i : perfdata->count.index_range()) {
2903 std::cout << perfdata->count_name[i] << " = " << perfdata->count[i] << "\n";
2904 }
2905 for (int i : perfdata->max.index_range()) {
2906 std::cout << perfdata->max_name[i] << " = " << perfdata->max[i] << "\n";
2907 }
2908 delete perfdata;
2909 }
2910 # endif
2911
2912 } // namespace blender::meshintersect
2913
2914 #endif // WITH_GMP
2915