Blender  V3.3
mesh_intersect.cc
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 /* The #blender::meshintersect API needs GMP. */
8 #ifdef WITH_GMP
9 
10 # include <algorithm>
11 # include <fstream>
12 # include <iostream>
13 # include <memory>
14 
15 # include "BLI_allocator.hh"
16 # include "BLI_array.hh"
17 # include "BLI_assert.h"
18 # include "BLI_delaunay_2d.h"
19 # include "BLI_hash.hh"
20 # include "BLI_kdopbvh.h"
21 # include "BLI_map.hh"
22 # include "BLI_math_boolean.hh"
23 # include "BLI_math_mpq.hh"
24 # include "BLI_math_vec_mpq_types.hh"
25 # include "BLI_math_vec_types.hh"
26 # include "BLI_math_vector.h"
27 # include "BLI_polyfill_2d.h"
28 # include "BLI_set.hh"
29 # include "BLI_sort.hh"
30 # include "BLI_span.hh"
31 # include "BLI_task.h"
32 # include "BLI_task.hh"
33 # include "BLI_threads.h"
34 # include "BLI_vector.hh"
35 # include "BLI_vector_set.hh"
36 
37 # include "PIL_time.h"
38 
39 # include "BLI_mesh_intersect.hh"
40 
41 // # define PERFDEBUG
42 
43 namespace blender::meshintersect {
44 
45 # ifdef PERFDEBUG
46 static void perfdata_init();
47 static void incperfcount(int countnum);
48 static void bumpperfcount(int countnum, int amt);
49 static void doperfmax(int maxnum, int val);
50 static void dump_perfdata();
51 # endif
52 
54 static constexpr bool intersect_use_threading = true;
55 
56 Vert::Vert(const mpq3 &mco, const double3 &dco, int id, int orig)
57  : co_exact(mco), co(dco), id(id), orig(orig)
58 {
59 }
60 
61 bool Vert::operator==(const Vert &other) const
62 {
63  return this->co_exact == other.co_exact;
64 }
65 
66 uint64_t Vert::hash() const
67 {
68  uint64_t x = *reinterpret_cast<const uint64_t *>(&co.x);
69  uint64_t y = *reinterpret_cast<const uint64_t *>(&co.y);
70  uint64_t z = *reinterpret_cast<const uint64_t *>(&co.z);
71  x = (x >> 56) ^ (x >> 46) ^ x;
72  y = (y >> 55) ^ (y >> 45) ^ y;
73  z = (z >> 54) ^ (z >> 44) ^ z;
74  return x ^ y ^ z;
75 }
76 
77 std::ostream &operator<<(std::ostream &os, const Vert *v)
78 {
79  constexpr int dbg_level = 0;
80  os << "v" << v->id;
81  if (v->orig != NO_INDEX) {
82  os << "o" << v->orig;
83  }
84  os << v->co;
85  if (dbg_level > 0) {
86  os << "=" << v->co_exact;
87  }
88  return os;
89 }
90 
91 bool Plane::operator==(const Plane &other) const
92 {
93  return norm_exact == other.norm_exact && d_exact == other.d_exact;
94 }
95 
96 void Plane::make_canonical()
97 {
98  if (norm_exact[0] != 0) {
99  mpq_class den = norm_exact[0];
100  norm_exact = mpq3(1, norm_exact[1] / den, norm_exact[2] / den);
101  d_exact = d_exact / den;
102  }
103  else if (norm_exact[1] != 0) {
104  mpq_class den = norm_exact[1];
105  norm_exact = mpq3(0, 1, norm_exact[2] / den);
106  d_exact = d_exact / den;
107  }
108  else {
109  if (norm_exact[2] != 0) {
110  mpq_class den = norm_exact[2];
111  norm_exact = mpq3(0, 0, 1);
112  d_exact = d_exact / den;
113  }
114  else {
115  /* A degenerate plane. */
116  d_exact = 0;
117  }
118  }
119  norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d());
120  d = d_exact.get_d();
121 }
122 
123 Plane::Plane(const mpq3 &norm_exact, const mpq_class &d_exact)
124  : norm_exact(norm_exact), d_exact(d_exact)
125 {
126  norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d());
127  d = d_exact.get_d();
128 }
129 
130 Plane::Plane(const double3 &norm, const double d) : norm(norm), d(d)
131 {
132  norm_exact = mpq3(0, 0, 0); /* Marks as "exact not yet populated". */
133 }
134 
135 bool Plane::exact_populated() const
136 {
137  return norm_exact[0] != 0 || norm_exact[1] != 0 || norm_exact[2] != 0;
138 }
139 
140 uint64_t Plane::hash() const
141 {
142  uint64_t x = *reinterpret_cast<const uint64_t *>(&this->norm.x);
143  uint64_t y = *reinterpret_cast<const uint64_t *>(&this->norm.y);
144  uint64_t z = *reinterpret_cast<const uint64_t *>(&this->norm.z);
145  uint64_t d = *reinterpret_cast<const uint64_t *>(&this->d);
146  x = (x >> 56) ^ (x >> 46) ^ x;
147  y = (y >> 55) ^ (y >> 45) ^ y;
148  z = (z >> 54) ^ (z >> 44) ^ z;
149  d = (d >> 53) ^ (d >> 43) ^ d;
150  return x ^ y ^ z ^ d;
151 }
152 
153 std::ostream &operator<<(std::ostream &os, const Plane *plane)
154 {
155  os << "[" << plane->norm << ";" << plane->d << "]";
156  return os;
157 }
158 
159 Face::Face(
160  Span<const Vert *> verts, int id, int orig, Span<int> edge_origs, Span<bool> is_intersect)
161  : vert(verts), edge_orig(edge_origs), is_intersect(is_intersect), id(id), orig(orig)
162 {
163 }
164 
165 Face::Face(Span<const Vert *> verts, int id, int orig) : vert(verts), id(id), orig(orig)
166 {
167 }
168 
169 void Face::populate_plane(bool need_exact)
170 {
171  if (plane != nullptr) {
172  if (!need_exact || plane->exact_populated()) {
173  return;
174  }
175  }
176  if (need_exact) {
177  mpq3 normal_exact;
178  if (vert.size() > 3) {
179  Array<mpq3> co(vert.size());
180  for (int i : index_range()) {
181  co[i] = vert[i]->co_exact;
182  }
183  normal_exact = math::cross_poly(co.as_span());
184  }
185  else {
186  mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact;
187  mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact;
188  normal_exact = math::cross(tr02, tr12);
189  }
190  mpq_class d_exact = -math::dot(normal_exact, vert[0]->co_exact);
191  plane = new Plane(normal_exact, d_exact);
192  }
193  else {
194  double3 normal;
195  if (vert.size() > 3) {
196  Array<double3> co(vert.size());
197  for (int i : index_range()) {
198  co[i] = vert[i]->co;
199  }
200  normal = math::cross_poly(co.as_span());
201  }
202  else {
203  double3 tr02 = vert[0]->co - vert[2]->co;
204  double3 tr12 = vert[1]->co - vert[2]->co;
205  normal = math::cross(tr02, tr12);
206  }
207  double d = -math::dot(normal, vert[0]->co);
208  plane = new Plane(normal, d);
209  }
210 }
211 
212 Face::~Face()
213 {
214  delete plane;
215 }
216 
217 bool Face::operator==(const Face &other) const
218 {
219  if (this->size() != other.size()) {
220  return false;
221  }
222  for (FacePos i : index_range()) {
223  /* Can test pointer equality since we will have
224  * unique vert pointers for unique co_equal's. */
225  if (this->vert[i] != other.vert[i]) {
226  return false;
227  }
228  }
229  return true;
230 }
231 
232 bool Face::cyclic_equal(const Face &other) const
233 {
234  if (this->size() != other.size()) {
235  return false;
236  }
237  int flen = this->size();
238  for (FacePos start : index_range()) {
239  for (FacePos start_other : index_range()) {
240  bool ok = true;
241  for (int i = 0; ok && i < flen; ++i) {
242  FacePos p = (start + i) % flen;
243  FacePos p_other = (start_other + i) % flen;
244  if (this->vert[p] != other.vert[p_other]) {
245  ok = false;
246  }
247  }
248  if (ok) {
249  return true;
250  }
251  }
252  }
253  return false;
254 }
255 
256 std::ostream &operator<<(std::ostream &os, const Face *f)
257 {
258  os << "f" << f->id << "o" << f->orig << "[";
259  for (const Vert *v : *f) {
260  os << v;
261  if (v != f->vert[f->size() - 1]) {
262  os << " ";
263  }
264  }
265  os << "]";
266  if (f->orig != NO_INDEX) {
267  os << "o" << f->orig;
268  }
269  os << " e_orig[";
270  for (int i : f->index_range()) {
271  os << f->edge_orig[i];
272  if (f->is_intersect[i]) {
273  os << "#";
274  }
275  if (i != f->size() - 1) {
276  os << " ";
277  }
278  }
279  os << "]";
280  return os;
281 }
282 
289 // #define USE_SPINLOCK
290 
298 class IMeshArena::IMeshArenaImpl : NonCopyable, NonMovable {
299 
305  struct VSetKey {
306  Vert *vert;
307 
308  VSetKey(Vert *p) : vert(p)
309  {
310  }
311 
312  uint64_t hash() const
313  {
314  return vert->hash();
315  }
316 
317  bool operator==(const VSetKey &other) const
318  {
319  return *this->vert == *other.vert;
320  }
321  };
322 
323  Set<VSetKey> vset_;
324 
330  Vector<std::unique_ptr<Vert>> allocated_verts_;
331  Vector<std::unique_ptr<Face>> allocated_faces_;
332 
333  /* Use these to allocate ids when Verts and Faces are allocated. */
334  int next_vert_id_ = 0;
335  int next_face_id_ = 0;
336 
337  /* Need a lock when multi-threading to protect allocation of new elements. */
338 # ifdef USE_SPINLOCK
339  SpinLock lock_;
340 # else
341  ThreadMutex *mutex_;
342 # endif
343 
344  public:
345  IMeshArenaImpl()
346  {
347  if (intersect_use_threading) {
348 # ifdef USE_SPINLOCK
349  BLI_spin_init(&lock_);
350 # else
351  mutex_ = BLI_mutex_alloc();
352 # endif
353  }
354  }
355  ~IMeshArenaImpl()
356  {
357  if (intersect_use_threading) {
358 # ifdef USE_SPINLOCK
359  BLI_spin_end(&lock_);
360 # else
361  BLI_mutex_free(mutex_);
362 # endif
363  }
364  }
365 
366  void reserve(int vert_num_hint, int face_num_hint)
367  {
368  vset_.reserve(vert_num_hint);
369  allocated_verts_.reserve(vert_num_hint);
370  allocated_faces_.reserve(face_num_hint);
371  }
372 
373  int tot_allocated_verts() const
374  {
375  return allocated_verts_.size();
376  }
377 
378  int tot_allocated_faces() const
379  {
380  return allocated_faces_.size();
381  }
382 
383  const Vert *add_or_find_vert(const mpq3 &co, int orig)
384  {
385  double3 dco(co[0].get_d(), co[1].get_d(), co[2].get_d());
386  return add_or_find_vert(co, dco, orig);
387  }
388 
389  const Vert *add_or_find_vert(const double3 &co, int orig)
390  {
391  mpq3 mco(co[0], co[1], co[2]);
392  return add_or_find_vert(mco, co, orig);
393  }
394 
395  const Vert *add_or_find_vert(Vert *vert)
396  {
397  return add_or_find_vert_(vert);
398  }
399 
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 
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 
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 
434  const Vert *find_vert(const mpq3 &co)
435  {
436  Vert vtry(co, double3(co[0].get_d(), co[1].get_d(), co[2].get_d()), NO_INDEX, NO_INDEX);
437  VSetKey vskey(&vtry);
438  if (intersect_use_threading) {
439 # ifdef USE_SPINLOCK
440  BLI_spin_lock(&lock_);
441 # else
442  BLI_mutex_lock(mutex_);
443 # endif
444  }
445  const VSetKey *lookup = vset_.lookup_key_ptr(vskey);
446  if (intersect_use_threading) {
447 # ifdef USE_SPINLOCK
448  BLI_spin_unlock(&lock_);
449 # else
450  BLI_mutex_unlock(mutex_);
451 # endif
452  }
453  if (!lookup) {
454  return nullptr;
455  }
456  return lookup->vert;
457  }
458 
464  const Face *find_face(Span<const Vert *> vs)
465  {
466  Array<int> eorig(vs.size(), NO_INDEX);
467  Array<bool> is_intersect(vs.size(), false);
468  Face ftry(vs, NO_INDEX, NO_INDEX, eorig, is_intersect);
469  for (const int i : allocated_faces_.index_range()) {
470  if (ftry.cyclic_equal(*allocated_faces_[i])) {
471  return allocated_faces_[i].get();
472  }
473  }
474  return nullptr;
475  }
476 
477  private:
478  const Vert *add_or_find_vert(const mpq3 &mco, const double3 &dco, int orig)
479  {
480  Vert *vtry = new Vert(mco, dco, NO_INDEX, NO_INDEX);
481  const Vert *ans;
482  VSetKey vskey(vtry);
483  if (intersect_use_threading) {
484 # ifdef USE_SPINLOCK
485  BLI_spin_lock(&lock_);
486 # else
487  BLI_mutex_lock(mutex_);
488 # endif
489  }
490  const VSetKey *lookup = vset_.lookup_key_ptr(vskey);
491  if (!lookup) {
492  vtry->id = next_vert_id_++;
493  vtry->orig = orig;
494  vskey.vert = vtry; // new Vert(mco, dco, next_vert_id_++, orig);
495  vset_.add_new(vskey);
496  allocated_verts_.append(std::unique_ptr<Vert>(vskey.vert));
497  ans = vskey.vert;
498  }
499  else {
500  /* It was a duplicate, so return the existing one.
501  * Note that the returned Vert may have a different orig.
502  * This is the intended semantics: if the Vert already
503  * exists then we are merging verts and using the first-seen
504  * one as the canonical one. */
505  delete vtry;
506  ans = lookup->vert;
507  }
508  if (intersect_use_threading) {
509 # ifdef USE_SPINLOCK
510  BLI_spin_unlock(&lock_);
511 # else
512  BLI_mutex_unlock(mutex_);
513 # endif
514  }
515  return ans;
516  };
517 
518  const Vert *add_or_find_vert_(Vert *vtry)
519  {
520  const Vert *ans;
521  VSetKey vskey(vtry);
522  if (intersect_use_threading) {
523 # ifdef USE_SPINLOCK
524  BLI_spin_lock(&lock_);
525 # else
526  BLI_mutex_lock(mutex_);
527 # endif
528  }
529  const VSetKey *lookup = vset_.lookup_key_ptr(vskey);
530  if (!lookup) {
531  vtry->id = next_vert_id_++;
532  vskey.vert = vtry; // new Vert(mco, dco, next_vert_id_++, orig);
533  vset_.add_new(vskey);
534  allocated_verts_.append(std::unique_ptr<Vert>(vskey.vert));
535  ans = vskey.vert;
536  }
537  else {
538  /* It was a duplicate, so return the existing one.
539  * Note that the returned Vert may have a different orig.
540  * This is the intended semantics: if the Vert already
541  * exists then we are merging verts and using the first-seen
542  * one as the canonical one. */
543  delete vtry;
544  ans = lookup->vert;
545  }
546  if (intersect_use_threading) {
547 # ifdef USE_SPINLOCK
548  BLI_spin_unlock(&lock_);
549 # else
550  BLI_mutex_unlock(mutex_);
551 # endif
552  }
553  return ans;
554  };
555 };
556 
557 IMeshArena::IMeshArena()
558 {
559  pimpl_ = std::make_unique<IMeshArenaImpl>();
560 }
561 
562 IMeshArena::~IMeshArena() = default;
563 
564 void IMeshArena::reserve(int vert_num_hint, int face_num_hint)
565 {
566  pimpl_->reserve(vert_num_hint, face_num_hint);
567 }
568 
569 int IMeshArena::tot_allocated_verts() const
570 {
571  return pimpl_->tot_allocated_verts();
572 }
573 
574 int IMeshArena::tot_allocated_faces() const
575 {
576  return pimpl_->tot_allocated_faces();
577 }
578 
579 const Vert *IMeshArena::add_or_find_vert(const mpq3 &co, int orig)
580 {
581  return pimpl_->add_or_find_vert(co, orig);
582 }
583 
584 const Vert *IMeshArena::add_or_find_vert(Vert *vert)
585 {
586  return pimpl_->add_or_find_vert(vert);
587 }
588 
589 Face *IMeshArena::add_face(Span<const Vert *> verts,
590  int orig,
591  Span<int> edge_origs,
592  Span<bool> is_intersect)
593 {
594  return pimpl_->add_face(verts, orig, edge_origs, is_intersect);
595 }
596 
597 Face *IMeshArena::add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs)
598 {
599  return pimpl_->add_face(verts, orig, edge_origs);
600 }
601 
602 Face *IMeshArena::add_face(Span<const Vert *> verts, int orig)
603 {
604  return pimpl_->add_face(verts, orig);
605 }
606 
607 const Vert *IMeshArena::add_or_find_vert(const double3 &co, int orig)
608 {
609  return pimpl_->add_or_find_vert(co, orig);
610 }
611 
612 const Vert *IMeshArena::find_vert(const mpq3 &co) const
613 {
614  return pimpl_->find_vert(co);
615 }
616 
617 const Face *IMeshArena::find_face(Span<const Vert *> verts) const
618 {
619  return pimpl_->find_face(verts);
620 }
621 
622 void IMesh::set_faces(Span<Face *> faces)
623 {
624  face_ = faces;
625 }
626 
627 int IMesh::lookup_vert(const Vert *v) const
628 {
629  BLI_assert(vert_populated_);
630  return vert_to_index_.lookup_default(v, NO_INDEX);
631 }
632 
633 void IMesh::populate_vert()
634 {
635  /* This is likely an overestimate, since verts are shared between
636  * faces. It is ok if estimate is over or even under. */
637  constexpr int ESTIMATE_VERTS_PER_FACE = 4;
638  int estimate_verts_num = ESTIMATE_VERTS_PER_FACE * face_.size();
639  populate_vert(estimate_verts_num);
640 }
641 
642 void IMesh::populate_vert(int max_verts)
643 {
644  if (vert_populated_) {
645  return;
646  }
647  vert_to_index_.reserve(max_verts);
648  int next_allocate_index = 0;
649  for (const Face *f : face_) {
650  for (const Vert *v : *f) {
651  if (v->id == 1) {
652  }
653  int index = vert_to_index_.lookup_default(v, NO_INDEX);
654  if (index == NO_INDEX) {
655  BLI_assert(next_allocate_index < UINT_MAX - 2);
656  vert_to_index_.add(v, next_allocate_index++);
657  }
658  }
659  }
660  int tot_v = next_allocate_index;
661  vert_ = Array<const Vert *>(tot_v);
662  for (auto item : vert_to_index_.items()) {
663  int index = item.value;
664  BLI_assert(index < tot_v);
665  vert_[index] = item.key;
666  }
667  /* Easier debugging (at least when there are no merged input verts)
668  * if output vert order is same as input, with new verts at the end.
669  * TODO: when all debugged, set fix_order = false. */
670  const bool fix_order = true;
671  if (fix_order) {
672  blender::parallel_sort(vert_.begin(), vert_.end(), [](const Vert *a, const Vert *b) {
673  if (a->orig != NO_INDEX && b->orig != NO_INDEX) {
674  return a->orig < b->orig;
675  }
676  if (a->orig != NO_INDEX) {
677  return true;
678  }
679  if (b->orig != NO_INDEX) {
680  return false;
681  }
682  return a->id < b->id;
683  });
684  for (int i : vert_.index_range()) {
685  const Vert *v = vert_[i];
686  vert_to_index_.add_overwrite(v, i);
687  }
688  }
689  vert_populated_ = true;
690 }
691 
692 bool IMesh::erase_face_positions(int f_index, Span<bool> face_pos_erase, IMeshArena *arena)
693 {
694  const Face *cur_f = this->face(f_index);
695  int cur_len = cur_f->size();
696  int to_erase_num = 0;
697  for (int i : cur_f->index_range()) {
698  if (face_pos_erase[i]) {
699  ++to_erase_num;
700  }
701  }
702  if (to_erase_num == 0) {
703  return false;
704  }
705  int new_len = cur_len - to_erase_num;
706  if (new_len < 3) {
707  /* This erase causes removal of whole face.
708  * Because this may be called from a loop over the face array,
709  * we don't want to compress that array right here; instead will
710  * mark with null pointer and caller should call remove_null_faces().
711  * the loop is done.
712  */
713  this->face_[f_index] = NULL;
714  return true;
715  }
716  Array<const Vert *> new_vert(new_len);
717  Array<int> new_edge_orig(new_len);
718  Array<bool> new_is_intersect(new_len);
719  int new_index = 0;
720  for (int i : cur_f->index_range()) {
721  if (!face_pos_erase[i]) {
722  new_vert[new_index] = (*cur_f)[i];
723  new_edge_orig[new_index] = cur_f->edge_orig[i];
724  new_is_intersect[new_index] = cur_f->is_intersect[i];
725  ++new_index;
726  }
727  }
728  BLI_assert(new_index == new_len);
729  this->face_[f_index] = arena->add_face(new_vert, cur_f->orig, new_edge_orig, new_is_intersect);
730  return false;
731 }
732 
733 void IMesh::remove_null_faces()
734 {
735  int64_t nullcount = 0;
736  for (Face *f : this->face_) {
737  if (f == nullptr) {
738  ++nullcount;
739  }
740  }
741  if (nullcount == 0) {
742  return;
743  }
744  int64_t new_size = this->face_.size() - nullcount;
745  int64_t copy_to_index = 0;
746  int64_t copy_from_index = 0;
747  Array<Face *> new_face(new_size);
748  while (copy_from_index < face_.size()) {
749  Face *f_from = face_[copy_from_index++];
750  if (f_from) {
751  new_face[copy_to_index++] = f_from;
752  }
753  }
754  this->face_ = new_face;
755 }
756 
757 std::ostream &operator<<(std::ostream &os, const IMesh &mesh)
758 {
759  if (mesh.has_verts()) {
760  os << "Verts:\n";
761  int i = 0;
762  for (const Vert *v : mesh.vertices()) {
763  os << i << ": " << v << "\n";
764  ++i;
765  }
766  }
767  os << "\nFaces:\n";
768  int i = 0;
769  for (const Face *f : mesh.faces()) {
770  os << i << ": " << f << "\n";
771  if (f->plane != nullptr) {
772  os << " plane=" << f->plane << " eorig=[";
773  for (Face::FacePos p = 0; p < f->size(); ++p) {
774  os << f->edge_orig[p] << " ";
775  }
776  os << "]\n";
777  }
778  ++i;
779  }
780  return os;
781 }
782 
783 bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b)
784 {
785  return isect_aabb_aabb_v3(bb_a.min, bb_a.max, bb_b.min, bb_b.max);
786 }
787 
794 struct BBChunkData {
795  double max_abs_val = 0.0;
796 };
797 
798 struct BBCalcData {
799  const IMesh &im;
800  Array<BoundingBox> *face_bounding_box;
801 
802  BBCalcData(const IMesh &im, Array<BoundingBox> *fbb) : im(im), face_bounding_box(fbb){};
803 };
804 
805 static void calc_face_bb_range_func(void *__restrict userdata,
806  const int iter,
807  const TaskParallelTLS *__restrict tls)
808 {
809  BBCalcData *bbdata = static_cast<BBCalcData *>(userdata);
810  double max_abs = 0.0;
811  const Face &face = *bbdata->im.face(iter);
812  BoundingBox &bb = (*bbdata->face_bounding_box)[iter];
813  for (const Vert *v : face) {
814  bb.combine(v->co);
815  for (int i = 0; i < 3; ++i) {
816  max_abs = max_dd(max_abs, fabs(v->co[i]));
817  }
818  }
819  BBChunkData *chunk = static_cast<BBChunkData *>(tls->userdata_chunk);
820  chunk->max_abs_val = max_dd(max_abs, chunk->max_abs_val);
821 }
822 
823 struct BBPadData {
824  Array<BoundingBox> *face_bounding_box;
825  double pad;
826 
827  BBPadData(Array<BoundingBox> *fbb, double pad) : face_bounding_box(fbb), pad(pad){};
828 };
829 
830 static void pad_face_bb_range_func(void *__restrict userdata,
831  const int iter,
832  const TaskParallelTLS *__restrict UNUSED(tls))
833 {
834  BBPadData *pad_data = static_cast<BBPadData *>(userdata);
835  (*pad_data->face_bounding_box)[iter].expand(pad_data->pad);
836 }
837 
838 static void calc_face_bb_reduce(const void *__restrict UNUSED(userdata),
839  void *__restrict chunk_join,
840  void *__restrict chunk)
841 {
842  BBChunkData *bbchunk_join = static_cast<BBChunkData *>(chunk_join);
843  BBChunkData *bbchunk = static_cast<BBChunkData *>(chunk);
844  bbchunk_join->max_abs_val = max_dd(bbchunk_join->max_abs_val, bbchunk->max_abs_val);
845 }
846 
852 static Array<BoundingBox> calc_face_bounding_boxes(const IMesh &m)
853 {
854  int n = m.face_size();
855  Array<BoundingBox> ans(n);
856  TaskParallelSettings settings;
857  BBCalcData data(m, &ans);
858  BBChunkData chunk_data;
860  settings.userdata_chunk = &chunk_data;
861  settings.userdata_chunk_size = sizeof(chunk_data);
862  settings.func_reduce = calc_face_bb_reduce;
863  settings.min_iter_per_thread = 1000;
864  settings.use_threading = intersect_use_threading;
865  BLI_task_parallel_range(0, n, &data, calc_face_bb_range_func, &settings);
866  double max_abs_val = chunk_data.max_abs_val;
867  constexpr float pad_factor = 10.0f;
868  float pad = max_abs_val == 0.0f ? FLT_EPSILON : 2 * FLT_EPSILON * max_abs_val;
869  pad *= pad_factor; /* For extra safety. */
870  TaskParallelSettings pad_settings;
872  settings.min_iter_per_thread = 1000;
873  settings.use_threading = intersect_use_threading;
874  BBPadData pad_data(&ans, pad);
875  BLI_task_parallel_range(0, n, &pad_data, pad_face_bb_range_func, &pad_settings);
876  return ans;
877 }
878 
888 class CoplanarCluster {
889  Vector<int> tris_;
890  BoundingBox bb_;
891 
892  public:
893  CoplanarCluster() = default;
894  CoplanarCluster(int t, const BoundingBox &bb)
895  {
896  this->add_tri(t, bb);
897  }
898 
899  /* Assume that caller knows this will not be a duplicate. */
900  void add_tri(int t, const BoundingBox &bb)
901  {
902  tris_.append(t);
903  bb_.combine(bb);
904  }
905  int tot_tri() const
906  {
907  return tris_.size();
908  }
909  int tri(int index) const
910  {
911  return tris_[index];
912  }
913  const int *begin() const
914  {
915  return tris_.begin();
916  }
917  const int *end() const
918  {
919  return tris_.end();
920  }
921 
922  const BoundingBox &bounding_box() const
923  {
924  return bb_;
925  }
926 };
927 
934 class CoplanarClusterInfo {
935  Vector<CoplanarCluster> clusters_;
936  Array<int> tri_cluster_;
937 
938  public:
939  CoplanarClusterInfo() = default;
940  explicit CoplanarClusterInfo(int numtri) : tri_cluster_(Array<int>(numtri))
941  {
942  tri_cluster_.fill(-1);
943  }
944 
945  int tri_cluster(int t) const
946  {
947  BLI_assert(t < tri_cluster_.size());
948  return tri_cluster_[t];
949  }
950 
951  int add_cluster(CoplanarCluster cl)
952  {
953  int c_index = clusters_.append_and_get_index(cl);
954  for (int t : cl) {
955  BLI_assert(t < tri_cluster_.size());
956  tri_cluster_[t] = c_index;
957  }
958  return c_index;
959  }
960 
961  int tot_cluster() const
962  {
963  return clusters_.size();
964  }
965 
966  const CoplanarCluster *begin()
967  {
968  return clusters_.begin();
969  }
970 
971  const CoplanarCluster *end()
972  {
973  return clusters_.end();
974  }
975 
976  IndexRange index_range() const
977  {
978  return clusters_.index_range();
979  }
980 
981  const CoplanarCluster &cluster(int index) const
982  {
983  BLI_assert(index < clusters_.size());
984  return clusters_[index];
985  }
986 };
987 
988 static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl);
989 
990 static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo);
991 
992 enum ITT_value_kind { INONE, IPOINT, ISEGMENT, ICOPLANAR };
993 
994 struct ITT_value {
995  mpq3 p1; /* Only relevant for IPOINT and ISEGMENT kind. */
996  mpq3 p2; /* Only relevant for ISEGMENT kind. */
997  int t_source = -1; /* Index of the source triangle that intersected the target one. */
998  enum ITT_value_kind kind = INONE;
999 
1000  ITT_value() = default;
1001  explicit ITT_value(ITT_value_kind k) : kind(k)
1002  {
1003  }
1004  ITT_value(ITT_value_kind k, int tsrc) : t_source(tsrc), kind(k)
1005  {
1006  }
1007  ITT_value(ITT_value_kind k, const mpq3 &p1) : p1(p1), kind(k)
1008  {
1009  }
1010  ITT_value(ITT_value_kind k, const mpq3 &p1, const mpq3 &p2) : p1(p1), p2(p2), kind(k)
1011  {
1012  }
1013 };
1014 
1015 static std::ostream &operator<<(std::ostream &os, const ITT_value &itt);
1016 
1022 static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
1023 {
1024  mpq2 p2d;
1025  switch (proj_axis) {
1026  case (0): {
1027  p2d[0] = p3d[1];
1028  p2d[1] = p3d[2];
1029  break;
1030  }
1031  case (1): {
1032  p2d[0] = p3d[0];
1033  p2d[1] = p3d[2];
1034  break;
1035  }
1036  case (2): {
1037  p2d[0] = p3d[0];
1038  p2d[1] = p3d[1];
1039  break;
1040  }
1041  default:
1042  BLI_assert(false);
1043  }
1044  return p2d;
1045 }
1046 
1077 static double supremum_dot_cross(const double3 &a, const double3 &b)
1078 {
1079  double3 abs_a = math::abs(a);
1080  double3 abs_b = math::abs(b);
1081  double3 c;
1082  /* This is dot(cross(a, b), cross(a,b)) but using absolute values for a and b
1083  * and always using + when operation is + or -. */
1084  c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1];
1085  c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2];
1086  c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0];
1087  return math::dot(c, c);
1088 }
1089 
1090 /* The index of dot when inputs are plane_coords with index 1 is much higher.
1091  * Plane coords have index 6.
1092  */
1093 constexpr int index_dot_plane_coords = 15;
1094 
1105 constexpr int index_dot_cross = 11;
1106 
1115 constexpr int index_plane_side = 3 + 2 * index_dot_plane_coords;
1116 
1117 static int filter_plane_side(const double3 &p,
1118  const double3 &plane_p,
1119  const double3 &plane_no,
1120  const double3 &abs_p,
1121  const double3 &abs_plane_p,
1122  const double3 &abs_plane_no)
1123 {
1124  double d = math::dot(p - plane_p, plane_no);
1125  if (d == 0.0) {
1126  return 0;
1127  }
1128  double supremum = math::dot(abs_p + abs_plane_p, abs_plane_no);
1129  double err_bound = supremum * index_plane_side * DBL_EPSILON;
1130  if (fabs(d) > err_bound) {
1131  return d > 0 ? 1 : -1;
1132  }
1133  return 0;
1134 }
1135 
1136 /*
1137  * #intersect_tri_tri and helper functions.
1138  * This code uses the algorithm of Guigue and Devillers, as described
1139  * in "Faster Triangle-Triangle Intersection Tests".
1140  * Adapted from code by Eric Haines:
1141  * https://github.com/erich666/jgt-code/tree/master/Volume_08/Number_1/Guigue2003
1142  */
1143 
1152 static inline mpq3 tti_interp(
1153  const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &n, mpq3 &ab, mpq3 &ac, mpq3 &dotbuf)
1154 {
1155  ab = a;
1156  ab -= b;
1157  ac = a;
1158  ac -= c;
1159  mpq_class den = math::dot_with_buffer(ab, n, dotbuf);
1160  BLI_assert(den != 0);
1161  mpq_class alpha = math::dot_with_buffer(ac, n, dotbuf) / den;
1162  return a - alpha * ab;
1163 }
1164 
1172 static inline int tti_above(const mpq3 &a,
1173  const mpq3 &b,
1174  const mpq3 &c,
1175  const mpq3 &ad,
1176  mpq3 &ba,
1177  mpq3 &ca,
1178  mpq3 &n,
1179  mpq3 &dotbuf)
1180 {
1181  ba = b;
1182  ba -= a;
1183  ca = c;
1184  ca -= a;
1185 
1186  n.x = ba.y * ca.z - ba.z * ca.y;
1187  n.y = ba.z * ca.x - ba.x * ca.z;
1188  n.z = ba.x * ca.y - ba.y * ca.x;
1189 
1190  return sgn(math::dot_with_buffer(ad, n, dotbuf));
1191 }
1192 
1206 static ITT_value itt_canon2(const mpq3 &p1,
1207  const mpq3 &q1,
1208  const mpq3 &r1,
1209  const mpq3 &p2,
1210  const mpq3 &q2,
1211  const mpq3 &r2,
1212  const mpq3 &n1,
1213  const mpq3 &n2)
1214 {
1215  constexpr int dbg_level = 0;
1216  if (dbg_level > 0) {
1217  std::cout << "\ntri_tri_intersect_canon:\n";
1218  std::cout << "p1=" << p1 << " q1=" << q1 << " r1=" << r1 << "\n";
1219  std::cout << "p2=" << p2 << " q2=" << q2 << " r2=" << r2 << "\n";
1220  std::cout << "n1=" << n1 << " n2=" << n2 << "\n";
1221  std::cout << "approximate values:\n";
1222  std::cout << "p1=(" << p1[0].get_d() << "," << p1[1].get_d() << "," << p1[2].get_d() << ")\n";
1223  std::cout << "q1=(" << q1[0].get_d() << "," << q1[1].get_d() << "," << q1[2].get_d() << ")\n";
1224  std::cout << "r1=(" << r1[0].get_d() << "," << r1[1].get_d() << "," << r1[2].get_d() << ")\n";
1225  std::cout << "p2=(" << p2[0].get_d() << "," << p2[1].get_d() << "," << p2[2].get_d() << ")\n";
1226  std::cout << "q2=(" << q2[0].get_d() << "," << q2[1].get_d() << "," << q2[2].get_d() << ")\n";
1227  std::cout << "r2=(" << r2[0].get_d() << "," << r2[1].get_d() << "," << r2[2].get_d() << ")\n";
1228  std::cout << "n1=(" << n1[0].get_d() << "," << n1[1].get_d() << "," << n1[2].get_d() << ")\n";
1229  std::cout << "n2=(" << n2[0].get_d() << "," << n2[1].get_d() << "," << n2[2].get_d() << ")\n";
1230  }
1231  mpq3 p1p2 = p2 - p1;
1232  mpq3 intersect_1;
1233  mpq3 intersect_2;
1234  mpq3 buf[4];
1235  bool no_overlap = false;
1236  /* Top test in classification tree. */
1237  if (tti_above(p1, q1, r2, p1p2, buf[0], buf[1], buf[2], buf[3]) > 0) {
1238  /* Middle right test in classification tree. */
1239  if (tti_above(p1, r1, r2, p1p2, buf[0], buf[1], buf[2], buf[3]) <= 0) {
1240  /* Bottom right test in classification tree. */
1241  if (tti_above(p1, r1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) > 0) {
1242  /* Overlap is [k [i l] j]. */
1243  if (dbg_level > 0) {
1244  std::cout << "overlap [k [i l] j]\n";
1245  }
1246  /* i is intersect with p1r1. l is intersect with p2r2. */
1247  intersect_1 = tti_interp(p1, r1, p2, n2, buf[0], buf[1], buf[2]);
1248  intersect_2 = tti_interp(p2, r2, p1, n1, buf[0], buf[1], buf[2]);
1249  }
1250  else {
1251  /* Overlap is [i [k l] j]. */
1252  if (dbg_level > 0) {
1253  std::cout << "overlap [i [k l] j]\n";
1254  }
1255  /* k is intersect with p2q2. l is intersect is p2r2. */
1256  intersect_1 = tti_interp(p2, q2, p1, n1, buf[0], buf[1], buf[2]);
1257  intersect_2 = tti_interp(p2, r2, p1, n1, buf[0], buf[1], buf[2]);
1258  }
1259  }
1260  else {
1261  /* No overlap: [k l] [i j]. */
1262  if (dbg_level > 0) {
1263  std::cout << "no overlap: [k l] [i j]\n";
1264  }
1265  no_overlap = true;
1266  }
1267  }
1268  else {
1269  /* Middle left test in classification tree. */
1270  if (tti_above(p1, q1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) < 0) {
1271  /* No overlap: [i j] [k l]. */
1272  if (dbg_level > 0) {
1273  std::cout << "no overlap: [i j] [k l]\n";
1274  }
1275  no_overlap = true;
1276  }
1277  else {
1278  /* Bottom left test in classification tree. */
1279  if (tti_above(p1, r1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) >= 0) {
1280  /* Overlap is [k [i j] l]. */
1281  if (dbg_level > 0) {
1282  std::cout << "overlap [k [i j] l]\n";
1283  }
1284  /* i is intersect with p1r1. j is intersect with p1q1. */
1285  intersect_1 = tti_interp(p1, r1, p2, n2, buf[0], buf[1], buf[2]);
1286  intersect_2 = tti_interp(p1, q1, p2, n2, buf[0], buf[1], buf[2]);
1287  }
1288  else {
1289  /* Overlap is [i [k j] l]. */
1290  if (dbg_level > 0) {
1291  std::cout << "overlap [i [k j] l]\n";
1292  }
1293  /* k is intersect with p2q2. j is intersect with p1q1. */
1294  intersect_1 = tti_interp(p2, q2, p1, n1, buf[0], buf[1], buf[2]);
1295  intersect_2 = tti_interp(p1, q1, p2, n2, buf[0], buf[1], buf[2]);
1296  }
1297  }
1298  }
1299  if (no_overlap) {
1300  return ITT_value(INONE);
1301  }
1302  if (intersect_1 == intersect_2) {
1303  if (dbg_level > 0) {
1304  std::cout << "single intersect: " << intersect_1 << "\n";
1305  }
1306  return ITT_value(IPOINT, intersect_1);
1307  }
1308  if (dbg_level > 0) {
1309  std::cout << "intersect segment: " << intersect_1 << ", " << intersect_2 << "\n";
1310  }
1311  return ITT_value(ISEGMENT, intersect_1, intersect_2);
1312 }
1313 
1314 /* Helper function for intersect_tri_tri. Arguments have been canonicalized for triangle 1. */
1315 
1316 static ITT_value itt_canon1(const mpq3 &p1,
1317  const mpq3 &q1,
1318  const mpq3 &r1,
1319  const mpq3 &p2,
1320  const mpq3 &q2,
1321  const mpq3 &r2,
1322  const mpq3 &n1,
1323  const mpq3 &n2,
1324  int sp2,
1325  int sq2,
1326  int sr2)
1327 {
1328  constexpr int dbg_level = 0;
1329  if (sp2 > 0) {
1330  if (sq2 > 0) {
1331  return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2);
1332  }
1333  if (sr2 > 0) {
1334  return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2);
1335  }
1336  return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2);
1337  }
1338  if (sp2 < 0) {
1339  if (sq2 < 0) {
1340  return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2);
1341  }
1342  if (sr2 < 0) {
1343  return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2);
1344  }
1345  return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2);
1346  }
1347  if (sq2 < 0) {
1348  if (sr2 >= 0) {
1349  return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2);
1350  }
1351  return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2);
1352  }
1353  if (sq2 > 0) {
1354  if (sr2 > 0) {
1355  return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2);
1356  }
1357  return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2);
1358  }
1359  if (sr2 > 0) {
1360  return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2);
1361  }
1362  if (sr2 < 0) {
1363  return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2);
1364  }
1365  if (dbg_level > 0) {
1366  std::cout << "triangles are co-planar\n";
1367  }
1368  return ITT_value(ICOPLANAR);
1369 }
1370 
1371 static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
1372 {
1373  constexpr int dbg_level = 0;
1374 # ifdef PERFDEBUG
1375  incperfcount(1); /* Intersect_tri_tri calls. */
1376 # endif
1377  const Face &tri1 = *tm.face(t1);
1378  const Face &tri2 = *tm.face(t2);
1379  BLI_assert(tri1.plane_populated() && tri2.plane_populated());
1380  const Vert *vp1 = tri1[0];
1381  const Vert *vq1 = tri1[1];
1382  const Vert *vr1 = tri1[2];
1383  const Vert *vp2 = tri2[0];
1384  const Vert *vq2 = tri2[1];
1385  const Vert *vr2 = tri2[2];
1386  if (dbg_level > 0) {
1387  std::cout << "\nINTERSECT_TRI_TRI t1=" << t1 << ", t2=" << t2 << "\n";
1388  std::cout << " p1 = " << vp1 << "\n";
1389  std::cout << " q1 = " << vq1 << "\n";
1390  std::cout << " r1 = " << vr1 << "\n";
1391  std::cout << " p2 = " << vp2 << "\n";
1392  std::cout << " q2 = " << vq2 << "\n";
1393  std::cout << " r2 = " << vr2 << "\n";
1394  }
1395 
1396  /* Get signs of t1's vertices' distances to plane of t2 and vice versa. */
1397 
1398  /* Try first getting signs with double arithmetic, with error bounds.
1399  * If the signs calculated in this section are not 0, they are the same
1400  * as what they would be using exact arithmetic. */
1401  const double3 &d_p1 = vp1->co;
1402  const double3 &d_q1 = vq1->co;
1403  const double3 &d_r1 = vr1->co;
1404  const double3 &d_p2 = vp2->co;
1405  const double3 &d_q2 = vq2->co;
1406  const double3 &d_r2 = vr2->co;
1407  const double3 &d_n2 = tri2.plane->norm;
1408 
1409  const double3 &abs_d_p1 = math::abs(d_p1);
1410  const double3 &abs_d_q1 = math::abs(d_q1);
1411  const double3 &abs_d_r1 = math::abs(d_r1);
1412  const double3 &abs_d_r2 = math::abs(d_r2);
1413  const double3 &abs_d_n2 = math::abs(d_n2);
1414 
1415  int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2);
1416  int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2);
1417  int sr1 = filter_plane_side(d_r1, d_r2, d_n2, abs_d_r1, abs_d_r2, abs_d_n2);
1418  if ((sp1 > 0 && sq1 > 0 && sr1 > 0) || (sp1 < 0 && sq1 < 0 && sr1 < 0)) {
1419 # ifdef PERFDEBUG
1420  incperfcount(2); /* Triangle-triangle intersects decided by filter plane tests. */
1421 # endif
1422  if (dbg_level > 0) {
1423  std::cout << "no intersection, all t1's verts above or below t2\n";
1424  }
1425  return ITT_value(INONE);
1426  }
1427 
1428  const double3 &d_n1 = tri1.plane->norm;
1429  const double3 &abs_d_p2 = math::abs(d_p2);
1430  const double3 &abs_d_q2 = math::abs(d_q2);
1431  const double3 &abs_d_n1 = math::abs(d_n1);
1432 
1433  int sp2 = filter_plane_side(d_p2, d_r1, d_n1, abs_d_p2, abs_d_r1, abs_d_n1);
1434  int sq2 = filter_plane_side(d_q2, d_r1, d_n1, abs_d_q2, abs_d_r1, abs_d_n1);
1435  int sr2 = filter_plane_side(d_r2, d_r1, d_n1, abs_d_r2, abs_d_r1, abs_d_n1);
1436  if ((sp2 > 0 && sq2 > 0 && sr2 > 0) || (sp2 < 0 && sq2 < 0 && sr2 < 0)) {
1437 # ifdef PERFDEBUG
1438  incperfcount(2); /* Triangle-triangle intersects decided by filter plane tests. */
1439 # endif
1440  if (dbg_level > 0) {
1441  std::cout << "no intersection, all t2's verts above or below t1\n";
1442  }
1443  return ITT_value(INONE);
1444  }
1445 
1446  mpq3 buf[2];
1447  const mpq3 &p1 = vp1->co_exact;
1448  const mpq3 &q1 = vq1->co_exact;
1449  const mpq3 &r1 = vr1->co_exact;
1450  const mpq3 &p2 = vp2->co_exact;
1451  const mpq3 &q2 = vq2->co_exact;
1452  const mpq3 &r2 = vr2->co_exact;
1453 
1454  const mpq3 &n2 = tri2.plane->norm_exact;
1455  if (sp1 == 0) {
1456  buf[0] = p1;
1457  buf[0] -= r2;
1458  sp1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
1459  }
1460  if (sq1 == 0) {
1461  buf[0] = q1;
1462  buf[0] -= r2;
1463  sq1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
1464  }
1465  if (sr1 == 0) {
1466  buf[0] = r1;
1467  buf[0] -= r2;
1468  sr1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
1469  }
1470 
1471  if (dbg_level > 1) {
1472  std::cout << " sp1=" << sp1 << " sq1=" << sq1 << " sr1=" << sr1 << "\n";
1473  }
1474 
1475  if ((sp1 * sq1 > 0) && (sp1 * sr1 > 0)) {
1476  if (dbg_level > 0) {
1477  std::cout << "no intersection, all t1's verts above or below t2 (exact)\n";
1478  }
1479 # ifdef PERFDEBUG
1480  incperfcount(3); /* Triangle-triangle intersects decided by exact plane tests. */
1481 # endif
1482  return ITT_value(INONE);
1483  }
1484 
1485  /* Repeat for signs of t2's vertices with respect to plane of t1. */
1486  const mpq3 &n1 = tri1.plane->norm_exact;
1487  if (sp2 == 0) {
1488  buf[0] = p2;
1489  buf[0] -= r1;
1490  sp2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
1491  }
1492  if (sq2 == 0) {
1493  buf[0] = q2;
1494  buf[0] -= r1;
1495  sq2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
1496  }
1497  if (sr2 == 0) {
1498  buf[0] = r2;
1499  buf[0] -= r1;
1500  sr2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
1501  }
1502 
1503  if (dbg_level > 1) {
1504  std::cout << " sp2=" << sp2 << " sq2=" << sq2 << " sr2=" << sr2 << "\n";
1505  }
1506 
1507  if ((sp2 * sq2 > 0) && (sp2 * sr2 > 0)) {
1508  if (dbg_level > 0) {
1509  std::cout << "no intersection, all t2's verts above or below t1 (exact)\n";
1510  }
1511 # ifdef PERFDEBUG
1512  incperfcount(3); /* Triangle-triangle intersects decided by exact plane tests. */
1513 # endif
1514  return ITT_value(INONE);
1515  }
1516 
1517  /* Do rest of the work with vertices in a canonical order, where p1 is on
1518  * positive side of plane and q1, r1 are not, or p1 is on the plane and
1519  * q1 and r1 are off the plane on the same side. */
1520  ITT_value ans;
1521  if (sp1 > 0) {
1522  if (sq1 > 0) {
1523  ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1524  }
1525  else if (sr1 > 0) {
1526  ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1527  }
1528  else {
1529  ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1530  }
1531  }
1532  else if (sp1 < 0) {
1533  if (sq1 < 0) {
1534  ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1535  }
1536  else if (sr1 < 0) {
1537  ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1538  }
1539  else {
1540  ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1541  }
1542  }
1543  else {
1544  if (sq1 < 0) {
1545  if (sr1 >= 0) {
1546  ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1547  }
1548  else {
1549  ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1550  }
1551  }
1552  else if (sq1 > 0) {
1553  if (sr1 > 0) {
1554  ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1555  }
1556  else {
1557  ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1558  }
1559  }
1560  else {
1561  if (sr1 > 0) {
1562  ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2);
1563  }
1564  else if (sr1 < 0) {
1565  ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2);
1566  }
1567  else {
1568  if (dbg_level > 0) {
1569  std::cout << "triangles are co-planar\n";
1570  }
1571  ans = ITT_value(ICOPLANAR);
1572  }
1573  }
1574  }
1575  if (ans.kind == ICOPLANAR) {
1576  ans.t_source = t2;
1577  }
1578 
1579 # ifdef PERFDEBUG
1580  if (ans.kind != INONE) {
1581  incperfcount(4);
1582  }
1583 # endif
1584  return ans;
1585 }
1586 
1587 struct CDT_data {
1588  const Plane *t_plane;
1589  Vector<mpq2> vert;
1590  Vector<std::pair<int, int>> edge;
1591  Vector<Vector<int>> face;
1593  Vector<int> input_face;
1595  Vector<bool> is_reversed;
1597  CDT_result<mpq_class> cdt_out;
1600  Map<std::pair<int, int>, int> verts_to_edge;
1601  int proj_axis;
1602 };
1603 
1607 static int prepare_need_vert(CDT_data &cd, const mpq3 &p3d)
1608 {
1609  mpq2 p2d = project_3d_to_2d(p3d, cd.proj_axis);
1610  int v = cd.vert.append_and_get_index(p2d);
1611  return v;
1612 }
1613 
1621 static mpq3 unproject_cdt_vert(const CDT_data &cd, const mpq2 &p2d)
1622 {
1623  mpq3 p3d;
1624  BLI_assert(cd.t_plane->exact_populated());
1625  BLI_assert(cd.t_plane->norm_exact[cd.proj_axis] != 0);
1626  const mpq3 &n = cd.t_plane->norm_exact;
1627  const mpq_class &d = cd.t_plane->d_exact;
1628  switch (cd.proj_axis) {
1629  case (0): {
1630  mpq_class num = n[1] * p2d[0] + n[2] * p2d[1] + d;
1631  num = -num;
1632  p3d[0] = num / n[0];
1633  p3d[1] = p2d[0];
1634  p3d[2] = p2d[1];
1635  break;
1636  }
1637  case (1): {
1638  p3d[0] = p2d[0];
1639  mpq_class num = n[0] * p2d[0] + n[2] * p2d[1] + d;
1640  num = -num;
1641  p3d[1] = num / n[1];
1642  p3d[2] = p2d[1];
1643  break;
1644  }
1645  case (2): {
1646  p3d[0] = p2d[0];
1647  p3d[1] = p2d[1];
1648  mpq_class num = n[0] * p2d[0] + n[1] * p2d[1] + d;
1649  num = -num;
1650  p3d[2] = num / n[2];
1651  break;
1652  }
1653  default:
1654  BLI_assert(false);
1655  }
1656  return p3d;
1657 }
1658 
1659 static void prepare_need_edge(CDT_data &cd, const mpq3 &p1, const mpq3 &p2)
1660 {
1661  int v1 = prepare_need_vert(cd, p1);
1662  int v2 = prepare_need_vert(cd, p2);
1663  cd.edge.append(std::pair<int, int>(v1, v2));
1664 }
1665 
1666 static void prepare_need_tri(CDT_data &cd, const IMesh &tm, int t)
1667 {
1668  const Face &tri = *tm.face(t);
1669  int v0 = prepare_need_vert(cd, tri[0]->co_exact);
1670  int v1 = prepare_need_vert(cd, tri[1]->co_exact);
1671  int v2 = prepare_need_vert(cd, tri[2]->co_exact);
1672  bool rev;
1673  /* How to get CCW orientation of projected triangle? Note that when look down y axis
1674  * as opposed to x or z, the orientation of the other two axes is not right-and-up. */
1675  BLI_assert(cd.t_plane->exact_populated());
1676  if (tri.plane->norm_exact[cd.proj_axis] >= 0) {
1677  rev = cd.proj_axis == 1;
1678  }
1679  else {
1680  rev = cd.proj_axis != 1;
1681  }
1682  int cd_t = cd.face.append_and_get_index(Vector<int>());
1683  cd.face[cd_t].append(v0);
1684  if (rev) {
1685  cd.face[cd_t].append(v2);
1686  cd.face[cd_t].append(v1);
1687  }
1688  else {
1689  cd.face[cd_t].append(v1);
1690  cd.face[cd_t].append(v2);
1691  }
1692  cd.input_face.append(t);
1693  cd.is_reversed.append(rev);
1694 }
1695 
1696 static CDT_data prepare_cdt_input(const IMesh &tm, int t, const Vector<ITT_value> itts)
1697 {
1698  CDT_data ans;
1699  BLI_assert(tm.face(t)->plane_populated());
1700  ans.t_plane = tm.face(t)->plane;
1701  BLI_assert(ans.t_plane->exact_populated());
1702  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
1703  prepare_need_tri(ans, tm, t);
1704  for (const ITT_value &itt : itts) {
1705  switch (itt.kind) {
1706  case INONE:
1707  break;
1708  case IPOINT: {
1709  prepare_need_vert(ans, itt.p1);
1710  break;
1711  }
1712  case ISEGMENT: {
1713  prepare_need_edge(ans, itt.p1, itt.p2);
1714  break;
1715  }
1716  case ICOPLANAR: {
1717  prepare_need_tri(ans, tm, itt.t_source);
1718  break;
1719  }
1720  }
1721  }
1722  return ans;
1723 }
1724 
1725 static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm,
1726  const CoplanarClusterInfo &clinfo,
1727  int c,
1728  const Vector<ITT_value> itts)
1729 {
1730  CDT_data ans;
1731  BLI_assert(c < clinfo.tot_cluster());
1732  const CoplanarCluster &cl = clinfo.cluster(c);
1733  BLI_assert(cl.tot_tri() > 0);
1734  int t0 = cl.tri(0);
1735  BLI_assert(tm.face(t0)->plane_populated());
1736  ans.t_plane = tm.face(t0)->plane;
1737  BLI_assert(ans.t_plane->exact_populated());
1738  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
1739  for (const int t : cl) {
1740  prepare_need_tri(ans, tm, t);
1741  }
1742  for (const ITT_value &itt : itts) {
1743  switch (itt.kind) {
1744  case IPOINT: {
1745  prepare_need_vert(ans, itt.p1);
1746  break;
1747  }
1748  case ISEGMENT: {
1749  prepare_need_edge(ans, itt.p1, itt.p2);
1750  break;
1751  }
1752  default:
1753  break;
1754  }
1755  }
1756  return ans;
1757 }
1758 
1759 /* Return a copy of the argument with the integers ordered in ascending order. */
1760 static inline std::pair<int, int> sorted_int_pair(std::pair<int, int> pair)
1761 {
1762  if (pair.first <= pair.second) {
1763  return pair;
1764  }
1765  return std::pair<int, int>(pair.second, pair.first);
1766 }
1767 
1772 static void populate_cdt_edge_map(Map<std::pair<int, int>, int> &verts_to_edge,
1773  const CDT_result<mpq_class> &cdt_out)
1774 {
1775  verts_to_edge.reserve(cdt_out.edge.size());
1776  for (int e : cdt_out.edge.index_range()) {
1777  std::pair<int, int> vpair = sorted_int_pair(cdt_out.edge[e]);
1778  /* There should be only one edge for each vertex pair. */
1779  verts_to_edge.add(vpair, e);
1780  }
1781 }
1782 
1786 static void do_cdt(CDT_data &cd)
1787 {
1788  constexpr int dbg_level = 0;
1789  CDT_input<mpq_class> cdt_in;
1790  cdt_in.vert = Span<mpq2>(cd.vert);
1791  cdt_in.edge = Span<std::pair<int, int>>(cd.edge);
1792  cdt_in.face = Span<Vector<int>>(cd.face);
1793  if (dbg_level > 0) {
1794  std::cout << "CDT input\nVerts:\n";
1795  for (int i : cdt_in.vert.index_range()) {
1796  std::cout << "v" << i << ": " << cdt_in.vert[i] << "=(" << cdt_in.vert[i][0].get_d() << ","
1797  << cdt_in.vert[i][1].get_d() << ")\n";
1798  }
1799  std::cout << "Edges:\n";
1800  for (int i : cdt_in.edge.index_range()) {
1801  std::cout << "e" << i << ": (" << cdt_in.edge[i].first << ", " << cdt_in.edge[i].second
1802  << ")\n";
1803  }
1804  std::cout << "Tris\n";
1805  for (int f : cdt_in.face.index_range()) {
1806  std::cout << "f" << f << ": ";
1807  for (int j : cdt_in.face[f].index_range()) {
1808  std::cout << cdt_in.face[f][j] << " ";
1809  }
1810  std::cout << "\n";
1811  }
1812  }
1813  cdt_in.epsilon = 0; /* TODO: needs attention for non-exact T. */
1815  constexpr int make_edge_map_threshold = 15;
1816  if (cd.cdt_out.edge.size() >= make_edge_map_threshold) {
1817  populate_cdt_edge_map(cd.verts_to_edge, cd.cdt_out);
1818  }
1819  if (dbg_level > 0) {
1820  std::cout << "\nCDT result\nVerts:\n";
1821  for (int i : cd.cdt_out.vert.index_range()) {
1822  std::cout << "v" << i << ": " << cd.cdt_out.vert[i] << "=(" << cd.cdt_out.vert[i][0].get_d()
1823  << "," << cd.cdt_out.vert[i][1].get_d() << "\n";
1824  }
1825  std::cout << "Tris\n";
1826  for (int f : cd.cdt_out.face.index_range()) {
1827  std::cout << "f" << f << ": ";
1828  for (int j : cd.cdt_out.face[f].index_range()) {
1829  std::cout << cd.cdt_out.face[f][j] << " ";
1830  }
1831  std::cout << "orig: ";
1832  for (int j : cd.cdt_out.face_orig[f].index_range()) {
1833  std::cout << cd.cdt_out.face_orig[f][j] << " ";
1834  }
1835  std::cout << "\n";
1836  }
1837  std::cout << "Edges\n";
1838  for (int e : cd.cdt_out.edge.index_range()) {
1839  std::cout << "e" << e << ": (" << cd.cdt_out.edge[e].first << ", "
1840  << cd.cdt_out.edge[e].second << ") ";
1841  std::cout << "orig: ";
1842  for (int j : cd.cdt_out.edge_orig[e].index_range()) {
1843  std::cout << cd.cdt_out.edge_orig[e][j] << " ";
1844  }
1845  std::cout << "\n";
1846  }
1847  }
1848 }
1849 
1850 /* Find an original edge index that goes with the edge that the CDT output edge
1851  * that goes between verts i0 and i1 (in the CDT output vert indexing scheme).
1852  * There may be more than one: if so, prefer one that was originally a face edge.
1853  * The input to CDT for a triangle with some intersecting segments from other triangles
1854  * will have both edges and a face constraint for the main triangle (this is redundant
1855  * but allows us to discover which face edge goes with which output edges).
1856  * If there is any face edge, return one of those as the original.
1857  * If there is no face edge but there is another edge in the input problem, then that
1858  * edge must have come from intersection with another triangle, so set *r_is_intersect
1859  * to true in that case. */
1860 static int get_cdt_edge_orig(
1861  int i0, int i1, const CDT_data &cd, const IMesh &in_tm, bool *r_is_intersect)
1862 {
1863  int foff = cd.cdt_out.face_edge_offset;
1864  *r_is_intersect = false;
1865  int e = NO_INDEX;
1866  if (cd.verts_to_edge.size() > 0) {
1867  /* Use the populated map to find the edge, if any, between vertices i0 and i1. */
1868  std::pair<int, int> vpair = sorted_int_pair(std::pair<int, int>(i0, i1));
1869  e = cd.verts_to_edge.lookup_default(vpair, NO_INDEX);
1870  }
1871  else {
1872  for (int ee : cd.cdt_out.edge.index_range()) {
1873  std::pair<int, int> edge = cd.cdt_out.edge[ee];
1874  if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) {
1875  e = ee;
1876  break;
1877  }
1878  }
1879  }
1880  if (e == NO_INDEX) {
1881  return NO_INDEX;
1882  }
1883 
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 one part of the nary input,
1886  * then want to set *r_is_intersect to true. */
1887  int face_eorig = NO_INDEX;
1888  bool have_non_face_eorig = false;
1889  for (int orig_index : cd.cdt_out.edge_orig[e]) {
1890  /* orig_index encodes the triangle and pos within the triangle of the input edge. */
1891  if (orig_index >= foff) {
1892  if (face_eorig == NO_INDEX) {
1893  int in_face_index = (orig_index / foff) - 1;
1894  int pos = orig_index % foff;
1895  /* We need to retrieve the edge orig field from the Face used to populate the
1896  * in_face_index'th face of the CDT, at the pos'th position of the face. */
1897  int in_tm_face_index = cd.input_face[in_face_index];
1898  BLI_assert(in_tm_face_index < in_tm.face_size());
1899  const Face *facep = in_tm.face(in_tm_face_index);
1900  BLI_assert(pos < facep->size());
1901  bool is_rev = cd.is_reversed[in_face_index];
1902  int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
1903  if (eorig != NO_INDEX) {
1904  face_eorig = eorig;
1905  }
1906  }
1907  }
1908  else {
1909  if (!have_non_face_eorig) {
1910  have_non_face_eorig = true;
1911  }
1912  if (face_eorig != NO_INDEX && have_non_face_eorig) {
1913  /* Only need at most one orig for each type. */
1914  break;
1915  }
1916  }
1917  }
1918  if (face_eorig != NO_INDEX) {
1919  return face_eorig;
1920  }
1921  if (have_non_face_eorig) {
1922  /* This must have been an input to the CDT problem that was an intersection edge. */
1923  /* TODO: maybe there is an orig index:
1924  * This happens if an input edge was formed by an input face having
1925  * an edge that is co-planar with the cluster, while the face as a whole is not. */
1926  *r_is_intersect = true;
1927  return NO_INDEX;
1928  }
1929  return NO_INDEX;
1930 }
1931 
1937 static Face *cdt_tri_as_imesh_face(
1938  int cdt_out_t, int cdt_in_t, const CDT_data &cd, const IMesh &tm, IMeshArena *arena)
1939 {
1940  const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
1941  int t_orig = tm.face(cd.input_face[cdt_in_t])->orig;
1942  BLI_assert(cdt_out.face[cdt_out_t].size() == 3);
1943  int i0 = cdt_out.face[cdt_out_t][0];
1944  int i1 = cdt_out.face[cdt_out_t][1];
1945  int i2 = cdt_out.face[cdt_out_t][2];
1946  mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]);
1947  mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]);
1948  mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]);
1949  /* No need to provide an original index: if coord matches
1950  * an original one, then it will already be in the arena
1951  * with the correct orig field. */
1952  const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX);
1953  const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX);
1954  const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX);
1955  Face *facep;
1956  bool is_isect0;
1957  bool is_isect1;
1958  bool is_isect2;
1959  if (cd.is_reversed[cdt_in_t]) {
1960  int oe0 = get_cdt_edge_orig(i0, i2, cd, tm, &is_isect0);
1961  int oe1 = get_cdt_edge_orig(i2, i1, cd, tm, &is_isect1);
1962  int oe2 = get_cdt_edge_orig(i1, i0, cd, tm, &is_isect2);
1963  facep = arena->add_face(
1964  {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
1965  }
1966  else {
1967  int oe0 = get_cdt_edge_orig(i0, i1, cd, tm, &is_isect0);
1968  int oe1 = get_cdt_edge_orig(i1, i2, cd, tm, &is_isect1);
1969  int oe2 = get_cdt_edge_orig(i2, i0, cd, tm, &is_isect2);
1970  facep = arena->add_face(
1971  {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
1972  }
1973  facep->populate_plane(false);
1974  return facep;
1975 }
1976 
1977 /* Like BLI_math's is_quad_flip_v3_first_third_fast_with_normal, with const double3's. */
1978 static bool is_quad_flip_first_third(const double3 &v1,
1979  const double3 &v2,
1980  const double3 &v3,
1981  const double3 &v4,
1982  const double3 &normal)
1983 {
1984  double3 dir_v3v1 = v3 - v1;
1985  double3 tangent = math::cross(dir_v3v1, normal);
1986  double dot = math::dot(v1, tangent);
1987  return (math::dot(v4, tangent) >= dot) || (math::dot(v2, tangent) <= dot);
1988 }
1989 
2002 static Array<Face *> polyfill_triangulate_poly(Face *f, IMeshArena *arena)
2003 {
2004  /* Similar to loop body in #BM_mesh_calc_tessellation. */
2005  int flen = f->size();
2006  BLI_assert(flen >= 4);
2007  if (!f->plane_populated()) {
2008  f->populate_plane(false);
2009  }
2010  const double3 &poly_normal = f->plane->norm;
2011  float no[3] = {float(poly_normal[0]), float(poly_normal[1]), float(poly_normal[2])};
2012  normalize_v3(no);
2013  if (flen == 4) {
2014  const Vert *v0 = (*f)[0];
2015  const Vert *v1 = (*f)[1];
2016  const Vert *v2 = (*f)[2];
2017  const Vert *v3 = (*f)[3];
2018  int eo_01 = f->edge_orig[0];
2019  int eo_12 = f->edge_orig[1];
2020  int eo_23 = f->edge_orig[2];
2021  int eo_30 = f->edge_orig[3];
2022  Face *f0, *f1;
2023  if (UNLIKELY(is_quad_flip_first_third(v0->co, v1->co, v2->co, v3->co, f->plane->norm))) {
2024  f0 = arena->add_face({v0, v1, v3}, f->orig, {eo_01, -1, eo_30}, {false, false, false});
2025  f1 = arena->add_face({v1, v2, v3}, f->orig, {eo_12, eo_23, -1}, {false, false, false});
2026  }
2027  else {
2028  f0 = arena->add_face({v0, v1, v2}, f->orig, {eo_01, eo_12, -1}, {false, false, false});
2029  f1 = arena->add_face({v0, v2, v3}, f->orig, {-1, eo_23, eo_30}, {false, false, false});
2030  }
2031  return Array<Face *>{f0, f1};
2032  }
2033  /* Project along negative face normal so (x,y) can be used in 2d. */ float axis_mat[3][3];
2034  float(*projverts)[2];
2035  unsigned int(*tris)[3];
2036  const int totfilltri = flen - 2;
2037  /* Prepare projected vertices and array to receive triangles in tessellation. */
2038  tris = static_cast<unsigned int(*)[3]>(MEM_malloc_arrayN(totfilltri, sizeof(*tris), __func__));
2039  projverts = static_cast<float(*)[2]>(MEM_malloc_arrayN(flen, sizeof(*projverts), __func__));
2040  axis_dominant_v3_to_m3_negate(axis_mat, no);
2041  for (int j = 0; j < flen; ++j) {
2042  const double3 &dco = (*f)[j]->co;
2043  float co[3] = {float(dco[0]), float(dco[1]), float(dco[2])};
2044  mul_v2_m3v3(projverts[j], axis_mat, co);
2045  }
2046  BLI_polyfill_calc(projverts, flen, 1, tris);
2047  /* Put tessellation triangles into Face form. Record original edges where they exist. */
2048  Array<Face *> ans(totfilltri);
2049  for (int t = 0; t < totfilltri; ++t) {
2050  unsigned int *tri = tris[t];
2051  int eo[3];
2052  const Vert *v[3];
2053  for (int k = 0; k < 3; k++) {
2054  BLI_assert(tri[k] < flen);
2055  v[k] = (*f)[tri[k]];
2056  /* If tri edge goes between two successive indices in
2057  * the original face, then it is an original edge. */
2058  if ((tri[k] + 1) % flen == tri[(k + 1) % 3]) {
2059  eo[k] = f->edge_orig[tri[k]];
2060  }
2061  else {
2062  eo[k] = NO_INDEX;
2063  }
2064  ans[t] = arena->add_face(
2065  {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false});
2066  }
2067  }
2068 
2069  MEM_freeN(tris);
2070  MEM_freeN(projverts);
2071 
2072  return ans;
2073 }
2074 
2090 static Array<Face *> exact_triangulate_poly(Face *f, IMeshArena *arena)
2091 {
2092  int flen = f->size();
2093  CDT_input<mpq_class> cdt_in;
2094  cdt_in.vert = Array<mpq2>(flen);
2095  cdt_in.face = Array<Vector<int>>(1);
2096  cdt_in.face[0].reserve(flen);
2097  for (int i : f->index_range()) {
2098  cdt_in.face[0].append(i);
2099  }
2100  /* Project poly along dominant axis of normal to get 2d coords. */
2101  if (!f->plane_populated()) {
2102  f->populate_plane(false);
2103  }
2104  const double3 &poly_normal = f->plane->norm;
2105  int axis = math::dominant_axis(poly_normal);
2106  /* If project down y axis as opposed to x or z, the orientation
2107  * of the polygon will be reversed.
2108  * Yet another reversal happens if the poly normal in the dominant
2109  * direction is opposite that of the positive dominant axis. */
2110  bool rev1 = (axis == 1);
2111  bool rev2 = poly_normal[axis] < 0;
2112  bool rev = rev1 ^ rev2;
2113  for (int i = 0; i < flen; ++i) {
2114  int ii = rev ? flen - i - 1 : i;
2115  mpq2 &p2d = cdt_in.vert[ii];
2116  int k = 0;
2117  for (int j = 0; j < 3; ++j) {
2118  if (j != axis) {
2119  p2d[k++] = (*f)[ii]->co_exact[j];
2120  }
2121  }
2122  }
2124  int n_tris = cdt_out.face.size();
2125  Array<Face *> ans(n_tris);
2126  for (int t = 0; t < n_tris; ++t) {
2127  int i_v_out[3];
2128  const Vert *v[3];
2129  int eo[3];
2130  bool needs_steiner = false;
2131  for (int i = 0; i < 3; ++i) {
2132  i_v_out[i] = cdt_out.face[t][i];
2133  if (cdt_out.vert_orig[i_v_out[i]].size() == 0) {
2134  needs_steiner = true;
2135  break;
2136  }
2137  v[i] = (*f)[cdt_out.vert_orig[i_v_out[i]][0]];
2138  }
2139  if (needs_steiner) {
2140  /* Fall back on the polyfill triangulator. */
2141  return polyfill_triangulate_poly(f, arena);
2142  }
2143  Map<std::pair<int, int>, int> verts_to_edge;
2144  populate_cdt_edge_map(verts_to_edge, cdt_out);
2145  int foff = cdt_out.face_edge_offset;
2146  for (int i = 0; i < 3; ++i) {
2147  std::pair<int, int> vpair(i_v_out[i], i_v_out[(i + 1) % 3]);
2148  std::pair<int, int> vpair_canon = sorted_int_pair(vpair);
2149  int e_out = verts_to_edge.lookup_default(vpair_canon, NO_INDEX);
2150  BLI_assert(e_out != NO_INDEX);
2151  eo[i] = NO_INDEX;
2152  for (int orig : cdt_out.edge_orig[e_out]) {
2153  if (orig >= foff) {
2154  int pos = orig % foff;
2155  BLI_assert(pos < f->size());
2156  eo[i] = f->edge_orig[pos];
2157  break;
2158  }
2159  }
2160  }
2161  if (rev) {
2162  ans[t] = arena->add_face(
2163  {v[0], v[2], v[1]}, f->orig, {eo[2], eo[1], eo[0]}, {false, false, false});
2164  }
2165  else {
2166  ans[t] = arena->add_face(
2167  {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false});
2168  }
2169  }
2170  return ans;
2171 }
2172 
2173 static bool face_is_degenerate(const Face *f)
2174 {
2175  const Face &face = *f;
2176  const Vert *v0 = face[0];
2177  const Vert *v1 = face[1];
2178  const Vert *v2 = face[2];
2179  if (v0 == v1 || v0 == v2 || v1 == v2) {
2180  return true;
2181  }
2182  double3 da = v2->co - v0->co;
2183  double3 db = v2->co - v1->co;
2184  double3 dab = math::cross(da, db);
2185  double dab_length_squared = math::length_squared(dab);
2186  double err_bound = supremum_dot_cross(dab, dab) * index_dot_cross * DBL_EPSILON;
2187  if (dab_length_squared > err_bound) {
2188  return false;
2189  }
2190  mpq3 a = v2->co_exact - v0->co_exact;
2191  mpq3 b = v2->co_exact - v1->co_exact;
2192  mpq3 ab = math::cross(a, b);
2193  if (ab.x == 0 && ab.y == 0 && ab.z == 0) {
2194  return true;
2195  }
2196 
2197  return false;
2198 }
2199 
2201 static bool any_degenerate_tris_fast(const Array<Face *> triangulation)
2202 {
2203  for (const Face *f : triangulation) {
2204  const Vert *v0 = (*f)[0];
2205  const Vert *v1 = (*f)[1];
2206  const Vert *v2 = (*f)[2];
2207  if (v0 == v1 || v0 == v2 || v1 == v2) {
2208  return true;
2209  }
2210  double3 da = v2->co - v0->co;
2211  double3 db = v2->co - v1->co;
2212  double da_length_squared = math::length_squared(da);
2213  double db_length_squared = math::length_squared(db);
2214  if (da_length_squared == 0.0 || db_length_squared == 0.0) {
2215  return true;
2216  }
2217  /* |da x db| = |da| |db| sin t, where t is angle between them.
2218  * The triangle is almost degenerate if sin t is almost 0.
2219  * sin^2 t = |da x db|^2 / (|da|^2 |db|^2)
2220  */
2221  double3 dab = math::cross(da, db);
2222  double dab_length_squared = math::length_squared(dab);
2223  double sin_squared_t = dab_length_squared / (da_length_squared * db_length_squared);
2224  if (sin_squared_t < 1e-8) {
2225  return true;
2226  }
2227  }
2228  return false;
2229 }
2230 
2239 static Array<Face *> triangulate_poly(Face *f, IMeshArena *arena)
2240 {
2241  /* Try the much faster method using Blender's BLI_polyfill_calc. */
2242  Array<Face *> ans = polyfill_triangulate_poly(f, arena);
2243 
2244  /* This may create degenerate triangles. If so, try the exact CDT-based triangulator. */
2245  if (any_degenerate_tris_fast(ans)) {
2246  return exact_triangulate_poly(f, arena);
2247  }
2248  return ans;
2249 }
2250 
2251 IMesh triangulate_polymesh(IMesh &imesh, IMeshArena *arena)
2252 {
2253  Vector<Face *> face_tris;
2254  constexpr int estimated_tris_per_face = 3;
2255  face_tris.reserve(estimated_tris_per_face * imesh.face_size());
2256  threading::parallel_for(imesh.face_index_range(), 2048, [&](IndexRange range) {
2257  for (int i : range) {
2258  Face *f = imesh.face(i);
2259  if (!f->plane_populated() && f->size() >= 4) {
2260  f->populate_plane(false);
2261  }
2262  }
2263  });
2264  for (Face *f : imesh.faces()) {
2265  /* Tessellate face f, following plan similar to #BM_face_calc_tessellation. */
2266  int flen = f->size();
2267  if (flen == 3) {
2268  face_tris.append(f);
2269  }
2270  else {
2271  Array<Face *> tris = triangulate_poly(f, arena);
2272  for (Face *tri : tris) {
2273  face_tris.append(tri);
2274  }
2275  }
2276  }
2277  return IMesh(face_tris);
2278 }
2279 
2284 static IMesh extract_subdivided_tri(const CDT_data &cd,
2285  const IMesh &in_tm,
2286  int t,
2287  IMeshArena *arena)
2288 {
2289  const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
2290  int t_in_cdt = -1;
2291  for (int i = 0; i < cd.input_face.size(); ++i) {
2292  if (cd.input_face[i] == t) {
2293  t_in_cdt = i;
2294  }
2295  }
2296  if (t_in_cdt == -1) {
2297  std::cout << "Could not find " << t << " in cdt input tris\n";
2298  BLI_assert(false);
2299  return IMesh();
2300  }
2301  constexpr int inline_buf_size = 20;
2302  Vector<Face *, inline_buf_size> faces;
2303  for (int f : cdt_out.face.index_range()) {
2304  if (cdt_out.face_orig[f].contains(t_in_cdt)) {
2305  Face *facep = cdt_tri_as_imesh_face(f, t_in_cdt, cd, in_tm, arena);
2306  faces.append(facep);
2307  }
2308  }
2309  return IMesh(faces);
2310 }
2311 
2312 static bool bvhtreeverlap_cmp(const BVHTreeOverlap &a, const BVHTreeOverlap &b)
2313 {
2314  if (a.indexA < b.indexA) {
2315  return true;
2316  }
2317  if ((a.indexA == b.indexA) & (a.indexB < b.indexB)) {
2318  return true;
2319  }
2320  return false;
2321 }
2322 class TriOverlaps {
2323  BVHTree *tree_{nullptr};
2324  BVHTree *tree_b_{nullptr};
2325  BVHTreeOverlap *overlap_{nullptr};
2326  Array<int> first_overlap_;
2327  uint overlap_num_{0};
2328 
2329  struct CBData {
2330  const IMesh &tm;
2331  std::function<int(int)> shape_fn;
2332  int nshapes;
2333  bool use_self;
2334  };
2335 
2336  public:
2337  TriOverlaps(const IMesh &tm,
2338  const Array<BoundingBox> &tri_bb,
2339  int nshapes,
2340  std::function<int(int)> shape_fn,
2341  bool use_self)
2342  {
2343  constexpr int dbg_level = 0;
2344  if (dbg_level > 0) {
2345  std::cout << "TriOverlaps construction\n";
2346  }
2347  /* Tree type is 8 => octree; axis = 6 => using XYZ axes only. */
2348  tree_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6);
2349  /* In the common case of a binary boolean and no self intersection in
2350  * each shape, we will use two trees and simple bounding box overlap. */
2351  bool two_trees_no_self = nshapes == 2 && !use_self;
2352  if (two_trees_no_self) {
2353  tree_b_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6);
2354  }
2355 
2356  /* Create a Vector containing face shape. */
2357  Vector<int> shapes;
2358  shapes.resize(tm.face_size());
2359  threading::parallel_for(tm.face_index_range(), 2048, [&](IndexRange range) {
2360  for (int t : range) {
2361  shapes[t] = shape_fn(tm.face(t)->orig);
2362  }
2363  });
2364 
2365  float bbpts[6];
2366  for (int t : tm.face_index_range()) {
2367  const BoundingBox &bb = tri_bb[t];
2368  copy_v3_v3(bbpts, bb.min);
2369  copy_v3_v3(bbpts + 3, bb.max);
2370  int shape = shapes[t];
2371  if (two_trees_no_self) {
2372  if (shape == 0) {
2373  BLI_bvhtree_insert(tree_, t, bbpts, 2);
2374  }
2375  else if (shape == 1) {
2376  BLI_bvhtree_insert(tree_b_, t, bbpts, 2);
2377  }
2378  }
2379  else {
2380  if (shape != -1) {
2381  BLI_bvhtree_insert(tree_, t, bbpts, 2);
2382  }
2383  }
2384  }
2385  BLI_bvhtree_balance(tree_);
2386  if (two_trees_no_self) {
2387  BLI_bvhtree_balance(tree_b_);
2388  /* Don't expect a lot of trivial intersects in this case. */
2389  overlap_ = BLI_bvhtree_overlap(tree_, tree_b_, &overlap_num_, nullptr, nullptr);
2390  }
2391  else {
2392  CBData cbdata{tm, shape_fn, nshapes, use_self};
2393  if (nshapes == 1) {
2394  overlap_ = BLI_bvhtree_overlap(tree_, tree_, &overlap_num_, nullptr, nullptr);
2395  }
2396  else {
2397  overlap_ = BLI_bvhtree_overlap(
2398  tree_, tree_, &overlap_num_, only_different_shapes, &cbdata);
2399  }
2400  }
2401  /* The rest of the code is simpler and easier to parallelize if, in the two-trees case,
2402  * we repeat the overlaps with indexA and indexB reversed. It is important that
2403  * in the repeated part, sorting will then bring things with indexB together. */
2404  if (two_trees_no_self) {
2405  overlap_ = static_cast<BVHTreeOverlap *>(
2406  MEM_reallocN(overlap_, 2 * overlap_num_ * sizeof(overlap_[0])));
2407  for (uint i = 0; i < overlap_num_; ++i) {
2408  overlap_[overlap_num_ + i].indexA = overlap_[i].indexB;
2409  overlap_[overlap_num_ + i].indexB = overlap_[i].indexA;
2410  }
2411  overlap_num_ += overlap_num_;
2412  }
2413  /* Sort the overlaps to bring all the intersects with a given indexA together. */
2414  std::sort(overlap_, overlap_ + overlap_num_, bvhtreeverlap_cmp);
2415  if (dbg_level > 0) {
2416  std::cout << overlap_num_ << " overlaps found:\n";
2417  for (BVHTreeOverlap ov : overlap()) {
2418  std::cout << "A: " << ov.indexA << ", B: " << ov.indexB << "\n";
2419  }
2420  }
2421  first_overlap_ = Array<int>(tm.face_size(), -1);
2422  for (int i = 0; i < static_cast<int>(overlap_num_); ++i) {
2423  int t = overlap_[i].indexA;
2424  if (first_overlap_[t] == -1) {
2425  first_overlap_[t] = i;
2426  }
2427  }
2428  }
2429 
2430  ~TriOverlaps()
2431  {
2432  if (tree_) {
2433  BLI_bvhtree_free(tree_);
2434  }
2435  if (tree_b_) {
2436  BLI_bvhtree_free(tree_b_);
2437  }
2438  if (overlap_) {
2439  MEM_freeN(overlap_);
2440  }
2441  }
2442 
2443  Span<BVHTreeOverlap> overlap() const
2444  {
2445  return Span<BVHTreeOverlap>(overlap_, overlap_num_);
2446  }
2447 
2448  int first_overlap_index(int t) const
2449  {
2450  return first_overlap_[t];
2451  }
2452 
2453  private:
2454  static bool only_different_shapes(void *userdata, int index_a, int index_b, int UNUSED(thread))
2455  {
2456  CBData *cbdata = static_cast<CBData *>(userdata);
2457  return cbdata->tm.face(index_a)->orig != cbdata->tm.face(index_b)->orig;
2458  }
2459 };
2460 
2464 struct OverlapIttsData {
2465  Vector<std::pair<int, int>> intersect_pairs;
2466  Map<std::pair<int, int>, ITT_value> &itt_map;
2467  const IMesh &tm;
2468  IMeshArena *arena;
2469 
2470  OverlapIttsData(Map<std::pair<int, int>, ITT_value> &itt_map, const IMesh &tm, IMeshArena *arena)
2471  : itt_map(itt_map), tm(tm), arena(arena)
2472  {
2473  }
2474 };
2475 
2480 static std::pair<int, int> canon_int_pair(int a, int b)
2481 {
2482  if (a > b) {
2483  std::swap(a, b);
2484  }
2485  return std::pair<int, int>(a, b);
2486 }
2487 
2488 static void calc_overlap_itts_range_func(void *__restrict userdata,
2489  const int iter,
2490  const TaskParallelTLS *__restrict UNUSED(tls))
2491 {
2492  constexpr int dbg_level = 0;
2493  OverlapIttsData *data = static_cast<OverlapIttsData *>(userdata);
2494  std::pair<int, int> tri_pair = data->intersect_pairs[iter];
2495  int a = tri_pair.first;
2496  int b = tri_pair.second;
2497  if (dbg_level > 0) {
2498  std::cout << "calc_overlap_itts_range_func a=" << a << ", b=" << b << "\n";
2499  }
2500  ITT_value itt = intersect_tri_tri(data->tm, a, b);
2501  if (dbg_level > 0) {
2502  std::cout << "result of intersecting " << a << " and " << b << " = " << itt << "\n";
2503  }
2504  BLI_assert(data->itt_map.contains(tri_pair));
2505  data->itt_map.add_overwrite(tri_pair, itt);
2506 }
2507 
2512 static void calc_overlap_itts(Map<std::pair<int, int>, ITT_value> &itt_map,
2513  const IMesh &tm,
2514  const TriOverlaps &ov,
2515  IMeshArena *arena)
2516 {
2517  OverlapIttsData data(itt_map, tm, arena);
2518  /* Put dummy values in `itt_map` initially,
2519  * so map entries will exist when doing the range function.
2520  * This means we won't have to protect the `itt_map.add_overwrite` function with a lock. */
2521  for (const BVHTreeOverlap &olap : ov.overlap()) {
2522  std::pair<int, int> key = canon_int_pair(olap.indexA, olap.indexB);
2523  if (!itt_map.contains(key)) {
2524  itt_map.add_new(key, ITT_value());
2525  data.intersect_pairs.append(key);
2526  }
2527  }
2528  int tot_intersect_pairs = data.intersect_pairs.size();
2529  TaskParallelSettings settings;
2531  settings.min_iter_per_thread = 1000;
2532  settings.use_threading = intersect_use_threading;
2533  BLI_task_parallel_range(0, tot_intersect_pairs, &data, calc_overlap_itts_range_func, &settings);
2534 }
2535 
2542 static void calc_subdivided_non_cluster_tris(Array<IMesh> &r_tri_subdivided,
2543  const IMesh &tm,
2544  const Map<std::pair<int, int>, ITT_value> &itt_map,
2545  const CoplanarClusterInfo &clinfo,
2546  const TriOverlaps &ov,
2547  IMeshArena *arena)
2548 {
2549  const int dbg_level = 0;
2550  if (dbg_level > 0) {
2551  std::cout << "\nCALC_SUBDIVIDED_TRIS\n\n";
2552  }
2553  Span<BVHTreeOverlap> overlap = ov.overlap();
2554  struct OverlapTriRange {
2555  int tri_index;
2556  int overlap_start;
2557  int len;
2558  };
2559  Vector<OverlapTriRange> overlap_tri_range;
2560  int overlap_num = overlap.size();
2561  overlap_tri_range.reserve(overlap_num);
2562  int overlap_index = 0;
2563  while (overlap_index < overlap_num) {
2564  int t = overlap[overlap_index].indexA;
2565  int i = overlap_index;
2566  while (i + 1 < overlap_num && overlap[i + 1].indexA == t) {
2567  ++i;
2568  }
2569  /* Now overlap[overlap_index] to overlap[i] have indexA == t.
2570  * We only record ranges for triangles that are not in clusters,
2571  * because the ones in clusters are handled separately. */
2572  if (clinfo.tri_cluster(t) == NO_INDEX) {
2573  int len = i - overlap_index + 1;
2574  if (!(len == 1 && overlap[overlap_index].indexB == t)) {
2575  OverlapTriRange range = {t, overlap_index, len};
2576  overlap_tri_range.append(range);
2577 # ifdef PERFDEBUG
2578  bumpperfcount(0, len); /* Non-cluster overlaps. */
2579 # endif
2580  }
2581  }
2582  overlap_index = i + 1;
2583  }
2584  int overlap_tri_range_num = overlap_tri_range.size();
2585  Array<CDT_data> cd_data(overlap_tri_range_num);
2586  int grain_size = 64;
2587  threading::parallel_for(overlap_tri_range.index_range(), grain_size, [&](IndexRange range) {
2588  for (int otr_index : range) {
2589  OverlapTriRange &otr = overlap_tri_range[otr_index];
2590  int t = otr.tri_index;
2591  if (dbg_level > 0) {
2592  std::cout << "handling overlap range\nt=" << t << " start=" << otr.overlap_start
2593  << " len=" << otr.len << "\n";
2594  }
2595  constexpr int inline_capacity = 100;
2596  Vector<ITT_value, inline_capacity> itts(otr.len);
2597  for (int j = otr.overlap_start; j < otr.overlap_start + otr.len; ++j) {
2598  int t_other = overlap[j].indexB;
2599  std::pair<int, int> key = canon_int_pair(t, t_other);
2600  ITT_value itt;
2601  if (itt_map.contains(key)) {
2602  itt = itt_map.lookup(key);
2603  }
2604  if (itt.kind != INONE) {
2605  itts.append(itt);
2606  }
2607  if (dbg_level > 0) {
2608  std::cout << " tri t" << t_other << "; result = " << itt << "\n";
2609  }
2610  }
2611  if (itts.size() > 0) {
2612  cd_data[otr_index] = prepare_cdt_input(tm, t, itts);
2613  do_cdt(cd_data[otr_index]);
2614  }
2615  }
2616  });
2617  /* Extract the new faces serially, so that Boolean is repeatable regardless of parallelism. */
2618  for (int otr_index : overlap_tri_range.index_range()) {
2619  CDT_data &cdd = cd_data[otr_index];
2620  if (cdd.vert.size() > 0) {
2621  int t = overlap_tri_range[otr_index].tri_index;
2622  r_tri_subdivided[t] = extract_subdivided_tri(cdd, tm, t, arena);
2623  if (dbg_level > 1) {
2624  std::cout << "subdivide output for tri " << t << " = " << r_tri_subdivided[t];
2625  }
2626  }
2627  }
2628  /* Now have to put in the triangles that are the same as the input ones, and not in clusters.
2629  */
2630  threading::parallel_for(tm.face_index_range(), 2048, [&](IndexRange range) {
2631  for (int t : range) {
2632  if (r_tri_subdivided[t].face_size() == 0 && clinfo.tri_cluster(t) == NO_INDEX) {
2633  r_tri_subdivided[t] = IMesh({tm.face(t)});
2634  }
2635  }
2636  });
2637 }
2638 
2646 static void calc_cluster_tris(Array<IMesh> &tri_subdivided,
2647  const IMesh &tm,
2648  const CoplanarClusterInfo &clinfo,
2649  const Array<CDT_data> &cluster_subdivided,
2650  IMeshArena *arena)
2651 {
2652  for (int c : clinfo.index_range()) {
2653  const CoplanarCluster &cl = clinfo.cluster(c);
2654  const CDT_data &cd = cluster_subdivided[c];
2655  /* Each triangle in cluster c should be an input triangle in cd.input_faces.
2656  * (See prepare_cdt_input_for_cluster.)
2657  * So accumulate a Vector of Face* for each input face by going through the
2658  * output faces and making a Face for each input face that it is part of.
2659  * (The Boolean algorithm wants duplicates if a given output triangle is part
2660  * of more than one input triangle.)
2661  */
2662  int n_cluster_tris = cl.tot_tri();
2663  const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
2664  BLI_assert(cd.input_face.size() == n_cluster_tris);
2665  Array<Vector<Face *>> face_vec(n_cluster_tris);
2666  for (int cdt_out_t : cdt_out.face.index_range()) {
2667  for (int cdt_in_t : cdt_out.face_orig[cdt_out_t]) {
2668  Face *f = cdt_tri_as_imesh_face(cdt_out_t, cdt_in_t, cd, tm, arena);
2669  face_vec[cdt_in_t].append(f);
2670  }
2671  }
2672  for (int cdt_in_t : cd.input_face.index_range()) {
2673  int tm_t = cd.input_face[cdt_in_t];
2674  BLI_assert(tri_subdivided[tm_t].face_size() == 0);
2675  tri_subdivided[tm_t] = IMesh(face_vec[cdt_in_t]);
2676  }
2677  }
2678 }
2679 
2680 static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo,
2681  int c,
2682  const IMesh &tm,
2683  const TriOverlaps &ov,
2684  const Map<std::pair<int, int>, ITT_value> &itt_map,
2685  IMeshArena *UNUSED(arena))
2686 {
2687  constexpr int dbg_level = 0;
2688  BLI_assert(c < clinfo.tot_cluster());
2689  const CoplanarCluster &cl = clinfo.cluster(c);
2690  /* Make a CDT input with triangles from C and intersects from other triangles in tm. */
2691  if (dbg_level > 0) {
2692  std::cout << "CALC_CLUSTER_SUBDIVIDED for cluster " << c << " = " << cl << "\n";
2693  }
2694  /* Get vector itts of all intersections of a triangle of cl with any triangle of tm not
2695  * in cl and not co-planar with it (for that latter, if there were an intersection,
2696  * it should already be in cluster cl). */
2697  Vector<ITT_value> itts;
2698  Span<BVHTreeOverlap> ovspan = ov.overlap();
2699  for (int t : cl) {
2700  if (dbg_level > 0) {
2701  std::cout << "find intersects with triangle " << t << " of cluster\n";
2702  }
2703  int first_i = ov.first_overlap_index(t);
2704  if (first_i == -1) {
2705  continue;
2706  }
2707  for (int i = first_i; i < ovspan.size() && ovspan[i].indexA == t; ++i) {
2708  int t_other = ovspan[i].indexB;
2709  if (clinfo.tri_cluster(t_other) != c) {
2710  if (dbg_level > 0) {
2711  std::cout << "use intersect(" << t << "," << t_other << "\n";
2712  }
2713  std::pair<int, int> key = canon_int_pair(t, t_other);
2714  if (itt_map.contains(key)) {
2715  ITT_value itt = itt_map.lookup(key);
2716  if (!ELEM(itt.kind, INONE, ICOPLANAR)) {
2717  itts.append(itt);
2718  if (dbg_level > 0) {
2719  std::cout << " itt = " << itt << "\n";
2720  }
2721  }
2722  }
2723  }
2724  }
2725  }
2726  /* Use CDT to subdivide the cluster triangles and the points and segments in itts. */
2727  CDT_data cd_data = prepare_cdt_input_for_cluster(tm, clinfo, c, itts);
2728  do_cdt(cd_data);
2729  return cd_data;
2730 }
2731 
2732 static IMesh union_tri_subdivides(const blender::Array<IMesh> &tri_subdivided)
2733 {
2734  int tot_tri = 0;
2735  for (const IMesh &m : tri_subdivided) {
2736  tot_tri += m.face_size();
2737  }
2738  Array<Face *> faces(tot_tri);
2739  int face_index = 0;
2740  for (const IMesh &m : tri_subdivided) {
2741  for (Face *f : m.faces()) {
2742  faces[face_index++] = f;
2743  }
2744  }
2745  return IMesh(faces);
2746 }
2747 
2748 static CoplanarClusterInfo find_clusters(const IMesh &tm,
2749  const Array<BoundingBox> &tri_bb,
2750  const Map<std::pair<int, int>, ITT_value> &itt_map)
2751 {
2752  constexpr int dbg_level = 0;
2753  if (dbg_level > 0) {
2754  std::cout << "FIND_CLUSTERS\n";
2755  }
2756  CoplanarClusterInfo ans(tm.face_size());
2757  /* Use a VectorSet to get stable order from run to run. */
2758  VectorSet<int> maybe_coplanar_tris;
2759  maybe_coplanar_tris.reserve(2 * itt_map.size());
2760  for (auto item : itt_map.items()) {
2761  if (item.value.kind == ICOPLANAR) {
2762  int t1 = item.key.first;
2763  int t2 = item.key.second;
2764  maybe_coplanar_tris.add_multiple({t1, t2});
2765  }
2766  }
2767  if (dbg_level > 0) {
2768  std::cout << "found " << maybe_coplanar_tris.size() << " possible coplanar tris\n";
2769  }
2770  if (maybe_coplanar_tris.size() == 0) {
2771  if (dbg_level > 0) {
2772  std::cout << "No possible coplanar tris, so no clusters\n";
2773  }
2774  return ans;
2775  }
2776  /* There can be more than one #CoplanarCluster per plane. Accumulate them in
2777  * a Vector. We will have to merge some elements of the Vector as we discover
2778  * triangles that form intersection bridges between two or more clusters. */
2779  Map<Plane, Vector<CoplanarCluster>> plane_cls;
2780  plane_cls.reserve(maybe_coplanar_tris.size());
2781  for (int t : maybe_coplanar_tris) {
2782  /* Use a canonical version of the plane for map index.
2783  * We can't just store the canonical version in the face
2784  * since canonicalizing loses the orientation of the normal. */
2785  Plane tplane = *tm.face(t)->plane;
2786  BLI_assert(tplane.exact_populated());
2787  tplane.make_canonical();
2788  if (dbg_level > 0) {
2789  std::cout << "plane for tri " << t << " = " << &tplane << "\n";
2790  }
2791  /* Assume all planes are in canonical from (see canon_plane()). */
2792  if (plane_cls.contains(tplane)) {
2793  Vector<CoplanarCluster> &curcls = plane_cls.lookup(tplane);
2794  if (dbg_level > 0) {
2795  std::cout << "already has " << curcls.size() << " clusters\n";
2796  }
2797  /* Partition `curcls` into those that intersect t non-trivially, and those that don't. */
2798  Vector<CoplanarCluster *> int_cls;
2799  Vector<CoplanarCluster *> no_int_cls;
2800  for (CoplanarCluster &cl : curcls) {
2801  if (dbg_level > 1) {
2802  std::cout << "consider intersecting with cluster " << cl << "\n";
2803  }
2804  if (bbs_might_intersect(tri_bb[t], cl.bounding_box())) {
2805  if (dbg_level > 1) {
2806  std::cout << "append to int_cls\n";
2807  }
2808  int_cls.append(&cl);
2809  }
2810  else {
2811  if (dbg_level > 1) {
2812  std::cout << "append to no_int_cls\n";
2813  }
2814  no_int_cls.append(&cl);
2815  }
2816  }
2817  if (int_cls.size() == 0) {
2818  /* t doesn't intersect any existing cluster in its plane, so make one just for it. */
2819  if (dbg_level > 1) {
2820  std::cout << "no intersecting clusters for t, make a new one\n";
2821  }
2822  curcls.append(CoplanarCluster(t, tri_bb[t]));
2823  }
2824  else if (int_cls.size() == 1) {
2825  /* t intersects exactly one existing cluster, so can add t to that cluster. */
2826  if (dbg_level > 1) {
2827  std::cout << "exactly one existing cluster, " << int_cls[0] << ", adding to it\n";
2828  }
2829  int_cls[0]->add_tri(t, tri_bb[t]);
2830  }
2831  else {
2832  /* t intersections 2 or more existing clusters: need to merge them and replace all the
2833  * originals with the merged one in `curcls`. */
2834  if (dbg_level > 1) {
2835  std::cout << "merging\n";
2836  }
2837  CoplanarCluster mergecl;
2838  mergecl.add_tri(t, tri_bb[t]);
2839  for (CoplanarCluster *cl : int_cls) {
2840  for (int t : *cl) {
2841  mergecl.add_tri(t, tri_bb[t]);
2842  }
2843  }
2844  Vector<CoplanarCluster> newvec;
2845  newvec.append(mergecl);
2846  for (CoplanarCluster *cl_no_int : no_int_cls) {
2847  newvec.append(*cl_no_int);
2848  }
2849  plane_cls.add_overwrite(tplane, newvec);
2850  }
2851  }
2852  else {
2853  if (dbg_level > 0) {
2854  std::cout << "first cluster for its plane\n";
2855  }
2856  plane_cls.add_new(tplane, Vector<CoplanarCluster>{CoplanarCluster(t, tri_bb[t])});
2857  }
2858  }
2859  /* Does this give deterministic order for cluster ids? I think so, since
2860  * hash for planes is on their values, not their addresses. */
2861  for (auto item : plane_cls.items()) {
2862  for (const CoplanarCluster &cl : item.value) {
2863  if (cl.tot_tri() > 1) {
2864  ans.add_cluster(cl);
2865  }
2866  }
2867  }
2868 
2869  return ans;
2870 }
2871 
2872 /* Data and functions to test triangle degeneracy in parallel. */
2873 struct DegenData {
2874  const IMesh &tm;
2875 };
2876 
2877 struct DegenChunkData {
2878  bool has_degenerate_tri = false;
2879 };
2880 
2881 static void degenerate_range_func(void *__restrict userdata,
2882  const int iter,
2883  const TaskParallelTLS *__restrict tls)
2884 {
2885  DegenData *data = static_cast<DegenData *>(userdata);
2886  DegenChunkData *chunk_data = static_cast<DegenChunkData *>(tls->userdata_chunk);
2887  const Face *f = data->tm.face(iter);
2888  bool is_degenerate = face_is_degenerate(f);
2889  chunk_data->has_degenerate_tri |= is_degenerate;
2890 }
2891 
2892 static void degenerate_reduce(const void *__restrict UNUSED(userdata),
2893  void *__restrict chunk_join,
2894  void *__restrict chunk)
2895 {
2896  DegenChunkData *degen_chunk_join = static_cast<DegenChunkData *>(chunk_join);
2897  DegenChunkData *degen_chunk = static_cast<DegenChunkData *>(chunk);
2898  degen_chunk_join->has_degenerate_tri |= degen_chunk->has_degenerate_tri;
2899 }
2900 
2901 /* Does triangle #IMesh tm have any triangles with zero area? */
2902 static bool has_degenerate_tris(const IMesh &tm)
2903 {
2904  DegenData degen_data = {tm};
2905  DegenChunkData degen_chunk_data;
2906  TaskParallelSettings settings;
2908  settings.userdata_chunk = &degen_chunk_data;
2909  settings.userdata_chunk_size = sizeof(degen_chunk_data);
2910  settings.func_reduce = degenerate_reduce;
2911  settings.min_iter_per_thread = 1000;
2912  settings.use_threading = intersect_use_threading;
2913  BLI_task_parallel_range(0, tm.face_size(), &degen_data, degenerate_range_func, &settings);
2914  return degen_chunk_data.has_degenerate_tri;
2915 }
2916 
2917 static IMesh remove_degenerate_tris(const IMesh &tm_in)
2918 {
2919  IMesh ans;
2920  Vector<Face *> new_faces;
2921  new_faces.reserve(tm_in.face_size());
2922  for (Face *f : tm_in.faces()) {
2923  if (!face_is_degenerate(f)) {
2924  new_faces.append(f);
2925  }
2926  }
2927  ans.set_faces(new_faces);
2928  return ans;
2929 }
2930 
2931 IMesh trimesh_self_intersect(const IMesh &tm_in, IMeshArena *arena)
2932 {
2933  return trimesh_nary_intersect(
2934  tm_in, 1, [](int UNUSED(t)) { return 0; }, true, arena);
2935 }
2936 
2937 IMesh trimesh_nary_intersect(const IMesh &tm_in,
2938  int nshapes,
2939  std::function<int(int)> shape_fn,
2940  bool use_self,
2941  IMeshArena *arena)
2942 {
2943  constexpr int dbg_level = 0;
2944  if (dbg_level > 0) {
2945  std::cout << "\nTRIMESH_NARY_INTERSECT nshapes=" << nshapes << " use_self=" << use_self
2946  << "\n";
2947  for (const Face *f : tm_in.faces()) {
2948  BLI_assert(f->is_tri());
2949  UNUSED_VARS_NDEBUG(f);
2950  }
2951  if (dbg_level > 1) {
2952  std::cout << "input mesh:\n" << tm_in;
2953  for (int t : tm_in.face_index_range()) {
2954  std::cout << "shape(" << t << ") = " << shape_fn(tm_in.face(t)->orig) << "\n";
2955  }
2956  write_obj_mesh(const_cast<IMesh &>(tm_in), "trimesh_input");
2957  }
2958  }
2959 # ifdef PERFDEBUG
2960  perfdata_init();
2961  double start_time = PIL_check_seconds_timer();
2962  std::cout << "trimesh_nary_intersect start\n";
2963 # endif
2964  /* Usually can use tm_in but if it has degenerate or illegal triangles,
2965  * then need to work on a copy of it without those triangles. */
2966  const IMesh *tm_clean = &tm_in;
2967  IMesh tm_cleaned;
2968  if (has_degenerate_tris(tm_in)) {
2969  if (dbg_level > 0) {
2970  std::cout << "cleaning degenerate triangles\n";
2971  }
2972  tm_cleaned = remove_degenerate_tris(tm_in);
2973  tm_clean = &tm_cleaned;
2974  if (dbg_level > 1) {
2975  std::cout << "cleaned input mesh:\n" << tm_cleaned;
2976  }
2977  }
2978 # ifdef PERFDEBUG
2979  double clean_time = PIL_check_seconds_timer();
2980  std::cout << "cleaned, time = " << clean_time - start_time << "\n";
2981 # endif
2982  Array<BoundingBox> tri_bb = calc_face_bounding_boxes(*tm_clean);
2983 # ifdef PERFDEBUG
2984  double bb_calc_time = PIL_check_seconds_timer();
2985  std::cout << "bbs calculated, time = " << bb_calc_time - clean_time << "\n";
2986 # endif
2987  TriOverlaps tri_ov(*tm_clean, tri_bb, nshapes, shape_fn, use_self);
2988 # ifdef PERFDEBUG
2989  double overlap_time = PIL_check_seconds_timer();
2990  std::cout << "intersect overlaps calculated, time = " << overlap_time - bb_calc_time << "\n";
2991 # endif
2992  Array<IMesh> tri_subdivided(tm_clean->face_size(), NoInitialization());
2993  threading::parallel_for(tm_clean->face_index_range(), 1024, [&](IndexRange range) {
2994  for (int t : range) {
2995  if (tri_ov.first_overlap_index(t) != -1) {
2996  tm_clean->face(t)->populate_plane(true);
2997  }
2998  new (static_cast<void *>(&tri_subdivided[t])) IMesh;
2999  }
3000  });
3001 # ifdef PERFDEBUG
3002  double plane_populate = PIL_check_seconds_timer();
3003  std::cout << "planes populated, time = " << plane_populate - overlap_time << "\n";
3004 # endif
3005  /* itt_map((a,b)) will hold the intersection value resulting from intersecting
3006  * triangles with indices a and b, where a < b. */
3007  Map<std::pair<int, int>, ITT_value> itt_map;
3008  itt_map.reserve(tri_ov.overlap().size());
3009  calc_overlap_itts(itt_map, *tm_clean, tri_ov, arena);
3010 # ifdef PERFDEBUG
3011  double itt_time = PIL_check_seconds_timer();
3012  std::cout << "itts found, time = " << itt_time - plane_populate << "\n";
3013 # endif
3014  CoplanarClusterInfo clinfo = find_clusters(*tm_clean, tri_bb, itt_map);
3015  if (dbg_level > 1) {
3016  std::cout << clinfo;
3017  }
3018 # ifdef PERFDEBUG
3019  double find_cluster_time = PIL_check_seconds_timer();
3020  std::cout << "clusters found, time = " << find_cluster_time - itt_time << "\n";
3021  doperfmax(0, tm_in.face_size());
3022  doperfmax(1, clinfo.tot_cluster());
3023  doperfmax(2, tri_ov.overlap().size());
3024 # endif
3025  calc_subdivided_non_cluster_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena);
3026 # ifdef PERFDEBUG
3027  double subdivided_tris_time = PIL_check_seconds_timer();
3028  std::cout << "subdivided non-cluster tris found, time = " << subdivided_tris_time - itt_time
3029  << "\n";
3030 # endif
3031  Array<CDT_data> cluster_subdivided(clinfo.tot_cluster());
3032  for (int c : clinfo.index_range()) {
3033  cluster_subdivided[c] = calc_cluster_subdivided(clinfo, c, *tm_clean, tri_ov, itt_map, arena);
3034  }
3035 # ifdef PERFDEBUG
3036  double cluster_subdivide_time = PIL_check_seconds_timer();
3037  std::cout << "subdivided clusters found, time = "
3038  << cluster_subdivide_time - subdivided_tris_time << "\n";
3039 # endif
3040  calc_cluster_tris(tri_subdivided, *tm_clean, clinfo, cluster_subdivided, arena);
3041 # ifdef PERFDEBUG
3042  double extract_time = PIL_check_seconds_timer();
3043  std::cout << "subdivided cluster tris found, time = " << extract_time - cluster_subdivide_time
3044  << "\n";
3045 # endif
3046  IMesh combined = union_tri_subdivides(tri_subdivided);
3047  if (dbg_level > 1) {
3048  std::cout << "TRIMESH_NARY_INTERSECT answer:\n";
3049  std::cout << combined;
3050  }
3051 # ifdef PERFDEBUG
3052  double end_time = PIL_check_seconds_timer();
3053  std::cout << "triangles combined, time = " << end_time - extract_time << "\n";
3054  std::cout << "trimesh_nary_intersect done, total time = " << end_time - start_time << "\n";
3055  dump_perfdata();
3056 # endif
3057  return combined;
3058 }
3059 
3060 static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl)
3061 {
3062  os << "cl(";
3063  bool first = true;
3064  for (const int t : cl) {
3065  if (first) {
3066  first = false;
3067  }
3068  else {
3069  os << ",";
3070  }
3071  os << t;
3072  }
3073  os << ")";
3074  return os;
3075 }
3076 
3077 static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo)
3078 {
3079  os << "Coplanar Cluster Info:\n";
3080  for (int c : clinfo.index_range()) {
3081  os << c << ": " << clinfo.cluster(c) << "\n";
3082  }
3083  return os;
3084 }
3085 
3086 static std::ostream &operator<<(std::ostream &os, const ITT_value &itt)
3087 {
3088  switch (itt.kind) {
3089  case INONE:
3090  os << "none";
3091  break;
3092  case IPOINT:
3093  os << "point " << itt.p1;
3094  break;
3095  case ISEGMENT:
3096  os << "segment " << itt.p1 << " " << itt.p2;
3097  break;
3098  case ICOPLANAR:
3099  os << "co-planar t" << itt.t_source;
3100  break;
3101  }
3102  return os;
3103 }
3104 
3105 void write_obj_mesh(IMesh &m, const std::string &objname)
3106 {
3107  /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library.
3108  * This is just for developer debugging anyway,
3109  * and should never be called in production Blender. */
3110 # ifdef _WIN_32
3111  const char *objdir = BLI_getenv("HOME");
3112 # else
3113  const char *objdir = "/tmp/";
3114 # endif
3115  if (m.face_size() == 0) {
3116  return;
3117  }
3118 
3119  std::string fname = std::string(objdir) + objname + std::string(".obj");
3120  std::ofstream f;
3121  f.open(fname);
3122  if (!f) {
3123  std::cout << "Could not open file " << fname << "\n";
3124  return;
3125  }
3126 
3127  if (!m.has_verts()) {
3128  m.populate_vert();
3129  }
3130  for (const Vert *v : m.vertices()) {
3131  const double3 dv = v->co;
3132  f << "v " << dv[0] << " " << dv[1] << " " << dv[2] << "\n";
3133  }
3134  int i = 0;
3135  for (const Face *face : m.faces()) {
3136  /* OBJ files use 1-indexing for vertices. */
3137  f << "f ";
3138  for (const Vert *v : *face) {
3139  int i = m.lookup_vert(v);
3140  BLI_assert(i != NO_INDEX);
3141  /* OBJ files use 1-indexing for vertices. */
3142  f << i + 1 << " ";
3143  }
3144  f << "\n";
3145  ++i;
3146  }
3147  f.close();
3148 }
3149 
3150 # ifdef PERFDEBUG
3151 struct PerfCounts {
3152  Vector<int> count;
3153  Vector<const char *> count_name;
3154  Vector<int> max;
3155  Vector<const char *> max_name;
3156 };
3157 
3158 static PerfCounts *perfdata = nullptr;
3159 
3160 static void perfdata_init()
3161 {
3162  perfdata = new PerfCounts;
3163 
3164  /* count 0. */
3165  perfdata->count.append(0);
3166  perfdata->count_name.append("Non-cluster overlaps");
3167 
3168  /* count 1. */
3169  perfdata->count.append(0);
3170  perfdata->count_name.append("intersect_tri_tri calls");
3171 
3172  /* count 2. */
3173  perfdata->count.append(0);
3174  perfdata->count_name.append("tri tri intersects decided by filter plane tests");
3175 
3176  /* count 3. */
3177  perfdata->count.append(0);
3178  perfdata->count_name.append("tri tri intersects decided by exact plane tests");
3179 
3180  /* count 4. */
3181  perfdata->count.append(0);
3182  perfdata->count_name.append("final non-NONE intersects");
3183 
3184  /* max 0. */
3185  perfdata->max.append(0);
3186  perfdata->max_name.append("total faces");
3187 
3188  /* max 1. */
3189  perfdata->max.append(0);
3190  perfdata->max_name.append("total clusters");
3191 
3192  /* max 2. */
3193  perfdata->max.append(0);
3194  perfdata->max_name.append("total overlaps");
3195 }
3196 
3197 static void incperfcount(int countnum)
3198 {
3199  perfdata->count[countnum]++;
3200 }
3201 
3202 static void bumpperfcount(int countnum, int amt)
3203 {
3204  perfdata->count[countnum] += amt;
3205 }
3206 
3207 static void doperfmax(int maxnum, int val)
3208 {
3209  perfdata->max[maxnum] = max_ii(perfdata->max[maxnum], val);
3210 }
3211 
3212 static void dump_perfdata()
3213 {
3214  std::cout << "\nPERFDATA\n";
3215  for (int i : perfdata->count.index_range()) {
3216  std::cout << perfdata->count_name[i] << " = " << perfdata->count[i] << "\n";
3217  }
3218  for (int i : perfdata->max.index_range()) {
3219  std::cout << perfdata->max_name[i] << " = " << perfdata->max[i] << "\n";
3220  }
3221  delete perfdata;
3222 }
3223 # endif
3224 
3225 } // namespace blender::meshintersect
3226 
3227 #endif // WITH_GMP
typedef float(TangentPoint)[2]
#define BLI_assert(a)
Definition: BLI_assert.h:46
@ CDT_INSIDE
BVHTreeOverlap * BLI_bvhtree_overlap(const BVHTree *tree1, const BVHTree *tree2, unsigned int *r_overlap_num, BVHTree_OverlapCallback callback, void *userdata)
Definition: BLI_kdopbvh.c:1363
void BLI_bvhtree_balance(BVHTree *tree)
Definition: BLI_kdopbvh.c:937
BVHTree * BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
Definition: BLI_kdopbvh.c:854
void BLI_bvhtree_free(BVHTree *tree)
Definition: BLI_kdopbvh.c:926
void BLI_bvhtree_insert(BVHTree *tree, int index, const float co[3], int numpoints)
Definition: BLI_kdopbvh.c:979
MINLINE int max_ii(int a, int b)
MINLINE double max_dd(double a, double b)
Math vector functions needed specifically for mesh intersect and boolean.
bool isect_aabb_aabb_v3(const float min1[3], const float max1[3], const float min2[3], const float max2[3])
Definition: math_geom.c:3048
void axis_dominant_v3_to_m3_negate(float r_mat[3][3], const float normal[3])
Definition: math_geom.c:3544
void mul_v2_m3v3(float r[2], const float M[3][3], const float a[3])
Definition: math_matrix.c:917
MINLINE float normalize_v3(float r[3])
MINLINE void copy_v3_v3(float r[3], const float a[3])
const char * BLI_getenv(const char *env) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT
Definition: path_util.c:1168
void BLI_polyfill_calc(const float(*coords)[2], unsigned int coords_num, int coords_sign, unsigned int(*r_tris)[3])
Definition: polyfill_2d.c:875
unsigned int uint
Definition: BLI_sys_types.h:67
void BLI_task_parallel_range(int start, int stop, void *userdata, TaskParallelRangeFunc func, const TaskParallelSettings *settings)
Definition: task_range.cc:94
BLI_INLINE void BLI_parallel_range_settings_defaults(TaskParallelSettings *settings)
Definition: BLI_task.h:293
pthread_spinlock_t SpinLock
Definition: BLI_threads.h:110
void BLI_mutex_free(ThreadMutex *mutex)
Definition: threads.cc:400
ThreadMutex * BLI_mutex_alloc(void)
Definition: threads.cc:393
void BLI_mutex_lock(ThreadMutex *mutex)
Definition: threads.cc:373
void BLI_mutex_unlock(ThreadMutex *mutex)
Definition: threads.cc:378
void BLI_spin_init(SpinLock *spin)
Definition: threads.cc:419
void BLI_spin_unlock(SpinLock *spin)
Definition: threads.cc:452
void BLI_spin_lock(SpinLock *spin)
Definition: threads.cc:433
pthread_mutex_t ThreadMutex
Definition: BLI_threads.h:82
void BLI_spin_end(SpinLock *spin)
Definition: threads.cc:467
#define UNUSED_VARS_NDEBUG(...)
#define UNUSED(x)
#define UNLIKELY(x)
#define ELEM(...)
void swap(T &a, T &b)
Definition: Common.h:19
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble z
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint i1
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble y2 _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat y2 _GL_VOID_RET _GL_VOID GLint GLint GLint y2 _GL_VOID_RET _GL_VOID GLshort GLshort GLshort y2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLuint *buffer _GL_VOID_RET _GL_VOID GLdouble t _GL_VOID_RET _GL_VOID GLfloat t _GL_VOID_RET _GL_VOID GLint t _GL_VOID_RET _GL_VOID GLshort t _GL_VOID_RET _GL_VOID GLdouble t
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble v1
#define MEM_reallocN(vmemh, len)
in reality light always falls off quadratically Particle Retrieve the data of the particle that spawned the object for example to give variation to multiple instances of an object Point Retrieve information about points in a point cloud Retrieve the edges of an object as it appears to Cycles topology will always appear triangulated Convert a blackbody temperature to an RGB value Normal Map
Platform independent time functions.
int pad[32 - sizeof(int)]
ATTR_WARN_UNUSED_RESULT const BMVert * v2
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
ATTR_WARN_UNUSED_RESULT const BMVert * v
static DBVT_INLINE btScalar size(const btDbvtVolume &a)
Definition: btDbvt.cpp:52
void sort(btMatrix3x3 &U, btVector3 &sigma, btMatrix3x3 &V, int t)
Helper function of 3X3 SVD for sorting singular values.
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
Definition: thread.h:34
blender::meshintersect::CDT_result< double > delaunay_2d_calc(const CDT_input< double > &input, CDT_output_type output_type)
struct Vert Vert
int len
Definition: draw_manager.c:108
static float verts[][3]
uint pos
#define UINT_MAX
Definition: hash_md5.c:43
IconTextureDrawCall normal
int count
void *(* MEM_malloc_arrayN)(size_t len, size_t size, const char *str)
Definition: mallocn.c:34
void(* MEM_freeN)(void *vmemh)
Definition: mallocn.c:27
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
static char faces[256]
struct Face Face
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
GAttributeReader lookup(const void *owner, const AttributeIDRef &attribute_id)
std::ostream & operator<<(std::ostream &stream, const AssetCatalogPath &path_to_append)
bool operator==(const AttributeIDRef &a, const AttributeIDRef &b)
T dot(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
vec_base< T, 3 > cross(const vec_base< T, 3 > &a, const vec_base< T, 3 > &b)
vec_base< T, 3 > cross_poly(Span< vec_base< T, 3 >> poly)
T abs(const T &a)
T length_squared(const vec_base< T, Size > &a)
int dominant_axis(const vec_base< T, 3 > &a)
static meshintersect::CDT_result< double > do_cdt(const bke::CurvesGeometry &curves, const CDT_output_type output_type)
void parallel_for(IndexRange range, int64_t grain_size, const Function &function)
Definition: BLI_task.hh:51
vec_base< double, 3 > double3
static int sgn(double x)
void parallel_sort(RandomAccessIterator begin, RandomAccessIterator end)
Definition: BLI_sort.hh:21
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
#define hash
Definition: noise.c:153
__int64 int64_t
Definition: stdint.h:89
unsigned __int64 uint64_t
Definition: stdint.h:90
float co[3]
Definition: bmesh_class.h:87
Point3i begin
Definition: GeoCommon.h:23
TaskParallelReduceFunc func_reduce
Definition: BLI_task.h:181
size_t userdata_chunk_size
Definition: BLI_task.h:169
double PIL_check_seconds_timer(void)
Definition: time.c:64
float max