Blender  V3.3
mesh_boolean.cc
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 #ifdef WITH_GMP
8 
9 # include <algorithm>
10 # include <atomic>
11 # include <fstream>
12 # include <iostream>
13 
14 # include "BLI_array.hh"
15 # include "BLI_assert.h"
16 # include "BLI_delaunay_2d.h"
17 # include "BLI_hash.hh"
18 # include "BLI_kdopbvh.h"
19 # include "BLI_map.hh"
20 # include "BLI_math.h"
21 # include "BLI_math_boolean.hh"
22 # include "BLI_math_geom.h"
23 # include "BLI_math_mpq.hh"
24 # include "BLI_math_vec_mpq_types.hh"
25 # include "BLI_math_vector.hh"
26 # include "BLI_mesh_intersect.hh"
27 # include "BLI_set.hh"
28 # include "BLI_span.hh"
29 # include "BLI_stack.hh"
30 # include "BLI_task.hh"
31 # include "BLI_vector.hh"
32 # include "BLI_vector_set.hh"
33 
34 # include "PIL_time.h"
35 
36 # include "BLI_mesh_boolean.hh"
37 
38 # ifdef WITH_TBB
39 # include <tbb/parallel_reduce.h>
40 # include <tbb/spin_mutex.h>
41 # endif
42 
43 // # define PERFDEBUG
44 
45 namespace blender::meshintersect {
46 
52 class Edge {
53  const Vert *v_[2]{nullptr, nullptr};
54 
55  public:
56  Edge() = default;
57  Edge(const Vert *v0, const Vert *v1)
58  {
59  if (v0->id <= v1->id) {
60  v_[0] = v0;
61  v_[1] = v1;
62  }
63  else {
64  v_[0] = v1;
65  v_[1] = v0;
66  }
67  }
68 
69  const Vert *v0() const
70  {
71  return v_[0];
72  }
73 
74  const Vert *v1() const
75  {
76  return v_[1];
77  }
78 
79  const Vert *operator[](int i) const
80  {
81  return v_[i];
82  }
83 
84  bool operator==(Edge other) const
85  {
86  return v_[0]->id == other.v_[0]->id && v_[1]->id == other.v_[1]->id;
87  }
88 
89  uint64_t hash() const
90  {
91  return get_default_hash_2(v_[0]->id, v_[1]->id);
92  }
93 };
94 
95 static std::ostream &operator<<(std::ostream &os, const Edge &e)
96 {
97  if (e.v0() == nullptr) {
98  BLI_assert(e.v1() == nullptr);
99  os << "(null,null)";
100  }
101  else {
102  os << "(" << e.v0() << "," << e.v1() << ")";
103  }
104  return os;
105 }
106 
107 static std::ostream &operator<<(std::ostream &os, const Span<int> &a)
108 {
109  for (int i : a.index_range()) {
110  os << a[i];
111  if (i != a.size() - 1) {
112  os << " ";
113  }
114  }
115  return os;
116 }
117 
118 static std::ostream &operator<<(std::ostream &os, const Array<int> &iarr)
119 {
120  os << Span<int>(iarr);
121  return os;
122 }
123 
125 class TriMeshTopology : NonCopyable {
127  Map<Edge, Vector<int> *> edge_tri_;
129  Map<const Vert *, Vector<Edge>> vert_edges_;
130 
131  public:
132  TriMeshTopology(const IMesh &tm);
133  ~TriMeshTopology();
134 
135  /* If e is manifold, return index of the other triangle (not t) that has it.
136  * Else return NO_INDEX. */
137  int other_tri_if_manifold(Edge e, int t) const
138  {
139  const auto *p = edge_tri_.lookup_ptr(e);
140  if (p != nullptr && (*p)->size() == 2) {
141  return ((**p)[0] == t) ? (**p)[1] : (**p)[0];
142  }
143  return NO_INDEX;
144  }
145 
146  /* Which triangles share edge e (in either orientation)? */
147  const Vector<int> *edge_tris(Edge e) const
148  {
149  return edge_tri_.lookup_default(e, nullptr);
150  }
151 
152  /* Which edges are incident on the given vertex?
153  * We assume v has some incident edges. */
154  const Vector<Edge> &vert_edges(const Vert *v) const
155  {
156  return vert_edges_.lookup(v);
157  }
158 
159  Map<Edge, Vector<int> *>::ItemIterator edge_tri_map_items() const
160  {
161  return edge_tri_.items();
162  }
163 };
164 
165 TriMeshTopology::TriMeshTopology(const IMesh &tm)
166 {
167  const int dbg_level = 0;
168  if (dbg_level > 0) {
169  std::cout << "TRIMESHTOPOLOGY CONSTRUCTION\n";
170  }
171  /* If everything were manifold, `F+V-E=2` and `E=3F/2`.
172  * So an likely overestimate, allowing for non-manifoldness, is `E=2F` and `V=F`. */
173  const int estimate_num_edges = 2 * tm.face_size();
174  const int estimate_verts_num = tm.face_size();
175  edge_tri_.reserve(estimate_num_edges);
176  vert_edges_.reserve(estimate_verts_num);
177  for (int t : tm.face_index_range()) {
178  const Face &tri = *tm.face(t);
179  BLI_assert(tri.is_tri());
180  for (int i = 0; i < 3; ++i) {
181  const Vert *v = tri[i];
182  const Vert *vnext = tri[(i + 1) % 3];
183  Edge e(v, vnext);
184  Vector<Edge> *edges = vert_edges_.lookup_ptr(v);
185  if (edges == nullptr) {
186  vert_edges_.add_new(v, Vector<Edge>());
187  edges = vert_edges_.lookup_ptr(v);
188  BLI_assert(edges != nullptr);
189  }
190  edges->append_non_duplicates(e);
191 
192  auto *p = edge_tri_.lookup_ptr(Edge(v, vnext));
193  if (p == nullptr) {
194  edge_tri_.add_new(e, new Vector<int>{t});
195  }
196  else {
197  (*p)->append_non_duplicates(t);
198  }
199  }
200  }
201  /* Debugging. */
202  if (dbg_level > 0) {
203  std::cout << "After TriMeshTopology construction\n";
204  for (auto item : edge_tri_.items()) {
205  std::cout << "tris for edge " << item.key << ": " << *item.value << "\n";
206  constexpr bool print_stats = false;
207  if (print_stats) {
208  edge_tri_.print_stats();
209  }
210  }
211  for (auto item : vert_edges_.items()) {
212  std::cout << "edges for vert " << item.key << ":\n";
213  for (const Edge &e : item.value) {
214  std::cout << " " << e << "\n";
215  }
216  std::cout << "\n";
217  }
218  }
219 }
220 
221 TriMeshTopology::~TriMeshTopology()
222 {
223  Vector<Vector<int> *> values;
224 
225  /* Deconstructing is faster in parallel, so it is worth building an array of things to delete. */
226  for (auto *item : edge_tri_.values()) {
227  values.append(item);
228  }
229 
230  threading::parallel_for(values.index_range(), 256, [&](IndexRange range) {
231  for (int i : range) {
232  delete values[i];
233  }
234  });
235 }
236 
238 class Patch {
239  Vector<int> tri_; /* Indices of triangles in the Patch. */
240 
241  public:
242  Patch() = default;
243 
244  void add_tri(int t)
245  {
246  tri_.append(t);
247  }
248 
249  int tot_tri() const
250  {
251  return tri_.size();
252  }
253 
254  int tri(int i) const
255  {
256  return tri_[i];
257  }
258 
259  IndexRange tri_range() const
260  {
261  return IndexRange(tri_.size());
262  }
263 
264  Span<int> tris() const
265  {
266  return Span<int>(tri_);
267  }
268 
269  int cell_above{NO_INDEX};
270  int cell_below{NO_INDEX};
271  int component{NO_INDEX};
272 };
273 
274 static std::ostream &operator<<(std::ostream &os, const Patch &patch)
275 {
276  os << "Patch " << patch.tris();
277  if (patch.cell_above != NO_INDEX) {
278  os << " cell_above=" << patch.cell_above;
279  }
280  else {
281  os << " cell_above not set";
282  }
283  if (patch.cell_below != NO_INDEX) {
284  os << " cell_below=" << patch.cell_below;
285  }
286  else {
287  os << " cell_below not set";
288  }
289  return os;
290 }
291 
292 class PatchesInfo {
294  Vector<Patch> patch_;
296  Array<int> tri_patch_;
298  Map<std::pair<int, int>, Edge> pp_edge_;
299 
300  public:
301  explicit PatchesInfo(int ntri)
302  {
303  constexpr int max_expected_patch_patch_incidences = 100;
304  tri_patch_ = Array<int>(ntri, NO_INDEX);
305  pp_edge_.reserve(max_expected_patch_patch_incidences);
306  }
307 
308  int tri_patch(int t) const
309  {
310  return tri_patch_[t];
311  }
312 
313  int add_patch()
314  {
315  int patch_index = patch_.append_and_get_index(Patch());
316  return patch_index;
317  }
318 
319  void grow_patch(int patch_index, int t)
320  {
321  tri_patch_[t] = patch_index;
322  patch_[patch_index].add_tri(t);
323  }
324 
325  bool tri_is_assigned(int t) const
326  {
327  return tri_patch_[t] != NO_INDEX;
328  }
329 
330  const Patch &patch(int patch_index) const
331  {
332  return patch_[patch_index];
333  }
334 
335  Patch &patch(int patch_index)
336  {
337  return patch_[patch_index];
338  }
339 
340  int tot_patch() const
341  {
342  return patch_.size();
343  }
344 
345  IndexRange index_range() const
346  {
347  return IndexRange(patch_.size());
348  }
349 
350  const Patch *begin() const
351  {
352  return patch_.begin();
353  }
354 
355  const Patch *end() const
356  {
357  return patch_.end();
358  }
359 
360  Patch *begin()
361  {
362  return patch_.begin();
363  }
364 
365  Patch *end()
366  {
367  return patch_.end();
368  }
369 
370  void add_new_patch_patch_edge(int p1, int p2, Edge e)
371  {
372  pp_edge_.add_new(std::pair<int, int>(p1, p2), e);
373  pp_edge_.add_new(std::pair<int, int>(p2, p1), e);
374  }
375 
376  Edge patch_patch_edge(int p1, int p2)
377  {
378  return pp_edge_.lookup_default(std::pair<int, int>(p1, p2), Edge());
379  }
380 
381  const Map<std::pair<int, int>, Edge> &patch_patch_edge_map()
382  {
383  return pp_edge_;
384  }
385 };
386 
387 static bool apply_bool_op(BoolOpType bool_optype, const Array<int> &winding);
388 
394 class Cell {
395  Set<int> patches_;
396  Array<int> winding_;
397  int merged_to_{NO_INDEX};
398  bool winding_assigned_{false};
399  /* in_output_volume_ will be true when this cell should be in the output volume. */
400  bool in_output_volume_{false};
401  /* zero_volume_ will be true when this is a zero-volume cell (inside a stack of identical
402  * triangles). */
403  bool zero_volume_{false};
404 
405  public:
406  Cell() = default;
407 
408  void add_patch(int p)
409  {
410  patches_.add(p);
411  zero_volume_ = false; /* If it was true before, it no longer is. */
412  }
413 
414  const Set<int> &patches() const
415  {
416  return patches_;
417  }
418 
420  int patch_other(int p) const
421  {
422  if (patches_.size() != 2) {
423  return NO_INDEX;
424  }
425  for (int pother : patches_) {
426  if (pother != p) {
427  return pother;
428  }
429  }
430  return NO_INDEX;
431  }
432 
433  Span<int> winding() const
434  {
435  return Span<int>(winding_);
436  }
437 
438  void init_winding(int winding_len)
439  {
440  winding_ = Array<int>(winding_len);
441  }
442 
443  void seed_ambient_winding()
444  {
445  winding_.fill(0);
446  winding_assigned_ = true;
447  }
448 
449  void set_winding_and_in_output_volume(const Cell &from_cell,
450  int shape,
451  int delta,
452  BoolOpType bool_optype)
453  {
454  std::copy(from_cell.winding().begin(), from_cell.winding().end(), winding_.begin());
455  if (shape >= 0) {
456  winding_[shape] += delta;
457  }
458  winding_assigned_ = true;
459  in_output_volume_ = apply_bool_op(bool_optype, winding_);
460  }
461 
462  bool in_output_volume() const
463  {
464  return in_output_volume_;
465  }
466 
467  bool winding_assigned() const
468  {
469  return winding_assigned_;
470  }
471 
472  bool zero_volume() const
473  {
474  return zero_volume_;
475  }
476 
477  int merged_to() const
478  {
479  return merged_to_;
480  }
481 
482  void set_merged_to(int c)
483  {
484  merged_to_ = c;
485  }
486 
491  void check_for_zero_volume(const PatchesInfo &pinfo, const IMesh &mesh);
492 };
493 
494 static std::ostream &operator<<(std::ostream &os, const Cell &cell)
495 {
496  os << "Cell patches";
497  for (int p : cell.patches()) {
498  std::cout << " " << p;
499  }
500  if (cell.winding().size() > 0) {
501  os << " winding=" << cell.winding();
502  os << " in_output_volume=" << cell.in_output_volume();
503  }
504  os << " zv=" << cell.zero_volume();
505  std::cout << "\n";
506  return os;
507 }
508 
509 static bool tris_have_same_verts(const IMesh &mesh, int t1, int t2)
510 {
511  const Face &tri1 = *mesh.face(t1);
512  const Face &tri2 = *mesh.face(t2);
513  BLI_assert(tri1.size() == 3 && tri2.size() == 3);
514  if (tri1.vert[0] == tri2.vert[0]) {
515  return ((tri1.vert[1] == tri2.vert[1] && tri1.vert[2] == tri2.vert[2]) ||
516  (tri1.vert[1] == tri2.vert[2] && tri1.vert[2] == tri2.vert[1]));
517  }
518  if (tri1.vert[0] == tri2.vert[1]) {
519  return ((tri1.vert[1] == tri2.vert[0] && tri1.vert[2] == tri2.vert[2]) ||
520  (tri1.vert[1] == tri2.vert[2] && tri1.vert[2] == tri2.vert[0]));
521  }
522  if (tri1.vert[0] == tri2.vert[2]) {
523  return ((tri1.vert[1] == tri2.vert[0] && tri1.vert[2] == tri2.vert[1]) ||
524  (tri1.vert[1] == tri2.vert[1] && tri1.vert[2] == tri2.vert[0]));
525  }
526  return false;
527 }
528 
534 void Cell::check_for_zero_volume(const PatchesInfo &pinfo, const IMesh &mesh)
535 {
536  if (patches_.size() == 2) {
537  int p1_index = NO_INDEX;
538  int p2_index = NO_INDEX;
539  for (int p : patches_) {
540  if (p1_index == NO_INDEX) {
541  p1_index = p;
542  }
543  else {
544  p2_index = p;
545  }
546  }
547  BLI_assert(p1_index != NO_INDEX && p2_index != NO_INDEX);
548  const Patch &p1 = pinfo.patch(p1_index);
549  const Patch &p2 = pinfo.patch(p2_index);
550  if (p1.tot_tri() == 1 && p2.tot_tri() == 1) {
551  if (tris_have_same_verts(mesh, p1.tri(0), p2.tri(0))) {
552  zero_volume_ = true;
553  }
554  }
555  }
556 }
557 
558 /* Information about all the Cells. */
559 class CellsInfo {
560  Vector<Cell> cell_;
561 
562  public:
563  CellsInfo() = default;
564 
565  int add_cell()
566  {
567  int index = cell_.append_and_get_index(Cell());
568  return index;
569  }
570 
571  Cell &cell(int c)
572  {
573  return cell_[c];
574  }
575 
576  const Cell &cell(int c) const
577  {
578  return cell_[c];
579  }
580 
581  int tot_cell() const
582  {
583  return cell_.size();
584  }
585 
586  IndexRange index_range() const
587  {
588  return cell_.index_range();
589  }
590 
591  const Cell *begin() const
592  {
593  return cell_.begin();
594  }
595 
596  const Cell *end() const
597  {
598  return cell_.end();
599  }
600 
601  Cell *begin()
602  {
603  return cell_.begin();
604  }
605 
606  Cell *end()
607  {
608  return cell_.end();
609  }
610 
611  void init_windings(int winding_len)
612  {
613  for (Cell &cell : cell_) {
614  cell.init_winding(winding_len);
615  }
616  }
617 };
618 
622 static void write_obj_cell_patch(const IMesh &m,
623  const CellsInfo &cinfo,
624  const PatchesInfo &pinfo,
625  bool cells_only,
626  const std::string &name)
627 {
628  /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library.
629  * This is just for developer debugging anyway,
630  * and should never be called in production Blender. */
631 # ifdef _WIN_32
632  const char *objdir = BLI_getenv("HOME");
633 # else
634  const char *objdir = "/tmp/";
635 # endif
636 
637  std::string fname = std::string(objdir) + name + std::string("_cellpatch.obj");
638  std::ofstream f;
639  f.open(fname);
640  if (!f) {
641  std::cout << "Could not open file " << fname << "\n";
642  return;
643  }
644 
645  /* Copy IMesh so can populate verts. */
646  IMesh mm = m;
647  mm.populate_vert();
648  f << "o cellpatch\n";
649  for (const Vert *v : mm.vertices()) {
650  const double3 dv = v->co;
651  f << "v " << dv[0] << " " << dv[1] << " " << dv[2] << "\n";
652  }
653  if (!cells_only) {
654  for (int p : pinfo.index_range()) {
655  f << "g patch" << p << "\n";
656  const Patch &patch = pinfo.patch(p);
657  for (int t : patch.tris()) {
658  const Face &tri = *mm.face(t);
659  f << "f ";
660  for (const Vert *v : tri) {
661  f << mm.lookup_vert(v) + 1 << " ";
662  }
663  f << "\n";
664  }
665  }
666  }
667  for (int c : cinfo.index_range()) {
668  f << "g cell" << c << "\n";
669  const Cell &cell = cinfo.cell(c);
670  for (int p : cell.patches()) {
671  const Patch &patch = pinfo.patch(p);
672  for (int t : patch.tris()) {
673  const Face &tri = *mm.face(t);
674  f << "f ";
675  for (const Vert *v : tri) {
676  f << mm.lookup_vert(v) + 1 << " ";
677  }
678  f << "\n";
679  }
680  }
681  }
682  f.close();
683 }
684 
685 static void merge_cells(int merge_to, int merge_from, CellsInfo &cinfo, PatchesInfo &pinfo)
686 {
687  if (merge_to == merge_from) {
688  return;
689  }
690  Cell &merge_from_cell = cinfo.cell(merge_from);
691  Cell &merge_to_cell = cinfo.cell(merge_to);
692  int final_merge_to = merge_to;
693  while (merge_to_cell.merged_to() != NO_INDEX) {
694  final_merge_to = merge_to_cell.merged_to();
695  merge_to_cell = cinfo.cell(final_merge_to);
696  }
697  for (int cell_p : merge_from_cell.patches()) {
698  merge_to_cell.add_patch(cell_p);
699  Patch &patch = pinfo.patch(cell_p);
700  if (patch.cell_above == merge_from) {
701  patch.cell_above = merge_to;
702  }
703  if (patch.cell_below == merge_from) {
704  patch.cell_below = merge_to;
705  }
706  }
707  merge_from_cell.set_merged_to(final_merge_to);
708 }
709 
713 static PatchesInfo find_patches(const IMesh &tm, const TriMeshTopology &tmtopo)
714 {
715  const int dbg_level = 0;
716  if (dbg_level > 0) {
717  std::cout << "\nFIND_PATCHES\n";
718  }
719  int ntri = tm.face_size();
720  PatchesInfo pinfo(ntri);
721  /* Algorithm: Grow patches across manifold edges as long as there are unassigned triangles. */
722  Stack<int> cur_patch_grow;
723 
724  /* Create an Array containing indices of adjacent faces. */
725  Array<std::array<int, 3>> t_others(tm.face_size());
726  threading::parallel_for(tm.face_index_range(), 2048, [&](IndexRange range) {
727  for (int t : range) {
728  const Face &tri = *tm.face(t);
729  for (int i = 0; i < 3; ++i) {
730  Edge e(tri[i], tri[(i + 1) % 3]);
731  t_others[t][i] = tmtopo.other_tri_if_manifold(e, t);
732  }
733  }
734  });
735  for (int t : tm.face_index_range()) {
736  if (pinfo.tri_patch(t) == -1) {
737  cur_patch_grow.push(t);
738  int cur_patch_index = pinfo.add_patch();
739  while (!cur_patch_grow.is_empty()) {
740  int tcand = cur_patch_grow.pop();
741  if (dbg_level > 1) {
742  std::cout << "pop tcand = " << tcand << "; assigned = " << pinfo.tri_is_assigned(tcand)
743  << "\n";
744  }
745  if (pinfo.tri_is_assigned(tcand)) {
746  continue;
747  }
748  if (dbg_level > 1) {
749  std::cout << "grow patch from seed tcand=" << tcand << "\n";
750  }
751  pinfo.grow_patch(cur_patch_index, tcand);
752  const Face &tri = *tm.face(tcand);
753  for (int i = 0; i < 3; ++i) {
754  Edge e(tri[i], tri[(i + 1) % 3]);
755  int t_other = t_others[tcand][i];
756  if (dbg_level > 1) {
757  std::cout << " edge " << e << " generates t_other=" << t_other << "\n";
758  }
759  if (t_other != NO_INDEX) {
760  if (!pinfo.tri_is_assigned(t_other)) {
761  if (dbg_level > 1) {
762  std::cout << " push t_other = " << t_other << "\n";
763  }
764  cur_patch_grow.push(t_other);
765  }
766  }
767  else {
768  /* e is non-manifold. Set any patch-patch incidences we can. */
769  if (dbg_level > 1) {
770  std::cout << " e non-manifold case\n";
771  }
772  const Vector<int> *etris = tmtopo.edge_tris(e);
773  if (etris != nullptr) {
774  for (int i : etris->index_range()) {
775  int t_other = (*etris)[i];
776  if (t_other != tcand && pinfo.tri_is_assigned(t_other)) {
777  int p_other = pinfo.tri_patch(t_other);
778  if (p_other == cur_patch_index) {
779  continue;
780  }
781  if (pinfo.patch_patch_edge(cur_patch_index, p_other).v0() == nullptr) {
782  pinfo.add_new_patch_patch_edge(cur_patch_index, p_other, e);
783  if (dbg_level > 1) {
784  std::cout << "added patch_patch_edge (" << cur_patch_index << "," << p_other
785  << ") = " << e << "\n";
786  }
787  }
788  }
789  }
790  }
791  }
792  }
793  }
794  }
795  }
796  if (dbg_level > 0) {
797  std::cout << "\nafter FIND_PATCHES: found " << pinfo.tot_patch() << " patches\n";
798  for (int p : pinfo.index_range()) {
799  std::cout << p << ": " << pinfo.patch(p) << "\n";
800  }
801  if (dbg_level > 1) {
802  std::cout << "\ntriangle map\n";
803  for (int t : tm.face_index_range()) {
804  std::cout << t << ": " << tm.face(t) << " patch " << pinfo.tri_patch(t) << "\n";
805  }
806  }
807  std::cout << "\npatch-patch incidences\n";
808  for (int p1 : pinfo.index_range()) {
809  for (int p2 : pinfo.index_range()) {
810  Edge e = pinfo.patch_patch_edge(p1, p2);
811  if (e.v0() != nullptr) {
812  std::cout << "p" << p1 << " and p" << p2 << " share edge " << e << "\n";
813  }
814  }
815  }
816  }
817  return pinfo;
818 }
819 
826 static const Vert *find_flap_vert(const Face &tri, const Edge e, bool *r_rev)
827 {
828  *r_rev = false;
829  const Vert *flapv;
830  if (tri[0] == e.v0()) {
831  if (tri[1] == e.v1()) {
832  *r_rev = false;
833  flapv = tri[2];
834  }
835  else {
836  if (tri[2] != e.v1()) {
837  return nullptr;
838  }
839  *r_rev = true;
840  flapv = tri[1];
841  }
842  }
843  else if (tri[1] == e.v0()) {
844  if (tri[2] == e.v1()) {
845  *r_rev = false;
846  flapv = tri[0];
847  }
848  else {
849  if (tri[0] != e.v1()) {
850  return nullptr;
851  }
852  *r_rev = true;
853  flapv = tri[2];
854  }
855  }
856  else {
857  if (tri[2] != e.v0()) {
858  return nullptr;
859  }
860  if (tri[0] == e.v1()) {
861  *r_rev = false;
862  flapv = tri[1];
863  }
864  else {
865  if (tri[1] != e.v1()) {
866  return nullptr;
867  }
868  *r_rev = true;
869  flapv = tri[0];
870  }
871  }
872  return flapv;
873 }
874 
889 static int sort_tris_class(const Face &tri, const Face &tri0, const Edge e)
890 {
891  const int dbg_level = 0;
892  if (dbg_level > 0) {
893  std::cout << "classify e = " << e << "\n";
894  }
895  mpq3 a0 = tri0[0]->co_exact;
896  mpq3 a1 = tri0[1]->co_exact;
897  mpq3 a2 = tri0[2]->co_exact;
898  bool rev;
899  bool rev0;
900  const Vert *flapv0 = find_flap_vert(tri0, e, &rev0);
901  const Vert *flapv = find_flap_vert(tri, e, &rev);
902  if (dbg_level > 0) {
903  std::cout << " t0 = " << tri0[0] << " " << tri0[1] << " " << tri0[2];
904  std::cout << " rev0 = " << rev0 << " flapv0 = " << flapv0 << "\n";
905  std::cout << " t = " << tri[0] << " " << tri[1] << " " << tri[2];
906  std::cout << " rev = " << rev << " flapv = " << flapv << "\n";
907  }
908  BLI_assert(flapv != nullptr && flapv0 != nullptr);
909  const mpq3 flap = flapv->co_exact;
910  /* orient will be positive if flap is below oriented plane of a0,a1,a2. */
911  int orient = orient3d(a0, a1, a2, flap);
912  int ans;
913  if (orient > 0) {
914  ans = rev0 ? 4 : 3;
915  }
916  else if (orient < 0) {
917  ans = rev0 ? 3 : 4;
918  }
919  else {
920  ans = flapv == flapv0 ? 1 : 2;
921  }
922  if (dbg_level > 0) {
923  std::cout << " orient = " << orient << " ans = " << ans << "\n";
924  }
925  return ans;
926 }
927 
928 constexpr int EXTRA_TRI_INDEX = INT_MAX;
929 
936 static void sort_by_signed_triangle_index(Vector<int> &g,
937  const Edge e,
938  const IMesh &tm,
939  const Face *extra_tri)
940 {
941  Array<int> signed_g(g.size());
942  for (int i : g.index_range()) {
943  const Face &tri = g[i] == EXTRA_TRI_INDEX ? *extra_tri : *tm.face(g[i]);
944  bool rev;
945  find_flap_vert(tri, e, &rev);
946  signed_g[i] = rev ? -g[i] : g[i];
947  }
948  std::sort(signed_g.begin(), signed_g.end());
949 
950  for (int i : g.index_range()) {
951  g[i] = abs(signed_g[i]);
952  }
953 }
954 
969 static Array<int> sort_tris_around_edge(
970  const IMesh &tm, const Edge e, const Span<int> tris, const int t0, const Face *extra_tri)
971 {
972  /* Divide and conquer, quick-sort-like sort.
973  * Pick a triangle t0, then partition into groups:
974  * (1) co-planar with t0 and on same side of e
975  * (2) co-planar with t0 and on opposite side of e
976  * (3) below plane of t0
977  * (4) above plane of t0
978  * Each group is sorted and then the sorts are merged to give the answer.
979  * We don't expect the input array to be very large - should typically
980  * be only 3 or 4 - so OK to make copies of arrays instead of swapping
981  * around in a single array. */
982  const int dbg_level = 0;
983  if (tris.size() == 0) {
984  return Array<int>();
985  }
986  if (dbg_level > 0) {
987  if (t0 == tris[0]) {
988  std::cout << "\n";
989  }
990  std::cout << "sort_tris_around_edge " << e << "\n";
991  std::cout << "tris = " << tris << "\n";
992  std::cout << "t0 = " << t0 << "\n";
993  }
994  Vector<int> g1{tris[0]};
995  Vector<int> g2;
996  Vector<int> g3;
997  Vector<int> g4;
998  std::array<Vector<int> *, 4> groups = {&g1, &g2, &g3, &g4};
999  const Face &triref = *tm.face(tris[0]);
1000  for (int i : tris.index_range()) {
1001  if (i == 0) {
1002  continue;
1003  }
1004  int t = tris[i];
1005  BLI_assert(t < tm.face_size() || (t == EXTRA_TRI_INDEX && extra_tri != nullptr));
1006  const Face &tri = (t == EXTRA_TRI_INDEX) ? *extra_tri : *tm.face(t);
1007  if (dbg_level > 2) {
1008  std::cout << "classifying tri " << t << " with respect to " << tris[0] << "\n";
1009  }
1010  int group_num = sort_tris_class(tri, triref, e);
1011  if (dbg_level > 2) {
1012  std::cout << " classify result : " << group_num << "\n";
1013  }
1014  groups[group_num - 1]->append(t);
1015  }
1016  if (dbg_level > 1) {
1017  std::cout << "g1 = " << g1 << "\n";
1018  std::cout << "g2 = " << g2 << "\n";
1019  std::cout << "g3 = " << g3 << "\n";
1020  std::cout << "g4 = " << g4 << "\n";
1021  }
1022  if (g1.size() > 1) {
1023  sort_by_signed_triangle_index(g1, e, tm, extra_tri);
1024  if (dbg_level > 1) {
1025  std::cout << "g1 sorted: " << g1 << "\n";
1026  }
1027  }
1028  if (g2.size() > 1) {
1029  sort_by_signed_triangle_index(g2, e, tm, extra_tri);
1030  if (dbg_level > 1) {
1031  std::cout << "g2 sorted: " << g2 << "\n";
1032  }
1033  }
1034  if (g3.size() > 1) {
1035  Array<int> g3sorted = sort_tris_around_edge(tm, e, g3, t0, extra_tri);
1036  std::copy(g3sorted.begin(), g3sorted.end(), g3.begin());
1037  if (dbg_level > 1) {
1038  std::cout << "g3 sorted: " << g3 << "\n";
1039  }
1040  }
1041  if (g4.size() > 1) {
1042  Array<int> g4sorted = sort_tris_around_edge(tm, e, g4, t0, extra_tri);
1043  std::copy(g4sorted.begin(), g4sorted.end(), g4.begin());
1044  if (dbg_level > 1) {
1045  std::cout << "g4 sorted: " << g4 << "\n";
1046  }
1047  }
1048  int group_tot_size = g1.size() + g2.size() + g3.size() + g4.size();
1049  Array<int> ans(group_tot_size);
1050  int *p = ans.begin();
1051  if (tris[0] == t0) {
1052  p = std::copy(g1.begin(), g1.end(), p);
1053  p = std::copy(g4.begin(), g4.end(), p);
1054  p = std::copy(g2.begin(), g2.end(), p);
1055  std::copy(g3.begin(), g3.end(), p);
1056  }
1057  else {
1058  p = std::copy(g3.begin(), g3.end(), p);
1059  p = std::copy(g1.begin(), g1.end(), p);
1060  p = std::copy(g4.begin(), g4.end(), p);
1061  std::copy(g2.begin(), g2.end(), p);
1062  }
1063  if (dbg_level > 0) {
1064  std::cout << "sorted tris = " << ans << "\n";
1065  }
1066  return ans;
1067 }
1068 
1075 static void find_cells_from_edge(const IMesh &tm,
1076  const TriMeshTopology &tmtopo,
1077  PatchesInfo &pinfo,
1078  CellsInfo &cinfo,
1079  const Edge e)
1080 {
1081  const int dbg_level = 0;
1082  if (dbg_level > 0) {
1083  std::cout << "FIND_CELLS_FROM_EDGE " << e << "\n";
1084  }
1085  const Vector<int> *edge_tris = tmtopo.edge_tris(e);
1086  BLI_assert(edge_tris != nullptr);
1087  Array<int> sorted_tris = sort_tris_around_edge(
1088  tm, e, Span<int>(*edge_tris), (*edge_tris)[0], nullptr);
1089 
1090  int n_edge_tris = edge_tris->size();
1091  Array<int> edge_patches(n_edge_tris);
1092  for (int i = 0; i < n_edge_tris; ++i) {
1093  edge_patches[i] = pinfo.tri_patch(sorted_tris[i]);
1094  if (dbg_level > 1) {
1095  std::cout << "edge_patches[" << i << "] = " << edge_patches[i] << "\n";
1096  }
1097  }
1098  for (int i = 0; i < n_edge_tris; ++i) {
1099  int inext = (i + 1) % n_edge_tris;
1100  int r_index = edge_patches[i];
1101  int rnext_index = edge_patches[inext];
1102  Patch &r = pinfo.patch(r_index);
1103  Patch &rnext = pinfo.patch(rnext_index);
1104  bool r_flipped;
1105  bool rnext_flipped;
1106  find_flap_vert(*tm.face(sorted_tris[i]), e, &r_flipped);
1107  find_flap_vert(*tm.face(sorted_tris[inext]), e, &rnext_flipped);
1108  int *r_follow_cell = r_flipped ? &r.cell_below : &r.cell_above;
1109  int *rnext_prev_cell = rnext_flipped ? &rnext.cell_above : &rnext.cell_below;
1110  if (dbg_level > 0) {
1111  std::cout << "process patch pair " << r_index << " " << rnext_index << "\n";
1112  std::cout << " r_flipped = " << r_flipped << " rnext_flipped = " << rnext_flipped << "\n";
1113  std::cout << " r_follow_cell (" << (r_flipped ? "below" : "above")
1114  << ") = " << *r_follow_cell << "\n";
1115  std::cout << " rnext_prev_cell (" << (rnext_flipped ? "above" : "below")
1116  << ") = " << *rnext_prev_cell << "\n";
1117  }
1118  if (*r_follow_cell == NO_INDEX && *rnext_prev_cell == NO_INDEX) {
1119  /* Neither is assigned: make a new cell. */
1120  int c = cinfo.add_cell();
1121  *r_follow_cell = c;
1122  *rnext_prev_cell = c;
1123  Cell &cell = cinfo.cell(c);
1124  cell.add_patch(r_index);
1125  cell.add_patch(rnext_index);
1126  cell.check_for_zero_volume(pinfo, tm);
1127  if (dbg_level > 0) {
1128  std::cout << " made new cell " << c << "\n";
1129  std::cout << " p" << r_index << "." << (r_flipped ? "cell_below" : "cell_above") << " = c"
1130  << c << "\n";
1131  std::cout << " p" << rnext_index << "." << (rnext_flipped ? "cell_above" : "cell_below")
1132  << " = c" << c << "\n";
1133  }
1134  }
1135  else if (*r_follow_cell != NO_INDEX && *rnext_prev_cell == NO_INDEX) {
1136  int c = *r_follow_cell;
1137  *rnext_prev_cell = c;
1138  Cell &cell = cinfo.cell(c);
1139  cell.add_patch(rnext_index);
1140  cell.check_for_zero_volume(pinfo, tm);
1141  if (dbg_level > 0) {
1142  std::cout << " reuse r_follow: p" << rnext_index << "."
1143  << (rnext_flipped ? "cell_above" : "cell_below") << " = c" << c << "\n";
1144  }
1145  }
1146  else if (*r_follow_cell == NO_INDEX && *rnext_prev_cell != NO_INDEX) {
1147  int c = *rnext_prev_cell;
1148  *r_follow_cell = c;
1149  Cell &cell = cinfo.cell(c);
1150  cell.add_patch(r_index);
1151  cell.check_for_zero_volume(pinfo, tm);
1152  if (dbg_level > 0) {
1153  std::cout << " reuse rnext prev: rprev_p" << r_index << "."
1154  << (r_flipped ? "cell_below" : "cell_above") << " = c" << c << "\n";
1155  }
1156  }
1157  else {
1158  if (*r_follow_cell != *rnext_prev_cell) {
1159  int follow_cell_num_patches = cinfo.cell(*r_follow_cell).patches().size();
1160  int prev_cell_num_patches = cinfo.cell(*rnext_prev_cell).patches().size();
1161  if (follow_cell_num_patches >= prev_cell_num_patches) {
1162  if (dbg_level > 0) {
1163  std::cout << " merge cell " << *rnext_prev_cell << " into cell " << *r_follow_cell
1164  << "\n";
1165  }
1166  merge_cells(*r_follow_cell, *rnext_prev_cell, cinfo, pinfo);
1167  }
1168  }
1169  else {
1170  if (dbg_level > 0) {
1171  std::cout << " merge cell " << *r_follow_cell << " into cell " << *rnext_prev_cell
1172  << "\n";
1173  }
1174  merge_cells(*rnext_prev_cell, *r_follow_cell, cinfo, pinfo);
1175  }
1176  }
1177  }
1178 }
1179 
1184 static CellsInfo find_cells(const IMesh &tm, const TriMeshTopology &tmtopo, PatchesInfo &pinfo)
1185 {
1186  const int dbg_level = 0;
1187  if (dbg_level > 0) {
1188  std::cout << "\nFIND_CELLS\n";
1189  }
1190  CellsInfo cinfo;
1191  /* For each unique edge shared between patch pairs, process it. */
1192  Set<Edge> processed_edges;
1193  for (const auto item : pinfo.patch_patch_edge_map().items()) {
1194  int p = item.key.first;
1195  int q = item.key.second;
1196  if (p < q) {
1197  const Edge &e = item.value;
1198  if (!processed_edges.contains(e)) {
1199  processed_edges.add_new(e);
1200  find_cells_from_edge(tm, tmtopo, pinfo, cinfo, e);
1201  }
1202  }
1203  }
1204  /* Some patches may have no cells at this point. These are either:
1205  * (a) a closed manifold patch only incident on itself (sphere, torus, klein bottle, etc.).
1206  * (b) an open manifold patch only incident on itself (has non-manifold boundaries).
1207  * Make above and below cells for these patches. This will create a disconnected patch-cell
1208  * bipartite graph, which will have to be fixed later. */
1209  for (int p : pinfo.index_range()) {
1210  Patch &patch = pinfo.patch(p);
1211  if (patch.cell_above == NO_INDEX) {
1212  int c = cinfo.add_cell();
1213  patch.cell_above = c;
1214  Cell &cell = cinfo.cell(c);
1215  cell.add_patch(p);
1216  }
1217  if (patch.cell_below == NO_INDEX) {
1218  int c = cinfo.add_cell();
1219  patch.cell_below = c;
1220  Cell &cell = cinfo.cell(c);
1221  cell.add_patch(p);
1222  }
1223  }
1224  if (dbg_level > 0) {
1225  std::cout << "\nFIND_CELLS found " << cinfo.tot_cell() << " cells\nCells\n";
1226  for (int i : cinfo.index_range()) {
1227  std::cout << i << ": " << cinfo.cell(i) << "\n";
1228  }
1229  std::cout << "Patches\n";
1230  for (int i : pinfo.index_range()) {
1231  std::cout << i << ": " << pinfo.patch(i) << "\n";
1232  }
1233  if (dbg_level > 1) {
1234  write_obj_cell_patch(tm, cinfo, pinfo, false, "postfindcells");
1235  }
1236  }
1237  return cinfo;
1238 }
1239 
1245 static Vector<Vector<int>> find_patch_components(const CellsInfo &cinfo, PatchesInfo &pinfo)
1246 {
1247  constexpr int dbg_level = 0;
1248  if (dbg_level > 0) {
1249  std::cout << "FIND_PATCH_COMPONENTS\n";
1250  }
1251  if (pinfo.tot_patch() == 0) {
1252  return Vector<Vector<int>>();
1253  }
1254  int current_component = 0;
1255  Array<bool> cell_processed(cinfo.tot_cell(), false);
1256  Stack<int> stack; /* Patch indices to visit. */
1257  Vector<Vector<int>> ans;
1258  for (int pstart : pinfo.index_range()) {
1259  Patch &patch_pstart = pinfo.patch(pstart);
1260  if (patch_pstart.component != NO_INDEX) {
1261  continue;
1262  }
1263  ans.append(Vector<int>());
1264  ans[current_component].append(pstart);
1265  stack.push(pstart);
1266  patch_pstart.component = current_component;
1267  while (!stack.is_empty()) {
1268  int p = stack.pop();
1269  Patch &patch = pinfo.patch(p);
1270  BLI_assert(patch.component == current_component);
1271  for (int c : {patch.cell_above, patch.cell_below}) {
1272  if (cell_processed[c]) {
1273  continue;
1274  }
1275  cell_processed[c] = true;
1276  for (int pn : cinfo.cell(c).patches()) {
1277  Patch &patch_neighbor = pinfo.patch(pn);
1278  if (patch_neighbor.component == NO_INDEX) {
1279  patch_neighbor.component = current_component;
1280  stack.push(pn);
1281  ans[current_component].append(pn);
1282  }
1283  }
1284  }
1285  }
1286  ++current_component;
1287  }
1288  if (dbg_level > 0) {
1289  std::cout << "found " << ans.size() << " components\n";
1290  for (int comp : ans.index_range()) {
1291  std::cout << comp << ": " << ans[comp] << "\n";
1292  }
1293  }
1294  return ans;
1295 }
1296 
1301 static bool patch_cell_graph_ok(const CellsInfo &cinfo, const PatchesInfo &pinfo)
1302 {
1303  for (int c : cinfo.index_range()) {
1304  const Cell &cell = cinfo.cell(c);
1305  if (cell.merged_to() != NO_INDEX) {
1306  continue;
1307  }
1308  if (cell.patches().size() == 0) {
1309  std::cout << "Patch/Cell graph disconnected at Cell " << c << " with no patches\n";
1310  return false;
1311  }
1312  for (int p : cell.patches()) {
1313  if (p >= pinfo.tot_patch()) {
1314  std::cout << "Patch/Cell graph has bad patch index at Cell " << c << "\n";
1315  return false;
1316  }
1317  }
1318  }
1319  for (int p : pinfo.index_range()) {
1320  const Patch &patch = pinfo.patch(p);
1321  if (patch.cell_above == NO_INDEX || patch.cell_below == NO_INDEX) {
1322  std::cout << "Patch/Cell graph disconnected at Patch " << p
1323  << " with one or two missing cells\n";
1324  return false;
1325  }
1326  if (patch.cell_above >= cinfo.tot_cell() || patch.cell_below >= cinfo.tot_cell()) {
1327  std::cout << "Patch/Cell graph has bad cell index at Patch " << p << "\n";
1328  return false;
1329  }
1330  }
1331  return true;
1332 }
1333 
1347 static bool is_pwn(const IMesh &tm, const TriMeshTopology &tmtopo)
1348 {
1349  constexpr int dbg_level = 0;
1350  std::atomic<bool> is_pwn = true;
1351  Vector<std::pair<Edge, Vector<int> *>> tris;
1352 
1353  for (auto item : tmtopo.edge_tri_map_items()) {
1354  tris.append(std::pair<Edge, Vector<int> *>(item.key, item.value));
1355  }
1356 
1357  threading::parallel_for(tris.index_range(), 2048, [&](IndexRange range) {
1358  if (!is_pwn.load()) {
1359  /* Early out if mesh is already determined to be non-pwn. */
1360  return;
1361  }
1362 
1363  for (int j : range) {
1364  const Edge &edge = tris[j].first;
1365  int tot_orient = 0;
1366  /* For each face t attached to edge, add +1 if the edge
1367  * is positively in t, and -1 if negatively in t. */
1368  for (int t : *tris[j].second) {
1369  const Face &face = *tm.face(t);
1370  BLI_assert(face.size() == 3);
1371  for (int i : face.index_range()) {
1372  if (face[i] == edge.v0()) {
1373  if (face[(i + 1) % 3] == edge.v1()) {
1374  ++tot_orient;
1375  }
1376  else {
1377  BLI_assert(face[(i + 3 - 1) % 3] == edge.v1());
1378  --tot_orient;
1379  }
1380  }
1381  }
1382  }
1383  if (tot_orient != 0) {
1384  if (dbg_level > 0) {
1385  std::cout << "edge causing non-pwn: " << edge << "\n";
1386  }
1387  is_pwn = false;
1388  break;
1389  }
1390  }
1391  });
1392  return is_pwn.load();
1393 }
1394 
1402 static int find_cell_for_point_near_edge(mpq3 p,
1403  const Edge &e,
1404  const IMesh &tm,
1405  const TriMeshTopology &tmtopo,
1406  const PatchesInfo &pinfo,
1407  IMeshArena *arena)
1408 {
1409  constexpr int dbg_level = 0;
1410  if (dbg_level > 0) {
1411  std::cout << "FIND_CELL_FOR_POINT_NEAR_EDGE, p=" << p << " e=" << e << "\n";
1412  }
1413  const Vector<int> *etris = tmtopo.edge_tris(e);
1414  const Vert *dummy_vert = arena->add_or_find_vert(p, NO_INDEX);
1415  const Face *dummy_tri = arena->add_face({e.v0(), e.v1(), dummy_vert},
1416  NO_INDEX,
1417  {NO_INDEX, NO_INDEX, NO_INDEX},
1418  {false, false, false});
1419  BLI_assert(etris != nullptr);
1420  Array<int> edge_tris(etris->size() + 1);
1421  std::copy(etris->begin(), etris->end(), edge_tris.begin());
1422  edge_tris[edge_tris.size() - 1] = EXTRA_TRI_INDEX;
1423  Array<int> sorted_tris = sort_tris_around_edge(tm, e, edge_tris, edge_tris[0], dummy_tri);
1424  if (dbg_level > 0) {
1425  std::cout << "sorted tris = " << sorted_tris << "\n";
1426  }
1427  int *p_sorted_dummy = std::find(sorted_tris.begin(), sorted_tris.end(), EXTRA_TRI_INDEX);
1428  BLI_assert(p_sorted_dummy != sorted_tris.end());
1429  int dummy_index = p_sorted_dummy - sorted_tris.begin();
1430  int prev_tri = (dummy_index == 0) ? sorted_tris[sorted_tris.size() - 1] :
1431  sorted_tris[dummy_index - 1];
1432  if (dbg_level > 0) {
1433  int next_tri = (dummy_index == sorted_tris.size() - 1) ? sorted_tris[0] :
1434  sorted_tris[dummy_index + 1];
1435  std::cout << "prev tri to dummy = " << prev_tri << "; next tri to dummy = " << next_tri
1436  << "\n";
1437  }
1438  const Patch &prev_patch = pinfo.patch(pinfo.tri_patch(prev_tri));
1439  if (dbg_level > 0) {
1440  std::cout << "prev_patch = " << prev_patch << "\n";
1441  }
1442  bool prev_flipped;
1443  find_flap_vert(*tm.face(prev_tri), e, &prev_flipped);
1444  int c = prev_flipped ? prev_patch.cell_below : prev_patch.cell_above;
1445  if (dbg_level > 0) {
1446  std::cout << "find_cell_for_point_near_edge returns " << c << "\n";
1447  }
1448  return c;
1449 }
1450 
1465 static int find_ambient_cell(const IMesh &tm,
1466  const Vector<int> *component_patches,
1467  const TriMeshTopology &tmtopo,
1468  const PatchesInfo &pinfo,
1469  IMeshArena *arena)
1470 {
1471  int dbg_level = 0;
1472  if (dbg_level > 0) {
1473  std::cout << "FIND_AMBIENT_CELL\n";
1474  }
1475  /* First find a vertex with the maximum x value. */
1476  /* Prefer not to populate the verts in the #IMesh just for this. */
1477  const Vert *v_extreme;
1478  auto max_x_vert = [](const Vert *a, const Vert *b) {
1479  return (a->co_exact.x > b->co_exact.x) ? a : b;
1480  };
1481  if (component_patches == nullptr) {
1482  v_extreme = threading::parallel_reduce(
1483  tm.face_index_range(),
1484  2048,
1485  (*tm.face(0))[0],
1486  [&](IndexRange range, const Vert *init) {
1487  const Vert *ans = init;
1488  for (int i : range) {
1489  const Face *f = tm.face(i);
1490  for (const Vert *v : *f) {
1491  if (v->co_exact.x > ans->co_exact.x) {
1492  ans = v;
1493  }
1494  }
1495  }
1496  return ans;
1497  },
1498  max_x_vert);
1499  }
1500  else {
1501  if (dbg_level > 0) {
1502  std::cout << "restrict to patches " << *component_patches << "\n";
1503  }
1504  int p0 = (*component_patches)[0];
1505  v_extreme = threading::parallel_reduce(
1506  component_patches->index_range(),
1507  2048,
1508  (*tm.face(pinfo.patch(p0).tri(0)))[0],
1509  [&](IndexRange range, const Vert *init) {
1510  const Vert *ans = init;
1511  for (int pi : range) {
1512  int p = (*component_patches)[pi];
1513  const Vert *tris_ans = threading::parallel_reduce(
1514  IndexRange(pinfo.patch(p).tot_tri()),
1515  2048,
1516  init,
1517  [&](IndexRange tris_range, const Vert *t_init) {
1518  const Vert *v_ans = t_init;
1519  for (int i : tris_range) {
1520  int t = pinfo.patch(p).tri(i);
1521  const Face *f = tm.face(t);
1522  for (const Vert *v : *f) {
1523  if (v->co_exact.x > v_ans->co_exact.x) {
1524  v_ans = v;
1525  }
1526  }
1527  }
1528  return v_ans;
1529  },
1530  max_x_vert);
1531  if (tris_ans->co_exact.x > ans->co_exact.x) {
1532  ans = tris_ans;
1533  }
1534  }
1535  return ans;
1536  },
1537  max_x_vert);
1538  }
1539  if (dbg_level > 0) {
1540  std::cout << "v_extreme = " << v_extreme << "\n";
1541  }
1542  /* Find edge attached to v_extreme with max absolute slope
1543  * when projected onto the XY plane. That edge is guaranteed to
1544  * be on the convex hull of the mesh. */
1545  const Vector<Edge> &edges = tmtopo.vert_edges(v_extreme);
1546  const mpq_class &extreme_x = v_extreme->co_exact.x;
1547  const mpq_class &extreme_y = v_extreme->co_exact.y;
1548  Edge ehull;
1549  mpq_class max_abs_slope = -1;
1550  for (Edge e : edges) {
1551  const Vert *v_other = (e.v0() == v_extreme) ? e.v1() : e.v0();
1552  const mpq3 &co_other = v_other->co_exact;
1553  mpq_class delta_x = co_other.x - extreme_x;
1554  if (delta_x == 0) {
1555  /* Vertical slope. */
1556  ehull = e;
1557  break;
1558  }
1559  mpq_class abs_slope = abs((co_other.y - extreme_y) / delta_x);
1560  if (abs_slope > max_abs_slope) {
1561  ehull = e;
1562  max_abs_slope = abs_slope;
1563  }
1564  }
1565  if (dbg_level > 0) {
1566  std::cout << "ehull = " << ehull << " slope = " << max_abs_slope << "\n";
1567  }
1568  /* Sort triangles around ehull, including a dummy triangle that include a known point in
1569  * ambient cell. */
1570  mpq3 p_in_ambient = v_extreme->co_exact;
1571  p_in_ambient.x += 1;
1572  int c_ambient = find_cell_for_point_near_edge(p_in_ambient, ehull, tm, tmtopo, pinfo, arena);
1573  if (dbg_level > 0) {
1574  std::cout << "FIND_AMBIENT_CELL returns " << c_ambient << "\n";
1575  }
1576  return c_ambient;
1577 }
1578 
1588 static Edge find_good_sorting_edge(const Vert *testp,
1589  const Vert *closestp,
1590  const TriMeshTopology &tmtopo)
1591 {
1592  constexpr int dbg_level = 0;
1593  if (dbg_level > 0) {
1594  std::cout << "FIND_GOOD_SORTING_EDGE testp = " << testp << ", closestp = " << closestp << "\n";
1595  }
1596  /* We want to project the edges incident to closestp onto a plane
1597  * whose ordinate direction will be regarded as going from closestp to testp,
1598  * and whose abscissa direction is some perpendicular to that.
1599  * A perpendicular direction can be found by swapping two coordinates
1600  * and negating one, and zeroing out the third, being careful that one
1601  * of the swapped vertices is non-zero. */
1602  const mpq3 &co_closest = closestp->co_exact;
1603  const mpq3 &co_test = testp->co_exact;
1604  BLI_assert(co_test != co_closest);
1605  mpq3 abscissa = co_test - co_closest;
1606  /* Find a non-zero-component axis of abscissa. */
1607  int axis;
1608  for (axis = 0; axis < 3; ++axis) {
1609  if (abscissa[axis] != 0) {
1610  break;
1611  }
1612  }
1613  BLI_assert(axis < 3);
1614  int axis_next = (axis + 1) % 3;
1615  int axis_next_next = (axis_next + 1) % 3;
1616  mpq3 ordinate;
1617  ordinate[axis] = abscissa[axis_next];
1618  ordinate[axis_next] = -abscissa[axis];
1619  ordinate[axis_next_next] = 0;
1620  /* By construction, dot(abscissa, ordinate) == 0, so they are perpendicular. */
1621  mpq3 normal = math::cross(abscissa, ordinate);
1622  if (dbg_level > 0) {
1623  std::cout << "abscissa = " << abscissa << "\n";
1624  std::cout << "ordinate = " << ordinate << "\n";
1625  std::cout << "normal = " << normal << "\n";
1626  }
1627  mpq_class nlen2 = math::length_squared(normal);
1628  mpq_class max_abs_slope = -1;
1629  Edge esort;
1630  const Vector<Edge> &edges = tmtopo.vert_edges(closestp);
1631  for (Edge e : edges) {
1632  const Vert *v_other = (e.v0() == closestp) ? e.v1() : e.v0();
1633  const mpq3 &co_other = v_other->co_exact;
1634  mpq3 evec = co_other - co_closest;
1635  /* Get projection of evec onto plane of abscissa and ordinate. */
1636  mpq3 proj_evec = evec - (math::dot(evec, normal) / nlen2) * normal;
1637  /* The projection calculations along the abscissa and ordinate should
1638  * be scaled by 1/abscissa and 1/ordinate respectively,
1639  * but we can skip: it won't affect which `evec` has the maximum slope. */
1640  mpq_class evec_a = math::dot(proj_evec, abscissa);
1641  mpq_class evec_o = math::dot(proj_evec, ordinate);
1642  if (dbg_level > 0) {
1643  std::cout << "e = " << e << "\n";
1644  std::cout << "v_other = " << v_other << "\n";
1645  std::cout << "evec = " << evec << ", proj_evec = " << proj_evec << "\n";
1646  std::cout << "evec_a = " << evec_a << ", evec_o = " << evec_o << "\n";
1647  }
1648  if (evec_a == 0) {
1649  /* evec is perpendicular to abscissa. */
1650  esort = e;
1651  if (dbg_level > 0) {
1652  std::cout << "perpendicular esort is " << esort << "\n";
1653  }
1654  break;
1655  }
1656  mpq_class abs_slope = abs(evec_o / evec_a);
1657  if (abs_slope > max_abs_slope) {
1658  esort = e;
1659  max_abs_slope = abs_slope;
1660  if (dbg_level > 0) {
1661  std::cout << "with abs_slope " << abs_slope << " new esort is " << esort << "\n";
1662  }
1663  }
1664  }
1665  return esort;
1666 }
1667 
1680 static int find_containing_cell(const Vert *v,
1681  int t,
1682  int close_edge,
1683  int close_vert,
1684  const PatchesInfo &pinfo,
1685  const IMesh &tm,
1686  const TriMeshTopology &tmtopo,
1687  IMeshArena *arena)
1688 {
1689  constexpr int dbg_level = 0;
1690  if (dbg_level > 0) {
1691  std::cout << "FIND_CONTAINING_CELL v=" << v << ", t=" << t << "\n";
1692  }
1693  const Face &tri = *tm.face(t);
1694  Edge etest;
1695  if (close_edge == -1 && close_vert == -1) {
1696  /* Choose any edge if closest point is inside the triangle. */
1697  close_edge = 0;
1698  }
1699  if (close_edge != -1) {
1700  const Vert *v0 = tri[close_edge];
1701  const Vert *v1 = tri[(close_edge + 1) % 3];
1702  const Vector<Edge> &edges = tmtopo.vert_edges(v0);
1703  if (dbg_level > 0) {
1704  std::cout << "look for edge containing " << v0 << " and " << v1 << "\n";
1705  std::cout << " in edges: ";
1706  for (Edge e : edges) {
1707  std::cout << e << " ";
1708  }
1709  std::cout << "\n";
1710  }
1711  for (Edge e : edges) {
1712  if ((e.v0() == v0 && e.v1() == v1) || (e.v0() == v1 && e.v1() == v0)) {
1713  etest = e;
1714  break;
1715  }
1716  }
1717  }
1718  else {
1719  int cv = close_vert;
1720  const Vert *vert_cv = tri[cv];
1721  if (vert_cv == v) {
1722  /* Need to use another one to find sorting edge. */
1723  vert_cv = tri[(cv + 1) % 3];
1724  BLI_assert(vert_cv != v);
1725  }
1726  etest = find_good_sorting_edge(v, vert_cv, tmtopo);
1727  }
1728  BLI_assert(etest.v0() != nullptr);
1729  if (dbg_level > 0) {
1730  std::cout << "etest = " << etest << "\n";
1731  }
1732  int c = find_cell_for_point_near_edge(v->co_exact, etest, tm, tmtopo, pinfo, arena);
1733  if (dbg_level > 0) {
1734  std::cout << "find_containing_cell returns " << c << "\n";
1735  }
1736  return c;
1737 }
1738 
1752 static mpq_class closest_on_tri_to_point(const mpq3 &p,
1753  const mpq3 &a,
1754  const mpq3 &b,
1755  const mpq3 &c,
1756  mpq3 &ab,
1757  mpq3 &ac,
1758  mpq3 &ap,
1759  mpq3 &bp,
1760  mpq3 &cp,
1761  mpq3 &m,
1762  mpq3 &r,
1763  int *r_edge,
1764  int *r_vert)
1765 {
1766  constexpr int dbg_level = 0;
1767  if (dbg_level > 0) {
1768  std::cout << "CLOSEST_ON_TRI_TO_POINT p = " << p << "\n";
1769  std::cout << " a = " << a << ", b = " << b << ", c = " << c << "\n";
1770  }
1771  /* Check if p in vertex region outside a. */
1772  ab = b;
1773  ab -= a;
1774  ac = c;
1775  ac -= a;
1776  ap = p;
1777  ap -= a;
1778 
1779  mpq_class d1 = math::dot_with_buffer(ab, ap, m);
1780  mpq_class d2 = math::dot_with_buffer(ac, ap, m);
1781  if (d1 <= 0 && d2 <= 0) {
1782  /* Barycentric coordinates (1,0,0). */
1783  *r_edge = -1;
1784  *r_vert = 0;
1785  if (dbg_level > 0) {
1786  std::cout << " answer = a\n";
1787  }
1788  return math::distance_squared_with_buffer(p, a, m);
1789  }
1790  /* Check if p in vertex region outside b. */
1791  bp = p;
1792  bp -= b;
1793  mpq_class d3 = math::dot_with_buffer(ab, bp, m);
1794  mpq_class d4 = math::dot_with_buffer(ac, bp, m);
1795  if (d3 >= 0 && d4 <= d3) {
1796  /* Barycentric coordinates (0,1,0). */
1797  *r_edge = -1;
1798  *r_vert = 1;
1799  if (dbg_level > 0) {
1800  std::cout << " answer = b\n";
1801  }
1802  return math::distance_squared_with_buffer(p, b, m);
1803  }
1804  /* Check if p in region of ab. */
1805  mpq_class vc = d1 * d4 - d3 * d2;
1806  if (vc <= 0 && d1 >= 0 && d3 <= 0) {
1807  mpq_class v = d1 / (d1 - d3);
1808  /* Barycentric coordinates (1-v,v,0). */
1809  r = ab;
1810  r *= v;
1811  r += a;
1812  *r_vert = -1;
1813  *r_edge = 0;
1814  if (dbg_level > 0) {
1815  std::cout << " answer = on ab at " << r << "\n";
1816  }
1817  return math::distance_squared_with_buffer(p, r, m);
1818  }
1819  /* Check if p in vertex region outside c. */
1820  cp = p;
1821  cp -= c;
1822  mpq_class d5 = math::dot_with_buffer(ab, cp, m);
1823  mpq_class d6 = math::dot_with_buffer(ac, cp, m);
1824  if (d6 >= 0 && d5 <= d6) {
1825  /* Barycentric coordinates (0,0,1). */
1826  *r_edge = -1;
1827  *r_vert = 2;
1828  if (dbg_level > 0) {
1829  std::cout << " answer = c\n";
1830  }
1831  return math::distance_squared_with_buffer(p, c, m);
1832  }
1833  /* Check if p in edge region of ac. */
1834  mpq_class vb = d5 * d2 - d1 * d6;
1835  if (vb <= 0 && d2 >= 0 && d6 <= 0) {
1836  mpq_class w = d2 / (d2 - d6);
1837  /* Barycentric coordinates (1-w,0,w). */
1838  r = ac;
1839  r *= w;
1840  r += a;
1841  *r_vert = -1;
1842  *r_edge = 2;
1843  if (dbg_level > 0) {
1844  std::cout << " answer = on ac at " << r << "\n";
1845  }
1846  return math::distance_squared_with_buffer(p, r, m);
1847  }
1848  /* Check if p in edge region of bc. */
1849  mpq_class va = d3 * d6 - d5 * d4;
1850  if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
1851  mpq_class w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
1852  /* Barycentric coordinates (0,1-w,w). */
1853  r = c;
1854  r -= b;
1855  r *= w;
1856  r += b;
1857  *r_vert = -1;
1858  *r_edge = 1;
1859  if (dbg_level > 0) {
1860  std::cout << " answer = on bc at " << r << "\n";
1861  }
1862  return math::distance_squared_with_buffer(p, r, m);
1863  }
1864  /* p inside face region. Compute barycentric coordinates (u,v,w). */
1865  mpq_class denom = 1 / (va + vb + vc);
1866  mpq_class v = vb * denom;
1867  mpq_class w = vc * denom;
1868  ac *= w;
1869  r = ab;
1870  r *= v;
1871  r += a;
1872  r += ac;
1873  *r_vert = -1;
1874  *r_edge = -1;
1875  if (dbg_level > 0) {
1876  std::cout << " answer = inside at " << r << "\n";
1877  }
1878  return math::distance_squared_with_buffer(p, r, m);
1879 }
1880 
1881 static float closest_on_tri_to_point_float_dist_squared(const float3 &p,
1882  const double3 &a,
1883  const double3 &b,
1884  const double3 &c)
1885 {
1886  float3 fa, fb, fc, closest;
1887  copy_v3fl_v3db(fa, a);
1888  copy_v3fl_v3db(fb, b);
1889  copy_v3fl_v3db(fc, c);
1890  closest_on_tri_to_point_v3(closest, p, fa, fb, fc);
1891  return len_squared_v3v3(p, closest);
1892 }
1893 
1894 struct ComponentContainer {
1895  int containing_component{NO_INDEX};
1896  int nearest_cell{NO_INDEX};
1897  mpq_class dist_to_cell;
1898 
1899  ComponentContainer(int cc, int cell, mpq_class d)
1900  : containing_component(cc), nearest_cell(cell), dist_to_cell(d)
1901  {
1902  }
1903 };
1904 
1911 static Vector<ComponentContainer> find_component_containers(int comp,
1912  const Vector<Vector<int>> &components,
1913  const Array<int> &ambient_cell,
1914  const IMesh &tm,
1915  const PatchesInfo &pinfo,
1916  const TriMeshTopology &tmtopo,
1917  Array<BoundingBox> &comp_bb,
1918  IMeshArena *arena)
1919 {
1920  constexpr int dbg_level = 0;
1921  if (dbg_level > 0) {
1922  std::cout << "FIND_COMPONENT_CONTAINERS for comp " << comp << "\n";
1923  }
1924  Vector<ComponentContainer> ans;
1925  int test_p = components[comp][0];
1926  int test_t = pinfo.patch(test_p).tri(0);
1927  const Vert *test_v = tm.face(test_t)[0].vert[0];
1928  if (dbg_level > 0) {
1929  std::cout << "test vertex in comp: " << test_v << "\n";
1930  }
1931  const double3 &test_v_d = test_v->co;
1932  float3 test_v_f(test_v_d[0], test_v_d[1], test_v_d[2]);
1933 
1934  mpq3 buf[7];
1935 
1936  for (int comp_other : components.index_range()) {
1937  if (comp == comp_other) {
1938  continue;
1939  }
1940  if (dbg_level > 0) {
1941  std::cout << "comp_other = " << comp_other << "\n";
1942  }
1943  if (!bbs_might_intersect(comp_bb[comp], comp_bb[comp_other])) {
1944  if (dbg_level > 0) {
1945  std::cout << "bounding boxes don't overlap\n";
1946  }
1947  continue;
1948  }
1949  int nearest_tri = NO_INDEX;
1950  int nearest_tri_close_vert = -1;
1951  int nearest_tri_close_edge = -1;
1952  mpq_class nearest_tri_dist_squared;
1953  float nearest_tri_dist_squared_float = FLT_MAX;
1954  for (int p : components[comp_other]) {
1955  const Patch &patch = pinfo.patch(p);
1956  for (int t : patch.tris()) {
1957  const Face &tri = *tm.face(t);
1958  if (dbg_level > 1) {
1959  std::cout << "tri " << t << " = " << &tri << "\n";
1960  }
1961  int close_vert;
1962  int close_edge;
1963  /* Try a cheap float test first. */
1964  float d2_f = closest_on_tri_to_point_float_dist_squared(
1965  test_v_f, tri[0]->co, tri[1]->co, tri[2]->co);
1966  if (d2_f - FLT_EPSILON > nearest_tri_dist_squared_float) {
1967  continue;
1968  }
1969  mpq_class d2 = closest_on_tri_to_point(test_v->co_exact,
1970  tri[0]->co_exact,
1971  tri[1]->co_exact,
1972  tri[2]->co_exact,
1973  buf[0],
1974  buf[1],
1975  buf[2],
1976  buf[3],
1977  buf[4],
1978  buf[5],
1979  buf[6],
1980  &close_edge,
1981  &close_vert);
1982  if (dbg_level > 1) {
1983  std::cout << " close_edge=" << close_edge << " close_vert=" << close_vert
1984  << " dsquared=" << d2.get_d() << "\n";
1985  }
1986  if (nearest_tri == NO_INDEX || d2 < nearest_tri_dist_squared) {
1987  nearest_tri = t;
1988  nearest_tri_close_edge = close_edge;
1989  nearest_tri_close_vert = close_vert;
1990  nearest_tri_dist_squared = d2;
1991  nearest_tri_dist_squared_float = d2_f;
1992  }
1993  }
1994  }
1995  if (dbg_level > 0) {
1996  std::cout << "closest tri to comp=" << comp << " in comp_other=" << comp_other << " is t"
1997  << nearest_tri << "\n";
1998  const Face *tn = tm.face(nearest_tri);
1999  std::cout << "tri = " << tn << "\n";
2000  std::cout << " (" << tn->vert[0]->co << "," << tn->vert[1]->co << "," << tn->vert[2]->co
2001  << ")\n";
2002  }
2003  int containing_cell = find_containing_cell(test_v,
2004  nearest_tri,
2005  nearest_tri_close_edge,
2006  nearest_tri_close_vert,
2007 
2008  pinfo,
2009  tm,
2010  tmtopo,
2011  arena);
2012  if (dbg_level > 0) {
2013  std::cout << "containing cell = " << containing_cell << "\n";
2014  }
2015  if (containing_cell != ambient_cell[comp_other]) {
2016  ans.append(ComponentContainer(comp_other, containing_cell, nearest_tri_dist_squared));
2017  }
2018  }
2019  return ans;
2020 }
2021 
2027 static void populate_comp_bbs(const Vector<Vector<int>> &components,
2028  const PatchesInfo &pinfo,
2029  const IMesh &im,
2030  Array<BoundingBox> &comp_bb)
2031 {
2032  const int comp_grainsize = 16;
2033  /* To get a good expansion epsilon, we need to find the maximum
2034  * absolute value of any coordinate. Do it first per component,
2035  * then get the overall max. */
2036  Array<double> max_abs(components.size(), 0.0);
2037  threading::parallel_for(components.index_range(), comp_grainsize, [&](IndexRange comp_range) {
2038  for (int c : comp_range) {
2039  BoundingBox &bb = comp_bb[c];
2040  double &maxa = max_abs[c];
2041  for (int p : components[c]) {
2042  const Patch &patch = pinfo.patch(p);
2043  for (int t : patch.tris()) {
2044  const Face &tri = *im.face(t);
2045  for (const Vert *v : tri) {
2046  bb.combine(v->co);
2047  for (int i = 0; i < 3; ++i) {
2048  maxa = max_dd(maxa, fabs(v->co[i]));
2049  }
2050  }
2051  }
2052  }
2053  }
2054  });
2055  double all_max_abs = 0.0;
2056  for (int c : components.index_range()) {
2057  all_max_abs = max_dd(all_max_abs, max_abs[c]);
2058  }
2059  constexpr float pad_factor = 10.0f;
2060  float pad = all_max_abs == 0.0 ? FLT_EPSILON : 2 * FLT_EPSILON * all_max_abs;
2061  pad *= pad_factor;
2062  for (int c : components.index_range()) {
2063  comp_bb[c].expand(pad);
2064  }
2065 }
2066 
2074 static void finish_patch_cell_graph(const IMesh &tm,
2075  CellsInfo &cinfo,
2076  PatchesInfo &pinfo,
2077  const TriMeshTopology &tmtopo,
2078  IMeshArena *arena)
2079 {
2080  constexpr int dbg_level = 0;
2081  if (dbg_level > 0) {
2082  std::cout << "FINISH_PATCH_CELL_GRAPH\n";
2083  }
2084  Vector<Vector<int>> components = find_patch_components(cinfo, pinfo);
2085  if (components.size() <= 1) {
2086  if (dbg_level > 0) {
2087  std::cout << "one component so finish_patch_cell_graph does no work\n";
2088  }
2089  return;
2090  }
2091  if (dbg_level > 0) {
2092  std::cout << "components:\n";
2093  for (int comp : components.index_range()) {
2094  std::cout << comp << ": " << components[comp] << "\n";
2095  }
2096  }
2097  Array<int> ambient_cell(components.size());
2098  for (int comp : components.index_range()) {
2099  ambient_cell[comp] = find_ambient_cell(tm, &components[comp], tmtopo, pinfo, arena);
2100  }
2101  if (dbg_level > 0) {
2102  std::cout << "ambient cells:\n";
2103  for (int comp : ambient_cell.index_range()) {
2104  std::cout << comp << ": " << ambient_cell[comp] << "\n";
2105  }
2106  }
2107  int tot_components = components.size();
2108  Array<Vector<ComponentContainer>> comp_cont(tot_components);
2109  if (tot_components > 1) {
2110  Array<BoundingBox> comp_bb(tot_components);
2111  populate_comp_bbs(components, pinfo, tm, comp_bb);
2112  for (int comp : components.index_range()) {
2113  comp_cont[comp] = find_component_containers(
2114  comp, components, ambient_cell, tm, pinfo, tmtopo, comp_bb, arena);
2115  }
2116  if (dbg_level > 0) {
2117  std::cout << "component containers:\n";
2118  for (int comp : comp_cont.index_range()) {
2119  std::cout << comp << ": ";
2120  for (const ComponentContainer &cc : comp_cont[comp]) {
2121  std::cout << "[containing_comp=" << cc.containing_component
2122  << ", nearest_cell=" << cc.nearest_cell << ", d2=" << cc.dist_to_cell << "] ";
2123  }
2124  std::cout << "\n";
2125  }
2126  }
2127  }
2128  if (dbg_level > 1) {
2129  write_obj_cell_patch(tm, cinfo, pinfo, false, "beforemerge");
2130  }
2131  /* For nested components, merge their ambient cell with the nearest containing cell. */
2132  Vector<int> outer_components;
2133  for (int comp : comp_cont.index_range()) {
2134  if (comp_cont[comp].size() == 0) {
2135  outer_components.append(comp);
2136  }
2137  else {
2138  ComponentContainer &closest = comp_cont[comp][0];
2139  for (int i = 1; i < comp_cont[comp].size(); ++i) {
2140  if (comp_cont[comp][i].dist_to_cell < closest.dist_to_cell) {
2141  closest = comp_cont[comp][i];
2142  }
2143  }
2144  int comp_ambient = ambient_cell[comp];
2145  int cont_cell = closest.nearest_cell;
2146  if (dbg_level > 0) {
2147  std::cout << "merge comp " << comp << "'s ambient cell=" << comp_ambient << " to cell "
2148  << cont_cell << "\n";
2149  }
2150  merge_cells(cont_cell, comp_ambient, cinfo, pinfo);
2151  }
2152  }
2153  /* For outer components (not nested in any other component), merge their ambient cells. */
2154  if (outer_components.size() > 1) {
2155  int merged_ambient = ambient_cell[outer_components[0]];
2156  for (int i = 1; i < outer_components.size(); ++i) {
2157  if (dbg_level > 0) {
2158  std::cout << "merge comp " << outer_components[i]
2159  << "'s ambient cell=" << ambient_cell[outer_components[i]] << " to cell "
2160  << merged_ambient << "\n";
2161  }
2162  merge_cells(merged_ambient, ambient_cell[outer_components[i]], cinfo, pinfo);
2163  }
2164  }
2165  if (dbg_level > 0) {
2166  std::cout << "after FINISH_PATCH_CELL_GRAPH\nCells\n";
2167  for (int i : cinfo.index_range()) {
2168  if (cinfo.cell(i).merged_to() == NO_INDEX) {
2169  std::cout << i << ": " << cinfo.cell(i) << "\n";
2170  }
2171  }
2172  std::cout << "Patches\n";
2173  for (int i : pinfo.index_range()) {
2174  std::cout << i << ": " << pinfo.patch(i) << "\n";
2175  }
2176  if (dbg_level > 1) {
2177  write_obj_cell_patch(tm, cinfo, pinfo, false, "finished");
2178  }
2179  }
2180 }
2181 
2195 static void propagate_windings_and_in_output_volume(PatchesInfo &pinfo,
2196  CellsInfo &cinfo,
2197  int c_ambient,
2198  BoolOpType op,
2199  int nshapes,
2200  std::function<int(int)> shape_fn)
2201 {
2202  int dbg_level = 0;
2203  if (dbg_level > 0) {
2204  std::cout << "PROPAGATE_WINDINGS, ambient cell = " << c_ambient << "\n";
2205  }
2206  Cell &cell_ambient = cinfo.cell(c_ambient);
2207  cell_ambient.seed_ambient_winding();
2208  /* Use a vector as a queue. It can't grow bigger than number of cells. */
2209  Vector<int> queue;
2210  queue.reserve(cinfo.tot_cell());
2211  int queue_head = 0;
2212  queue.append(c_ambient);
2213  while (queue_head < queue.size()) {
2214  int c = queue[queue_head++];
2215  if (dbg_level > 1) {
2216  std::cout << "process cell " << c << "\n";
2217  }
2218  Cell &cell = cinfo.cell(c);
2219  for (int p : cell.patches()) {
2220  Patch &patch = pinfo.patch(p);
2221  bool p_above_c = patch.cell_below == c;
2222  int c_neighbor = p_above_c ? patch.cell_above : patch.cell_below;
2223  if (dbg_level > 1) {
2224  std::cout << " patch " << p << " p_above_c = " << p_above_c << "\n";
2225  std::cout << " c_neighbor = " << c_neighbor << "\n";
2226  }
2227  Cell &cell_neighbor = cinfo.cell(c_neighbor);
2228  if (!cell_neighbor.winding_assigned()) {
2229  int winding_delta = p_above_c ? -1 : 1;
2230  int t = patch.tri(0);
2231  int shape = shape_fn(t);
2232  BLI_assert(shape < nshapes);
2233  UNUSED_VARS_NDEBUG(nshapes);
2234  if (dbg_level > 1) {
2235  std::cout << " representative tri " << t << ": in shape " << shape << "\n";
2236  }
2237  cell_neighbor.set_winding_and_in_output_volume(cell, shape, winding_delta, op);
2238  if (dbg_level > 1) {
2239  std::cout << " now cell_neighbor = " << cell_neighbor << "\n";
2240  }
2241  queue.append(c_neighbor);
2242  BLI_assert(queue.size() <= cinfo.tot_cell());
2243  }
2244  }
2245  }
2246  if (dbg_level > 0) {
2247  std::cout << "\nPROPAGATE_WINDINGS result\n";
2248  for (int i = 0; i < cinfo.tot_cell(); ++i) {
2249  std::cout << i << ": " << cinfo.cell(i) << "\n";
2250  }
2251  }
2252 }
2253 
2263 static bool apply_bool_op(BoolOpType bool_optype, const Array<int> &winding)
2264 {
2265  int nw = winding.size();
2266  BLI_assert(nw > 0);
2267  switch (bool_optype) {
2268  case BoolOpType::Intersect: {
2269  for (int i = 0; i < nw; ++i) {
2270  if (winding[i] == 0) {
2271  return false;
2272  }
2273  }
2274  return true;
2275  }
2276  case BoolOpType::Union: {
2277  for (int i = 0; i < nw; ++i) {
2278  if (winding[i] != 0) {
2279  return true;
2280  }
2281  }
2282  return false;
2283  }
2284  case BoolOpType::Difference: {
2285  /* if nw > 2, make it shape 0 minus the union of the rest. */
2286  if (winding[0] == 0) {
2287  return false;
2288  }
2289  if (nw == 1) {
2290  return true;
2291  }
2292  for (int i = 1; i < nw; ++i) {
2293  if (winding[i] >= 1) {
2294  return false;
2295  }
2296  }
2297  return true;
2298  }
2299  default:
2300  return false;
2301  }
2302 }
2303 
2309 static void extract_zero_volume_cell_tris(Vector<Face *> &r_tris,
2310  const IMesh &tm_subdivided,
2311  const PatchesInfo &pinfo,
2312  const CellsInfo &cinfo,
2313  IMeshArena *arena)
2314 {
2315  int dbg_level = 0;
2316  if (dbg_level > 0) {
2317  std::cout << "extract_zero_volume_cell_tris\n";
2318  }
2319  /* Find patches that are adjacent to zero-volume cells. */
2320  Array<bool> adj_to_zv(pinfo.tot_patch());
2321  for (int p : pinfo.index_range()) {
2322  const Patch &patch = pinfo.patch(p);
2323  if (cinfo.cell(patch.cell_above).zero_volume() || cinfo.cell(patch.cell_below).zero_volume()) {
2324  adj_to_zv[p] = true;
2325  }
2326  else {
2327  adj_to_zv[p] = false;
2328  }
2329  }
2330  /* Partition the adj_to_zv patches into stacks. */
2331  Vector<Vector<int>> patch_stacks;
2332  Array<bool> allocated_to_stack(pinfo.tot_patch(), false);
2333  for (int p : pinfo.index_range()) {
2334  if (!adj_to_zv[p] || allocated_to_stack[p]) {
2335  continue;
2336  }
2337  int stack_index = patch_stacks.size();
2338  patch_stacks.append(Vector<int>{p});
2339  Vector<int> &stack = patch_stacks[stack_index];
2340  Vector<bool> flipped{false};
2341  allocated_to_stack[p] = true;
2342  /* We arbitrarily choose p's above and below directions as above and below for whole stack.
2343  * Triangles in the stack that don't follow that convention are marked with flipped = true.
2344  * The non-zero-volume cell above the whole stack, following this convention, is
2345  * above_stack_cell. The non-zero-volume cell below the whole stack is #below_stack_cell. */
2346  /* First, walk in the above_cell direction from p. */
2347  int pwalk = p;
2348  const Patch *pwalk_patch = &pinfo.patch(pwalk);
2349  int c = pwalk_patch->cell_above;
2350  const Cell *cell = &cinfo.cell(c);
2351  while (cell->zero_volume()) {
2352  /* In zero-volume cells, the cell should have exactly two patches. */
2353  BLI_assert(cell->patches().size() == 2);
2354  int pother = cell->patch_other(pwalk);
2355  bool flip = pinfo.patch(pother).cell_above == c;
2356  flipped.append(flip);
2357  stack.append(pother);
2358  allocated_to_stack[pother] = true;
2359  pwalk = pother;
2360  pwalk_patch = &pinfo.patch(pwalk);
2361  c = flip ? pwalk_patch->cell_below : pwalk_patch->cell_above;
2362  cell = &cinfo.cell(c);
2363  }
2364  const Cell *above_stack_cell = cell;
2365  /* Now walk in the below_cell direction from p. */
2366  pwalk = p;
2367  pwalk_patch = &pinfo.patch(pwalk);
2368  c = pwalk_patch->cell_below;
2369  cell = &cinfo.cell(c);
2370  while (cell->zero_volume()) {
2371  BLI_assert(cell->patches().size() == 2);
2372  int pother = cell->patch_other(pwalk);
2373  bool flip = pinfo.patch(pother).cell_below == c;
2374  flipped.append(flip);
2375  stack.append(pother);
2376  allocated_to_stack[pother] = true;
2377  pwalk = pother;
2378  pwalk_patch = &pinfo.patch(pwalk);
2379  c = flip ? pwalk_patch->cell_above : pwalk_patch->cell_below;
2380  cell = &cinfo.cell(c);
2381  }
2382  const Cell *below_stack_cell = cell;
2383  if (dbg_level > 0) {
2384  std::cout << "process zero-volume patch stack " << stack << "\n";
2385  std::cout << "flipped = ";
2386  for (bool b : flipped) {
2387  std::cout << b << " ";
2388  }
2389  std::cout << "\n";
2390  }
2391  if (above_stack_cell->in_output_volume() ^ below_stack_cell->in_output_volume()) {
2392  bool need_flipped_tri = above_stack_cell->in_output_volume();
2393  if (dbg_level > 0) {
2394  std::cout << "need tri " << (need_flipped_tri ? "flipped" : "") << "\n";
2395  }
2396  int t_to_add = NO_INDEX;
2397  for (int i : stack.index_range()) {
2398  if (flipped[i] == need_flipped_tri) {
2399  t_to_add = pinfo.patch(stack[i]).tri(0);
2400  if (dbg_level > 0) {
2401  std::cout << "using tri " << t_to_add << "\n";
2402  }
2403  r_tris.append(tm_subdivided.face(t_to_add));
2404  break;
2405  }
2406  }
2407  if (t_to_add == NO_INDEX) {
2408  const Face *f = tm_subdivided.face(pinfo.patch(p).tri(0));
2409  const Face &tri = *f;
2410  /* We need flipped version or else we would have found it above. */
2411  std::array<const Vert *, 3> flipped_vs = {tri[0], tri[2], tri[1]};
2412  std::array<int, 3> flipped_e_origs = {
2413  tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]};
2414  std::array<bool, 3> flipped_is_intersect = {
2415  tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]};
2416  Face *flipped_f = arena->add_face(
2417  flipped_vs, f->orig, flipped_e_origs, flipped_is_intersect);
2418  r_tris.append(flipped_f);
2419  }
2420  }
2421  }
2422 }
2423 
2435 static IMesh extract_from_in_output_volume_diffs(const IMesh &tm_subdivided,
2436  const PatchesInfo &pinfo,
2437  const CellsInfo &cinfo,
2438  IMeshArena *arena)
2439 {
2440  constexpr int dbg_level = 0;
2441  if (dbg_level > 0) {
2442  std::cout << "\nEXTRACT_FROM_FLAG_DIFFS\n";
2443  }
2444  Vector<Face *> out_tris;
2445  out_tris.reserve(tm_subdivided.face_size());
2446  bool any_zero_volume_cell = false;
2447  for (int t : tm_subdivided.face_index_range()) {
2448  int p = pinfo.tri_patch(t);
2449  const Patch &patch = pinfo.patch(p);
2450  const Cell &cell_above = cinfo.cell(patch.cell_above);
2451  const Cell &cell_below = cinfo.cell(patch.cell_below);
2452  if (dbg_level > 0) {
2453  std::cout << "tri " << t << ": cell_above=" << patch.cell_above
2454  << " cell_below=" << patch.cell_below << "\n";
2455  std::cout << " in_output_volume_above=" << cell_above.in_output_volume()
2456  << " in_output_volume_below=" << cell_below.in_output_volume() << "\n";
2457  }
2458  bool adjacent_zero_volume_cell = cell_above.zero_volume() || cell_below.zero_volume();
2459  any_zero_volume_cell |= adjacent_zero_volume_cell;
2460  if (cell_above.in_output_volume() ^ cell_below.in_output_volume() &&
2461  !adjacent_zero_volume_cell) {
2462  bool flip = cell_above.in_output_volume();
2463  if (dbg_level > 0) {
2464  std::cout << "need tri " << t << " flip=" << flip << "\n";
2465  }
2466  Face *f = tm_subdivided.face(t);
2467  if (flip) {
2468  Face &tri = *f;
2469  std::array<const Vert *, 3> flipped_vs = {tri[0], tri[2], tri[1]};
2470  std::array<int, 3> flipped_e_origs = {
2471  tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]};
2472  std::array<bool, 3> flipped_is_intersect = {
2473  tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]};
2474  Face *flipped_f = arena->add_face(
2475  flipped_vs, f->orig, flipped_e_origs, flipped_is_intersect);
2476  out_tris.append(flipped_f);
2477  }
2478  else {
2479  out_tris.append(f);
2480  }
2481  }
2482  }
2483  if (any_zero_volume_cell) {
2484  extract_zero_volume_cell_tris(out_tris, tm_subdivided, pinfo, cinfo, arena);
2485  }
2486  return IMesh(out_tris);
2487 }
2488 
2489 static const char *bool_optype_name(BoolOpType op)
2490 {
2491  switch (op) {
2492  case BoolOpType::None:
2493  return "none";
2494  case BoolOpType::Intersect:
2495  return "intersect";
2496  case BoolOpType::Union:
2497  return "union";
2498  case BoolOpType::Difference:
2499  return "difference";
2500  default:
2501  return "<unknown>";
2502  }
2503 }
2504 
2505 static double3 calc_point_inside_tri_db(const Face &tri)
2506 {
2507  const Vert *v0 = tri.vert[0];
2508  const Vert *v1 = tri.vert[1];
2509  const Vert *v2 = tri.vert[2];
2510  double3 ans = v0->co / 3 + v1->co / 3 + v2->co / 3;
2511  return ans;
2512 }
2513 class InsideShapeTestData {
2514  public:
2515  const IMesh &tm;
2516  std::function<int(int)> shape_fn;
2517  int nshapes;
2518  /* A per-shape vector of parity of hits of that shape. */
2519  Array<int> hit_parity;
2520 
2521  InsideShapeTestData(const IMesh &tm, std::function<int(int)> shape_fn, int nshapes)
2522  : tm(tm), shape_fn(shape_fn), nshapes(nshapes)
2523  {
2524  }
2525 };
2526 
2527 static void inside_shape_callback(void *userdata,
2528  int index,
2529  const BVHTreeRay *ray,
2530  BVHTreeRayHit *UNUSED(hit))
2531 {
2532  const int dbg_level = 0;
2533  if (dbg_level > 0) {
2534  std::cout << "inside_shape_callback, index = " << index << "\n";
2535  }
2536  InsideShapeTestData *data = static_cast<InsideShapeTestData *>(userdata);
2537  const Face &tri = *data->tm.face(index);
2538  int shape = data->shape_fn(tri.orig);
2539  if (shape == -1) {
2540  return;
2541  }
2542  float dist;
2543  float fv0[3];
2544  float fv1[3];
2545  float fv2[3];
2546  for (int i = 0; i < 3; ++i) {
2547  fv0[i] = float(tri.vert[0]->co[i]);
2548  fv1[i] = float(tri.vert[1]->co[i]);
2549  fv2[i] = float(tri.vert[2]->co[i]);
2550  }
2551  if (dbg_level > 0) {
2552  std::cout << " fv0=(" << fv0[0] << "," << fv0[1] << "," << fv0[2] << ")\n";
2553  std::cout << " fv1=(" << fv1[0] << "," << fv1[1] << "," << fv1[2] << ")\n";
2554  std::cout << " fv2=(" << fv2[0] << "," << fv2[1] << "," << fv2[2] << ")\n";
2555  }
2557  ray->origin, ray->direction, fv0, fv1, fv2, &dist, nullptr, FLT_EPSILON)) {
2558  /* Count parity as +1 if ray is in the same direction as tri's normal,
2559  * and -1 if the directions are opposite. */
2560  double3 o_db{double(ray->origin[0]), double(ray->origin[1]), double(ray->origin[2])};
2561  int parity = orient3d(tri.vert[0]->co, tri.vert[1]->co, tri.vert[2]->co, o_db);
2562  if (dbg_level > 0) {
2563  std::cout << "origin at " << o_db << ", parity = " << parity << "\n";
2564  }
2565  data->hit_parity[shape] += parity;
2566  }
2567 }
2568 
2578 static void test_tri_inside_shapes(const IMesh &tm,
2579  std::function<int(int)> shape_fn,
2580  int nshapes,
2581  int test_t_index,
2582  BVHTree *tree,
2583  Array<float> &in_shape)
2584 {
2585  const int dbg_level = 0;
2586  if (dbg_level > 0) {
2587  std::cout << "test_point_inside_shapes, t_index = " << test_t_index << "\n";
2588  }
2589  Face &tri_test = *tm.face(test_t_index);
2590  int shape = shape_fn(tri_test.orig);
2591  if (shape == -1) {
2592  in_shape.fill(0.0f);
2593  return;
2594  }
2595  double3 test_point = calc_point_inside_tri_db(tri_test);
2596  /* Offset the test point a tiny bit in the tri_test normal direction. */
2597  tri_test.populate_plane(false);
2598  double3 norm = math::normalize(tri_test.plane->norm);
2599  const double offset_amount = 1e-5;
2600  double3 offset_test_point = test_point + offset_amount * norm;
2601  if (dbg_level > 0) {
2602  std::cout << "test tri is in shape " << shape << "\n";
2603  std::cout << "test point = " << test_point << "\n";
2604  std::cout << "offset_test_point = " << offset_test_point << "\n";
2605  }
2606  /* Try six test rays almost along orthogonal axes.
2607  * Perturb their directions slightly to make it less likely to hit a seam.
2608  * Ray-cast assumes they have unit length, so use r1 near 1 and
2609  * ra near 0.5, and rb near .01, but normalized so `sqrt(r1^2 + ra^2 + rb^2) == 1`. */
2610  constexpr int rays_num = 6;
2611  constexpr float r1 = 0.9987025295199663f;
2612  constexpr float ra = 0.04993512647599832f;
2613  constexpr float rb = 0.009987025295199663f;
2614  const float test_rays[rays_num][3] = {
2615  {r1, ra, rb}, {-r1, -ra, -rb}, {rb, r1, ra}, {-rb, -r1, -ra}, {ra, rb, r1}, {-ra, -rb, -r1}};
2616  InsideShapeTestData data(tm, shape_fn, nshapes);
2617  data.hit_parity = Array<int>(nshapes, 0);
2618  Array<int> count_insides(nshapes, 0);
2619  const float co[3] = {
2620  float(offset_test_point[0]), float(offset_test_point[1]), float(offset_test_point[2])};
2621  for (int i = 0; i < rays_num; ++i) {
2622  if (dbg_level > 0) {
2623  std::cout << "shoot ray " << i << "(" << test_rays[i][0] << "," << test_rays[i][1] << ","
2624  << test_rays[i][2] << ")\n";
2625  }
2626  BLI_bvhtree_ray_cast_all(tree, co, test_rays[i], 0.0f, FLT_MAX, inside_shape_callback, &data);
2627  if (dbg_level > 0) {
2628  std::cout << "ray " << i << " result:";
2629  for (int j = 0; j < nshapes; ++j) {
2630  std::cout << " " << data.hit_parity[j];
2631  }
2632  std::cout << "\n";
2633  }
2634  for (int j = 0; j < nshapes; ++j) {
2635  if (j != shape && data.hit_parity[j] > 0) {
2636  ++count_insides[j];
2637  }
2638  }
2639  data.hit_parity.fill(0);
2640  }
2641  for (int j = 0; j < nshapes; ++j) {
2642  if (j == shape) {
2643  in_shape[j] = 1.0f; /* Let's say a shape is always inside itself. */
2644  }
2645  else {
2646  in_shape[j] = float(count_insides[j]) / float(rays_num);
2647  }
2648  if (dbg_level > 0) {
2649  std::cout << "shape " << j << " inside = " << in_shape[j] << "\n";
2650  }
2651  }
2652 }
2653 
2660 static BVHTree *raycast_tree(const IMesh &tm)
2661 {
2662  BVHTree *tree = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 4, 6);
2663  for (int i : tm.face_index_range()) {
2664  const Face *f = tm.face(i);
2665  float t_cos[9];
2666  for (int j = 0; j < 3; ++j) {
2667  const Vert *v = f->vert[j];
2668  for (int k = 0; k < 3; ++k) {
2669  t_cos[3 * j + k] = float(v->co[k]);
2670  }
2671  }
2672  BLI_bvhtree_insert(tree, i, t_cos, 3);
2673  }
2675  return tree;
2676 }
2677 
2682 static bool raycast_test_remove(BoolOpType op, Array<int> &winding, int shape, bool *r_do_flip)
2683 {
2684  constexpr int dbg_level = 0;
2685  /* Find out the "in the output volume" flag for each of the cases of winding[shape] == 0
2686  * and winding[shape] == 1. If the flags are different, this patch should be in the output.
2687  * Also, if this is a Difference and the shape isn't the first one, need to flip the normals.
2688  */
2689  winding[shape] = 0;
2690  bool in_output_volume_0 = apply_bool_op(op, winding);
2691  winding[shape] = 1;
2692  bool in_output_volume_1 = apply_bool_op(op, winding);
2693  bool do_remove = in_output_volume_0 == in_output_volume_1;
2694  bool do_flip = !do_remove && op == BoolOpType::Difference && shape != 0;
2695  if (dbg_level > 0) {
2696  std::cout << "winding = ";
2697  for (int i = 0; i < winding.size(); ++i) {
2698  std::cout << winding[i] << " ";
2699  }
2700  std::cout << "\niv0=" << in_output_volume_0 << ", iv1=" << in_output_volume_1 << "\n";
2701  std::cout << " remove=" << do_remove << ", flip=" << do_flip << "\n";
2702  }
2703  *r_do_flip = do_flip;
2704  return do_remove;
2705 }
2706 
2708 static void raycast_add_flipped(Vector<Face *> &out_faces, Face &tri, IMeshArena *arena)
2709 {
2710 
2711  Array<const Vert *> flipped_vs = {tri[0], tri[2], tri[1]};
2712  Array<int> flipped_e_origs = {tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]};
2713  Array<bool> flipped_is_intersect = {
2714  tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]};
2715  Face *flipped_f = arena->add_face(flipped_vs, tri.orig, flipped_e_origs, flipped_is_intersect);
2716  out_faces.append(flipped_f);
2717 }
2718 
2727 static IMesh raycast_tris_boolean(const IMesh &tm,
2728  BoolOpType op,
2729  int nshapes,
2730  std::function<int(int)> shape_fn,
2731  IMeshArena *arena)
2732 {
2733  constexpr int dbg_level = 0;
2734  if (dbg_level > 0) {
2735  std::cout << "RAYCAST_TRIS_BOOLEAN\n";
2736  }
2737  IMesh ans;
2738  BVHTree *tree = raycast_tree(tm);
2739  Vector<Face *> out_faces;
2740  out_faces.reserve(tm.face_size());
2741 # ifdef WITH_TBB
2742  tbb::spin_mutex mtx;
2743 # endif
2744  const int grainsize = 256;
2745  threading::parallel_for(IndexRange(tm.face_size()), grainsize, [&](IndexRange range) {
2746  Array<float> in_shape(nshapes, 0);
2747  Array<int> winding(nshapes, 0);
2748  for (int t : range) {
2749  Face &tri = *tm.face(t);
2750  int shape = shape_fn(tri.orig);
2751  if (dbg_level > 0) {
2752  std::cout << "process triangle " << t << " = " << &tri << "\n";
2753  std::cout << "shape = " << shape << "\n";
2754  }
2755  test_tri_inside_shapes(tm, shape_fn, nshapes, t, tree, in_shape);
2756  for (int other_shape = 0; other_shape < nshapes; ++other_shape) {
2757  if (other_shape == shape) {
2758  continue;
2759  }
2760  /* The in_shape array has a confidence value for "insideness".
2761  * For most operations, even a hint of being inside
2762  * gives good results, but when shape is a cutter in a Difference
2763  * operation, we want to be pretty sure that the point is inside other_shape.
2764  * E.g., T75827.
2765  * Also, when the operation is intersection, we also want high confidence.
2766  */
2767  bool need_high_confidence = (op == BoolOpType::Difference && shape != 0) ||
2768  op == BoolOpType::Intersect;
2769  bool inside = in_shape[other_shape] >= (need_high_confidence ? 0.5f : 0.1f);
2770  if (dbg_level > 0) {
2771  std::cout << "test point is " << (inside ? "inside" : "outside") << " other_shape "
2772  << other_shape << " val = " << in_shape[other_shape] << "\n";
2773  }
2774  winding[other_shape] = inside;
2775  }
2776  bool do_flip;
2777  bool do_remove = raycast_test_remove(op, winding, shape, &do_flip);
2778  {
2779 # ifdef WITH_TBB
2780  tbb::spin_mutex::scoped_lock lock(mtx);
2781 # endif
2782  if (!do_remove) {
2783  if (!do_flip) {
2784  out_faces.append(&tri);
2785  }
2786  else {
2787  raycast_add_flipped(out_faces, tri, arena);
2788  }
2789  }
2790  }
2791  }
2792  });
2794  ans.set_faces(out_faces);
2795  return ans;
2796 }
2797 
2798 /* This is (sometimes much faster) version of raycast boolean
2799  * that does it per patch rather than per triangle.
2800  * It may fail in cases where raycast_tri_boolean will succeed,
2801  * but the latter can be very slow on huge meshes. */
2802 static IMesh raycast_patches_boolean(const IMesh &tm,
2803  BoolOpType op,
2804  int nshapes,
2805  std::function<int(int)> shape_fn,
2806  const PatchesInfo &pinfo,
2807  IMeshArena *arena)
2808 {
2809  constexpr int dbg_level = 0;
2810  if (dbg_level > 0) {
2811  std::cout << "RAYCAST_PATCHES_BOOLEAN\n";
2812  }
2813  IMesh ans;
2814  BVHTree *tree = raycast_tree(tm);
2815  Vector<Face *> out_faces;
2816  out_faces.reserve(tm.face_size());
2817  Array<float> in_shape(nshapes, 0);
2818  Array<int> winding(nshapes, 0);
2819  for (int p : pinfo.index_range()) {
2820  const Patch &patch = pinfo.patch(p);
2821  /* For test triangle, choose one in the middle of patch list
2822  * as the ones near the beginning may be very near other patches. */
2823  int test_t_index = patch.tri(patch.tot_tri() / 2);
2824  Face &tri_test = *tm.face(test_t_index);
2825  /* Assume all triangles in a patch are in the same shape. */
2826  int shape = shape_fn(tri_test.orig);
2827  if (dbg_level > 0) {
2828  std::cout << "process patch " << p << " = " << patch << "\n";
2829  std::cout << "test tri = " << test_t_index << " = " << &tri_test << "\n";
2830  std::cout << "shape = " << shape << "\n";
2831  }
2832  if (shape == -1) {
2833  continue;
2834  }
2835  test_tri_inside_shapes(tm, shape_fn, nshapes, test_t_index, tree, in_shape);
2836  for (int other_shape = 0; other_shape < nshapes; ++other_shape) {
2837  if (other_shape == shape) {
2838  continue;
2839  }
2840  bool need_high_confidence = (op == BoolOpType::Difference && shape != 0) ||
2841  op == BoolOpType::Intersect;
2842  bool inside = in_shape[other_shape] >= (need_high_confidence ? 0.5f : 0.1f);
2843  if (dbg_level > 0) {
2844  std::cout << "test point is " << (inside ? "inside" : "outside") << " other_shape "
2845  << other_shape << " val = " << in_shape[other_shape] << "\n";
2846  }
2847  winding[other_shape] = inside;
2848  }
2849  bool do_flip;
2850  bool do_remove = raycast_test_remove(op, winding, shape, &do_flip);
2851  if (!do_remove) {
2852  for (int t : patch.tris()) {
2853  Face *f = tm.face(t);
2854  if (!do_flip) {
2855  out_faces.append(f);
2856  }
2857  else {
2858  raycast_add_flipped(out_faces, *f, arena);
2859  }
2860  }
2861  }
2862  }
2864 
2865  ans.set_faces(out_faces);
2866  return ans;
2867 }
2873 static std::pair<int, int> find_tris_common_edge(const Face &tri1, const Face &tri2)
2874 {
2875  for (int i = 0; i < 3; ++i) {
2876  for (int j = 0; j < 3; ++j) {
2877  if (tri1[(i + 1) % 3] == tri2[j] && tri1[i] == tri2[(j + 1) % 3]) {
2878  return std::pair<int, int>(i, j);
2879  }
2880  }
2881  }
2882  return std::pair<int, int>(-1, -1);
2883 }
2884 
2885 struct MergeEdge {
2887  double len_squared = 0.0;
2888  /* v1 and v2 are the ends of the edge, ordered so that `v1->id < v2->id` */
2889  const Vert *v1 = nullptr;
2890  const Vert *v2 = nullptr;
2891  /* left_face and right_face are indices into #FaceMergeState.face. */
2892  int left_face = -1;
2893  int right_face = -1;
2894  int orig = -1; /* An edge orig index that can be used for this edge. */
2896  bool dissolvable = false;
2898  bool is_intersect = false;
2899 
2900  MergeEdge() = default;
2901 
2902  MergeEdge(const Vert *va, const Vert *vb)
2903  {
2904  if (va->id < vb->id) {
2905  this->v1 = va;
2906  this->v2 = vb;
2907  }
2908  else {
2909  this->v1 = vb;
2910  this->v2 = va;
2911  }
2912  };
2913 };
2914 
2915 struct MergeFace {
2917  Vector<const Vert *> vert;
2919  Vector<int> edge;
2921  int merge_to = -1;
2923  int orig = -1;
2924 };
2925 struct FaceMergeState {
2929  Vector<MergeFace> face;
2934  Vector<MergeEdge> edge;
2939  Map<std::pair<int, int>, int> edge_map;
2940 };
2941 
2942 static std::ostream &operator<<(std::ostream &os, const FaceMergeState &fms)
2943 {
2944  os << "faces:\n";
2945  for (int f : fms.face.index_range()) {
2946  const MergeFace &mf = fms.face[f];
2947  std::cout << f << ": orig=" << mf.orig << " verts ";
2948  for (const Vert *v : mf.vert) {
2949  std::cout << v << " ";
2950  }
2951  std::cout << "\n";
2952  std::cout << " edges " << mf.edge << "\n";
2953  std::cout << " merge_to = " << mf.merge_to << "\n";
2954  }
2955  os << "\nedges:\n";
2956  for (int e : fms.edge.index_range()) {
2957  const MergeEdge &me = fms.edge[e];
2958  std::cout << e << ": (" << me.v1 << "," << me.v2 << ") left=" << me.left_face
2959  << " right=" << me.right_face << " dis=" << me.dissolvable << " orig=" << me.orig
2960  << " is_int=" << me.is_intersect << "\n";
2961  }
2962  return os;
2963 }
2964 
2975 static void init_face_merge_state(FaceMergeState *fms,
2976  const Vector<int> &tris,
2977  const IMesh &tm,
2978  const double3 &norm)
2979 {
2980  constexpr int dbg_level = 0;
2981  /* Reserve enough faces and edges so that neither will have to resize. */
2982  fms->face.reserve(tris.size() + 1);
2983  fms->edge.reserve(3 * tris.size());
2984  fms->edge_map.reserve(3 * tris.size());
2985  if (dbg_level > 0) {
2986  std::cout << "\nINIT_FACE_MERGE_STATE\n";
2987  }
2988  for (int t : tris.index_range()) {
2989  MergeFace mf;
2990  const Face &tri = *tm.face(tris[t]);
2991  if (dbg_level > 0) {
2992  std::cout << "process tri = " << &tri << "\n";
2993  }
2994  BLI_assert(tri.plane_populated());
2995  if (math::dot(norm, tri.plane->norm) <= 0.0) {
2996  if (dbg_level > 0) {
2997  std::cout << "triangle has wrong orientation, skipping\n";
2998  }
2999  continue;
3000  }
3001  mf.vert.append(tri[0]);
3002  mf.vert.append(tri[1]);
3003  mf.vert.append(tri[2]);
3004  mf.orig = tri.orig;
3005  int f = fms->face.append_and_get_index(mf);
3006  if (dbg_level > 1) {
3007  std::cout << "appended MergeFace for tri at f = " << f << "\n";
3008  }
3009  for (int i = 0; i < 3; ++i) {
3010  int inext = (i + 1) % 3;
3011  MergeEdge new_me(mf.vert[i], mf.vert[inext]);
3012  std::pair<int, int> canon_vs(new_me.v1->id, new_me.v2->id);
3013  int me_index = fms->edge_map.lookup_default(canon_vs, -1);
3014  if (dbg_level > 1) {
3015  std::cout << "new_me = canon_vs = " << new_me.v1 << ", " << new_me.v2 << "\n";
3016  std::cout << "me_index lookup = " << me_index << "\n";
3017  }
3018  if (me_index == -1) {
3019  double3 vec = new_me.v2->co - new_me.v1->co;
3020  new_me.len_squared = math::length_squared(vec);
3021  new_me.orig = tri.edge_orig[i];
3022  new_me.is_intersect = tri.is_intersect[i];
3023  new_me.dissolvable = (new_me.orig == NO_INDEX && !new_me.is_intersect);
3024  fms->edge.append(new_me);
3025  me_index = fms->edge.size() - 1;
3026  fms->edge_map.add_new(canon_vs, me_index);
3027  if (dbg_level > 1) {
3028  std::cout << "added new me with me_index = " << me_index << "\n";
3029  std::cout << " len_squared = " << new_me.len_squared << " orig = " << new_me.orig
3030  << ", is_intersect" << new_me.is_intersect
3031  << ", dissolvable = " << new_me.dissolvable << "\n";
3032  }
3033  }
3034  MergeEdge &me = fms->edge[me_index];
3035  if (dbg_level > 1) {
3036  std::cout << "retrieved me at index " << me_index << ":\n";
3037  std::cout << " v1 = " << me.v1 << " v2 = " << me.v2 << "\n";
3038  std::cout << " dis = " << me.dissolvable << " int = " << me.is_intersect << "\n";
3039  std::cout << " left_face = " << me.left_face << " right_face = " << me.right_face << "\n";
3040  }
3041  if (me.dissolvable && tri.edge_orig[i] != NO_INDEX) {
3042  if (dbg_level > 1) {
3043  std::cout << "reassigning orig to " << tri.edge_orig[i] << ", dissolvable = false\n";
3044  }
3045  me.dissolvable = false;
3046  me.orig = tri.edge_orig[i];
3047  }
3048  if (me.dissolvable && tri.is_intersect[i]) {
3049  if (dbg_level > 1) {
3050  std::cout << "reassigning dissolvable = false, is_intersect = true\n";
3051  }
3052  me.dissolvable = false;
3053  me.is_intersect = true;
3054  }
3055  /* This face is left or right depending on orientation of edge. */
3056  if (me.v1 == mf.vert[i]) {
3057  if (dbg_level > 1) {
3058  std::cout << "me.v1 == mf.vert[i] so set edge[" << me_index << "].left_face = " << f
3059  << "\n";
3060  }
3061  if (me.left_face != -1) {
3062  /* Unexpected in the normal case: this means more than one triangle shares this
3063  * edge in the same orientation. But be tolerant of this case. By making this
3064  * edge not dissolvable, we'll avoid future problems due to this non-manifold topology.
3065  */
3066  if (dbg_level > 1) {
3067  std::cout << "me.left_face was already occupied, so triangulation wasn't good\n";
3068  }
3069  me.dissolvable = false;
3070  }
3071  else {
3072  fms->edge[me_index].left_face = f;
3073  }
3074  }
3075  else {
3076  if (dbg_level > 1) {
3077  std::cout << "me.v1 != mf.vert[i] so set edge[" << me_index << "].right_face = " << f
3078  << "\n";
3079  }
3080  if (me.right_face != -1) {
3081  /* Unexpected, analogous to the me.left_face != -1 case above. */
3082  if (dbg_level > 1) {
3083  std::cout << "me.right_face was already occupied, so triangulation wasn't good\n";
3084  }
3085  me.dissolvable = false;
3086  }
3087  else {
3088  fms->edge[me_index].right_face = f;
3089  }
3090  }
3091  fms->face[f].edge.append(me_index);
3092  }
3093  }
3094  if (dbg_level > 0) {
3095  std::cout << *fms;
3096  }
3097 }
3098 
3105 static bool dissolve_leaves_valid_bmesh(FaceMergeState *fms,
3106  const MergeEdge &me,
3107  int me_index,
3108  const MergeFace &mf_left,
3109  const MergeFace &mf_right)
3110 {
3111  int a_edge_start = mf_left.edge.first_index_of_try(me_index);
3112  BLI_assert(a_edge_start != -1);
3113  int alen = mf_left.vert.size();
3114  int blen = mf_right.vert.size();
3115  int b_left_face = me.right_face;
3116  bool ok = true;
3117  /* Is there another edge, not me, in A's face, whose right face is B's left? */
3118  for (int a_e_index = (a_edge_start + 1) % alen; ok && a_e_index != a_edge_start;
3119  a_e_index = (a_e_index + 1) % alen) {
3120  const MergeEdge &a_me_cur = fms->edge[mf_left.edge[a_e_index]];
3121  if (a_me_cur.right_face == b_left_face) {
3122  ok = false;
3123  }
3124  }
3125  /* Is there a vert in A, not me.v1 or me.v2, that is also in B?
3126  * One could avoid this O(n^2) algorithm if had a structure
3127  * saying which faces a vertex touches. */
3128  for (int a_v_index = 0; ok && a_v_index < alen; ++a_v_index) {
3129  const Vert *a_v = mf_left.vert[a_v_index];
3130  if (!ELEM(a_v, me.v1, me.v2)) {
3131  for (int b_v_index = 0; b_v_index < blen; ++b_v_index) {
3132  const Vert *b_v = mf_right.vert[b_v_index];
3133  if (a_v == b_v) {
3134  ok = false;
3135  }
3136  }
3137  }
3138  }
3139  return ok;
3140 }
3141 
3149 static void splice_faces(
3150  FaceMergeState *fms, MergeEdge &me, int me_index, MergeFace &mf_left, MergeFace &mf_right)
3151 {
3152  int a_edge_start = mf_left.edge.first_index_of_try(me_index);
3153  int b_edge_start = mf_right.edge.first_index_of_try(me_index);
3154  BLI_assert(a_edge_start != -1 && b_edge_start != -1);
3155  int alen = mf_left.vert.size();
3156  int blen = mf_right.vert.size();
3157  Vector<const Vert *> splice_vert;
3158  Vector<int> splice_edge;
3159  splice_vert.reserve(alen + blen - 2);
3160  splice_edge.reserve(alen + blen - 2);
3161  int ai = 0;
3162  while (ai < a_edge_start) {
3163  splice_vert.append(mf_left.vert[ai]);
3164  splice_edge.append(mf_left.edge[ai]);
3165  ++ai;
3166  }
3167  int bi = b_edge_start + 1;
3168  while (bi != b_edge_start) {
3169  if (bi >= blen) {
3170  bi = 0;
3171  if (bi == b_edge_start) {
3172  break;
3173  }
3174  }
3175  splice_vert.append(mf_right.vert[bi]);
3176  splice_edge.append(mf_right.edge[bi]);
3177  if (mf_right.vert[bi] == fms->edge[mf_right.edge[bi]].v1) {
3178  fms->edge[mf_right.edge[bi]].left_face = me.left_face;
3179  }
3180  else {
3181  fms->edge[mf_right.edge[bi]].right_face = me.left_face;
3182  }
3183  ++bi;
3184  }
3185  ai = a_edge_start + 1;
3186  while (ai < alen) {
3187  splice_vert.append(mf_left.vert[ai]);
3188  splice_edge.append(mf_left.edge[ai]);
3189  ++ai;
3190  }
3191  mf_right.merge_to = me.left_face;
3192  mf_left.vert = splice_vert;
3193  mf_left.edge = splice_edge;
3194  me.left_face = -1;
3195  me.right_face = -1;
3196 }
3197 
3205 static void do_dissolve(FaceMergeState *fms)
3206 {
3207  const int dbg_level = 0;
3208  if (dbg_level > 1) {
3209  std::cout << "\nDO_DISSOLVE\n";
3210  }
3211  Vector<int> dissolve_edges;
3212  for (int e : fms->edge.index_range()) {
3213  if (fms->edge[e].dissolvable) {
3214  dissolve_edges.append(e);
3215  }
3216  }
3217  if (dissolve_edges.size() == 0) {
3218  return;
3219  }
3220  /* Things look nicer if we dissolve the longer edges first. */
3221  std::sort(
3222  dissolve_edges.begin(), dissolve_edges.end(), [fms](const int &a, const int &b) -> bool {
3223  return (fms->edge[a].len_squared > fms->edge[b].len_squared);
3224  });
3225  if (dbg_level > 0) {
3226  std::cout << "Sorted dissolvable edges: " << dissolve_edges << "\n";
3227  }
3228  for (int me_index : dissolve_edges) {
3229  MergeEdge &me = fms->edge[me_index];
3230  if (me.left_face == -1 || me.right_face == -1) {
3231  continue;
3232  }
3233  MergeFace &mf_left = fms->face[me.left_face];
3234  MergeFace &mf_right = fms->face[me.right_face];
3235  if (!dissolve_leaves_valid_bmesh(fms, me, me_index, mf_left, mf_right)) {
3236  continue;
3237  }
3238  if (dbg_level > 0) {
3239  std::cout << "Removing edge " << me_index << "\n";
3240  }
3241  splice_faces(fms, me, me_index, mf_left, mf_right);
3242  if (dbg_level > 1) {
3243  std::cout << "state after removal:\n";
3244  std::cout << *fms;
3245  }
3246  }
3247 }
3248 
3259 static Vector<Face *> merge_tris_for_face(Vector<int> tris,
3260  const IMesh &tm,
3261  const IMesh &imesh_in,
3262  IMeshArena *arena)
3263 {
3264  constexpr int dbg_level = 0;
3265  if (dbg_level > 0) {
3266  std::cout << "merge_tris_for_face\n";
3267  std::cout << "tris: " << tris << "\n";
3268  }
3269  Vector<Face *> ans;
3270  if (tris.size() <= 1) {
3271  if (tris.size() == 1) {
3272  ans.append(tm.face(tris[0]));
3273  }
3274  return ans;
3275  }
3276  bool done = false;
3277  double3 first_tri_normal = tm.face(tris[0])->plane->norm;
3278  double3 second_tri_normal = tm.face(tris[1])->plane->norm;
3279  if (tris.size() == 2 && math::dot(first_tri_normal, second_tri_normal) > 0.0) {
3280  /* Is this a case where quad with one diagonal remained unchanged?
3281  * Worth special handling because this case will be very common. */
3282  Face &tri1 = *tm.face(tris[0]);
3283  Face &tri2 = *tm.face(tris[1]);
3284  Face *in_face = imesh_in.face(tri1.orig);
3285  if (in_face->size() == 4) {
3286  std::pair<int, int> estarts = find_tris_common_edge(tri1, tri2);
3287  if (estarts.first != -1 && tri1.edge_orig[estarts.first] == NO_INDEX) {
3288  if (dbg_level > 0) {
3289  std::cout << "try recovering orig quad case\n";
3290  std::cout << "tri1 = " << &tri1 << "\n";
3291  std::cout << "tri1 = " << &tri2 << "\n";
3292  }
3293  int i0 = estarts.first;
3294  int i1 = (i0 + 1) % 3;
3295  int i2 = (i0 + 2) % 3;
3296  int j2 = (estarts.second + 2) % 3;
3297  Face tryface({tri1[i1], tri1[i2], tri1[i0], tri2[j2]}, -1, -1, {}, {});
3298  if (tryface.cyclic_equal(*in_face)) {
3299  if (dbg_level > 0) {
3300  std::cout << "inface = " << in_face << "\n";
3301  std::cout << "quad recovery worked\n";
3302  }
3303  ans.append(in_face);
3304  done = true;
3305  }
3306  }
3307  }
3308  }
3309  if (done) {
3310  return ans;
3311  }
3312 
3313  double3 first_tri_normal_rev = -first_tri_normal;
3314  for (const double3 &norm : {first_tri_normal, first_tri_normal_rev}) {
3315  FaceMergeState fms;
3316  init_face_merge_state(&fms, tris, tm, norm);
3317  do_dissolve(&fms);
3318  if (dbg_level > 0) {
3319  std::cout << "faces in merged result:\n";
3320  }
3321  for (const MergeFace &mf : fms.face) {
3322  if (mf.merge_to == -1) {
3323  Array<int> e_orig(mf.edge.size());
3324  Array<bool> is_intersect(mf.edge.size());
3325  for (int i : mf.edge.index_range()) {
3326  e_orig[i] = fms.edge[mf.edge[i]].orig;
3327  is_intersect[i] = fms.edge[mf.edge[i]].is_intersect;
3328  }
3329  Face *facep = arena->add_face(mf.vert, mf.orig, e_orig, is_intersect);
3330  ans.append(facep);
3331  if (dbg_level > 0) {
3332  std::cout << " " << facep << "\n";
3333  }
3334  }
3335  }
3336  }
3337  return ans;
3338 }
3339 
3340 static bool approx_in_line(const double3 &a, const double3 &b, const double3 &c)
3341 {
3342  double3 vec1 = b - a;
3343  double3 vec2 = c - b;
3344  double cos_ang = math::dot(math::normalize(vec1), math::normalize(vec2));
3345  return fabs(cos_ang - 1.0) < 1e-4;
3346 }
3347 
3354 static Array<bool> find_dissolve_verts(IMesh &imesh_out, int *r_count_dissolve)
3355 {
3356  imesh_out.populate_vert();
3357  /* dissolve[i] will say whether imesh_out.vert(i) can be dissolved. */
3358  Array<bool> dissolve(imesh_out.vert_size());
3359  for (int v_index : imesh_out.vert_index_range()) {
3360  const Vert &vert = *imesh_out.vert(v_index);
3361  dissolve[v_index] = (vert.orig == NO_INDEX);
3362  }
3363  /* neighbors[i] will be a pair giving the up-to-two neighboring vertices
3364  * of the vertex v in position i of imesh_out.vert.
3365  * If we encounter a third, then v will not be dissolvable. */
3366  Array<std::pair<const Vert *, const Vert *>> neighbors(
3367  imesh_out.vert_size(), std::pair<const Vert *, const Vert *>(nullptr, nullptr));
3368  for (int f : imesh_out.face_index_range()) {
3369  const Face &face = *imesh_out.face(f);
3370  for (int i : face.index_range()) {
3371  const Vert *v = face[i];
3372  int v_index = imesh_out.lookup_vert(v);
3373  BLI_assert(v_index != NO_INDEX);
3374  if (dissolve[v_index]) {
3375  const Vert *n1 = face[face.next_pos(i)];
3376  const Vert *n2 = face[face.prev_pos(i)];
3377  const Vert *f_n1 = neighbors[v_index].first;
3378  const Vert *f_n2 = neighbors[v_index].second;
3379  if (f_n1 != nullptr) {
3380  /* Already has a neighbor in another face; can't dissolve unless they are the same. */
3381  if (!((n1 == f_n2 && n2 == f_n1) || (n1 == f_n1 && n2 == f_n2))) {
3382  /* Different neighbors, so can't dissolve. */
3383  dissolve[v_index] = false;
3384  }
3385  }
3386  else {
3387  /* These are the first-seen neighbors. */
3388  neighbors[v_index] = std::pair<const Vert *, const Vert *>(n1, n2);
3389  }
3390  }
3391  }
3392  }
3393  int count = 0;
3394  for (int v_out : imesh_out.vert_index_range()) {
3395  if (dissolve[v_out]) {
3396  dissolve[v_out] = false; /* Will set back to true if final condition is satisfied. */
3397  const std::pair<const Vert *, const Vert *> &nbrs = neighbors[v_out];
3398  if (nbrs.first != nullptr) {
3399  BLI_assert(nbrs.second != nullptr);
3400  const Vert *v_v_out = imesh_out.vert(v_out);
3401  if (approx_in_line(nbrs.first->co, v_v_out->co, nbrs.second->co)) {
3402  dissolve[v_out] = true;
3403  ++count;
3404  }
3405  }
3406  }
3407  }
3408  if (r_count_dissolve != nullptr) {
3409  *r_count_dissolve = count;
3410  }
3411  return dissolve;
3412 }
3413 
3419 static void dissolve_verts(IMesh *imesh, const Array<bool> dissolve, IMeshArena *arena)
3420 {
3421  constexpr int inline_face_size = 100;
3422  Vector<bool, inline_face_size> face_pos_erase;
3423  bool any_faces_erased = false;
3424  for (int f : imesh->face_index_range()) {
3425  const Face &face = *imesh->face(f);
3426  face_pos_erase.clear();
3427  int erase_num = 0;
3428  for (const Vert *v : face) {
3429  int v_index = imesh->lookup_vert(v);
3430  BLI_assert(v_index != NO_INDEX);
3431  if (dissolve[v_index]) {
3432  face_pos_erase.append(true);
3433  ++erase_num;
3434  }
3435  else {
3436  face_pos_erase.append(false);
3437  }
3438  }
3439  if (erase_num > 0) {
3440  any_faces_erased |= imesh->erase_face_positions(f, face_pos_erase, arena);
3441  }
3442  }
3443  imesh->set_dirty_verts();
3444  if (any_faces_erased) {
3445  imesh->remove_null_faces();
3446  }
3447 }
3448 
3460 static IMesh polymesh_from_trimesh_with_dissolve(const IMesh &tm_out,
3461  const IMesh &imesh_in,
3462  IMeshArena *arena)
3463 {
3464  const int dbg_level = 0;
3465  if (dbg_level > 1) {
3466  std::cout << "\nPOLYMESH_FROM_TRIMESH_WITH_DISSOLVE\n";
3467  }
3468  /* For now: need plane normals for all triangles. */
3469  const int grainsize = 1024;
3470  threading::parallel_for(tm_out.face_index_range(), grainsize, [&](IndexRange range) {
3471  for (int i : range) {
3472  Face *tri = tm_out.face(i);
3473  tri->populate_plane(false);
3474  }
3475  });
3476  /* Gather all output triangles that are part of each input face.
3477  * face_output_tris[f] will be indices of triangles in tm_out
3478  * that have f as their original face. */
3479  int tot_in_face = imesh_in.face_size();
3480  Array<Vector<int>> face_output_tris(tot_in_face);
3481  for (int t : tm_out.face_index_range()) {
3482  const Face &tri = *tm_out.face(t);
3483  int in_face = tri.orig;
3484  face_output_tris[in_face].append(t);
3485  }
3486  if (dbg_level > 1) {
3487  std::cout << "face_output_tris:\n";
3488  for (int f : face_output_tris.index_range()) {
3489  std::cout << f << ": " << face_output_tris[f] << "\n";
3490  }
3491  }
3492 
3493  /* Merge triangles that we can from face_output_tri to make faces for output.
3494  * face_output_face[f] will be new original const Face *'s that
3495  * make up whatever part of the boolean output remains of input face f. */
3496  Array<Vector<Face *>> face_output_face(tot_in_face);
3497  int tot_out_face = 0;
3498  for (int in_f : imesh_in.face_index_range()) {
3499  if (dbg_level > 1) {
3500  std::cout << "merge tris for face " << in_f << "\n";
3501  }
3502  int out_tris_for_face_num = face_output_tris.size();
3503  if (out_tris_for_face_num == 0) {
3504  continue;
3505  }
3506  face_output_face[in_f] = merge_tris_for_face(face_output_tris[in_f], tm_out, imesh_in, arena);
3507  tot_out_face += face_output_face[in_f].size();
3508  }
3509  Array<Face *> face(tot_out_face);
3510  int out_f_index = 0;
3511  for (int in_f : imesh_in.face_index_range()) {
3512  const Vector<Face *> &f_faces = face_output_face[in_f];
3513  if (f_faces.size() > 0) {
3514  std::copy(f_faces.begin(), f_faces.end(), &face[out_f_index]);
3515  out_f_index += f_faces.size();
3516  }
3517  }
3518  IMesh imesh_out(face);
3519  /* Dissolve vertices that were (a) not original; and (b) now have valence 2 and
3520  * are between two other vertices that are exactly in line with them.
3521  * These were created because of triangulation edges that have been dissolved. */
3522  int count_dissolve;
3523  Array<bool> v_dissolve = find_dissolve_verts(imesh_out, &count_dissolve);
3524  if (count_dissolve > 0) {
3525  dissolve_verts(&imesh_out, v_dissolve, arena);
3526  }
3527  if (dbg_level > 1) {
3528  write_obj_mesh(imesh_out, "boolean_post_dissolve");
3529  }
3530 
3531  return imesh_out;
3532 }
3533 
3540 IMesh boolean_trimesh(IMesh &tm_in,
3541  BoolOpType op,
3542  int nshapes,
3543  std::function<int(int)> shape_fn,
3544  bool use_self,
3545  bool hole_tolerant,
3546  IMeshArena *arena)
3547 {
3548  constexpr int dbg_level = 0;
3549  if (dbg_level > 0) {
3550  std::cout << "BOOLEAN of " << nshapes << " operand" << (nshapes == 1 ? "" : "s")
3551  << " op=" << bool_optype_name(op) << "\n";
3552  if (dbg_level > 1) {
3553  tm_in.populate_vert();
3554  std::cout << "boolean_trimesh input:\n" << tm_in;
3555  write_obj_mesh(tm_in, "boolean_in");
3556  }
3557  }
3558  if (tm_in.face_size() == 0) {
3559  return IMesh(tm_in);
3560  }
3561 # ifdef PERFDEBUG
3562  double start_time = PIL_check_seconds_timer();
3563  std::cout << " boolean_trimesh, timing begins\n";
3564 # endif
3565 
3566  IMesh tm_si = trimesh_nary_intersect(tm_in, nshapes, shape_fn, use_self, arena);
3567  if (dbg_level > 1) {
3568  write_obj_mesh(tm_si, "boolean_tm_si");
3569  std::cout << "\nboolean_tm_input after intersection:\n" << tm_si;
3570  }
3571 # ifdef PERFDEBUG
3572  double intersect_time = PIL_check_seconds_timer();
3573  std::cout << " intersected, time = " << intersect_time - start_time << "\n";
3574 # endif
3575 
3576  /* It is possible for tm_si to be empty if all the input triangles are bogus/degenerate. */
3577  if (tm_si.face_size() == 0 || op == BoolOpType::None) {
3578  return tm_si;
3579  }
3580  auto si_shape_fn = [shape_fn, tm_si](int t) { return shape_fn(tm_si.face(t)->orig); };
3581  TriMeshTopology tm_si_topo(tm_si);
3582 # ifdef PERFDEBUG
3583  double topo_time = PIL_check_seconds_timer();
3584  std::cout << " topology built, time = " << topo_time - intersect_time << "\n";
3585 # endif
3586  bool pwn = is_pwn(tm_si, tm_si_topo);
3587 # ifdef PERFDEBUG
3588  double pwn_time = PIL_check_seconds_timer();
3589  std::cout << " pwn checked, time = " << pwn_time - topo_time << "\n";
3590 # endif
3591  IMesh tm_out;
3592  if (!pwn) {
3593  if (dbg_level > 0) {
3594  std::cout << "Input is not PWN, using raycast method\n";
3595  }
3596  if (hole_tolerant) {
3597  tm_out = raycast_tris_boolean(tm_si, op, nshapes, shape_fn, arena);
3598  }
3599  else {
3600  PatchesInfo pinfo = find_patches(tm_si, tm_si_topo);
3601  tm_out = raycast_patches_boolean(tm_si, op, nshapes, shape_fn, pinfo, arena);
3602  }
3603 # ifdef PERFDEBUG
3604  double raycast_time = PIL_check_seconds_timer();
3605  std::cout << " raycast_boolean done, time = " << raycast_time - pwn_time << "\n";
3606 # endif
3607  }
3608  else {
3609  PatchesInfo pinfo = find_patches(tm_si, tm_si_topo);
3610 # ifdef PERFDEBUG
3611  double patch_time = PIL_check_seconds_timer();
3612  std::cout << " patches found, time = " << patch_time - pwn_time << "\n";
3613 # endif
3614  CellsInfo cinfo = find_cells(tm_si, tm_si_topo, pinfo);
3615  if (dbg_level > 0) {
3616  std::cout << "Input is PWN\n";
3617  }
3618 # ifdef PERFDEBUG
3619  double cell_time = PIL_check_seconds_timer();
3620  std::cout << " cells found, time = " << cell_time - pwn_time << "\n";
3621 # endif
3622  finish_patch_cell_graph(tm_si, cinfo, pinfo, tm_si_topo, arena);
3623 # ifdef PERFDEBUG
3624  double finish_pc_time = PIL_check_seconds_timer();
3625  std::cout << " finished patch-cell graph, time = " << finish_pc_time - cell_time << "\n";
3626 # endif
3627  bool pc_ok = patch_cell_graph_ok(cinfo, pinfo);
3628  if (!pc_ok) {
3629  /* TODO: if bad input can lead to this, diagnose the problem. */
3630  std::cout << "Something funny about input or a bug in boolean\n";
3631  return IMesh(tm_in);
3632  }
3633  cinfo.init_windings(nshapes);
3634  int c_ambient = find_ambient_cell(tm_si, nullptr, tm_si_topo, pinfo, arena);
3635 # ifdef PERFDEBUG
3636  double amb_time = PIL_check_seconds_timer();
3637  std::cout << " ambient cell found, time = " << amb_time - finish_pc_time << "\n";
3638 # endif
3639  if (c_ambient == NO_INDEX) {
3640  /* TODO: find a way to propagate this error to user properly. */
3641  std::cout << "Could not find an ambient cell; input not valid?\n";
3642  return IMesh(tm_si);
3643  }
3644  propagate_windings_and_in_output_volume(pinfo, cinfo, c_ambient, op, nshapes, si_shape_fn);
3645 # ifdef PERFDEBUG
3646  double propagate_time = PIL_check_seconds_timer();
3647  std::cout << " windings propagated, time = " << propagate_time - amb_time << "\n";
3648 # endif
3649  tm_out = extract_from_in_output_volume_diffs(tm_si, pinfo, cinfo, arena);
3650 # ifdef PERFDEBUG
3651  double extract_time = PIL_check_seconds_timer();
3652  std::cout << " extracted, time = " << extract_time - propagate_time << "\n";
3653 # endif
3654  if (dbg_level > 0) {
3655  /* Check if output is PWN. */
3656  TriMeshTopology tm_out_topo(tm_out);
3657  if (!is_pwn(tm_out, tm_out_topo)) {
3658  std::cout << "OUTPUT IS NOT PWN!\n";
3659  }
3660  }
3661  }
3662  if (dbg_level > 1) {
3663  write_obj_mesh(tm_out, "boolean_tm_output");
3664  std::cout << "boolean tm output:\n" << tm_out;
3665  }
3666 # ifdef PERFDEBUG
3667  double end_time = PIL_check_seconds_timer();
3668  std::cout << " boolean_trimesh done, total time = " << end_time - start_time << "\n";
3669 # endif
3670  return tm_out;
3671 }
3672 
3673 static void dump_test_spec(IMesh &imesh)
3674 {
3675  std::cout << "test spec = " << imesh.vert_size() << " " << imesh.face_size() << "\n";
3676  for (const Vert *v : imesh.vertices()) {
3677  std::cout << v->co_exact[0] << " " << v->co_exact[1] << " " << v->co_exact[2] << " # "
3678  << v->co[0] << " " << v->co[1] << " " << v->co[2] << "\n";
3679  }
3680  for (const Face *f : imesh.faces()) {
3681  for (const Vert *fv : *f) {
3682  std::cout << imesh.lookup_vert(fv) << " ";
3683  }
3684  std::cout << "\n";
3685  }
3686 }
3687 
3688 IMesh boolean_mesh(IMesh &imesh,
3689  BoolOpType op,
3690  int nshapes,
3691  std::function<int(int)> shape_fn,
3692  bool use_self,
3693  bool hole_tolerant,
3694  IMesh *imesh_triangulated,
3695  IMeshArena *arena)
3696 {
3697  constexpr int dbg_level = 0;
3698  if (dbg_level > 0) {
3699  std::cout << "\nBOOLEAN_MESH\n"
3700  << nshapes << " operand" << (nshapes == 1 ? "" : "s")
3701  << " op=" << bool_optype_name(op) << "\n";
3702  if (dbg_level > 1) {
3703  write_obj_mesh(imesh, "boolean_mesh_in");
3704  std::cout << imesh;
3705  if (dbg_level > 2) {
3706  dump_test_spec(imesh);
3707  }
3708  }
3709  }
3710  IMesh *tm_in = imesh_triangulated;
3711  IMesh our_triangulation;
3712 # ifdef PERFDEBUG
3713  double start_time = PIL_check_seconds_timer();
3714  std::cout << "boolean_mesh, timing begins\n";
3715 # endif
3716  if (tm_in == nullptr) {
3717  our_triangulation = triangulate_polymesh(imesh, arena);
3718  tm_in = &our_triangulation;
3719  }
3720 # ifdef PERFDEBUG
3721  double tri_time = PIL_check_seconds_timer();
3722  std::cout << "triangulated, time = " << tri_time - start_time << "\n";
3723 # endif
3724  if (dbg_level > 1) {
3725  write_obj_mesh(*tm_in, "boolean_tm_in");
3726  }
3727  IMesh tm_out = boolean_trimesh(*tm_in, op, nshapes, shape_fn, use_self, hole_tolerant, arena);
3728 # ifdef PERFDEBUG
3729  double bool_tri_time = PIL_check_seconds_timer();
3730  std::cout << "boolean_trimesh done, time = " << bool_tri_time - tri_time << "\n";
3731 # endif
3732  if (dbg_level > 1) {
3733  std::cout << "bool_trimesh_output:\n" << tm_out;
3734  write_obj_mesh(tm_out, "bool_trimesh_output");
3735  }
3736  IMesh ans = polymesh_from_trimesh_with_dissolve(tm_out, imesh, arena);
3737 # ifdef PERFDEBUG
3738  double dissolve_time = PIL_check_seconds_timer();
3739  std::cout << "polymesh from dissolving, time = " << dissolve_time - bool_tri_time << "\n";
3740 # endif
3741  if (dbg_level > 0) {
3742  std::cout << "boolean_mesh output:\n" << ans;
3743  if (dbg_level > 2) {
3744  ans.populate_vert();
3745  dump_test_spec(ans);
3746  }
3747  }
3748 # ifdef PERFDEBUG
3749  double end_time = PIL_check_seconds_timer();
3750  std::cout << "boolean_mesh done, total time = " << end_time - start_time << "\n";
3751 # endif
3752  return ans;
3753 }
3754 
3755 } // namespace blender::meshintersect
3756 
3757 #endif // WITH_GMP
typedef float(TangentPoint)[2]
#define BLI_assert(a)
Definition: BLI_assert.h:46
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
void BLI_bvhtree_ray_cast_all(BVHTree *tree, const float co[3], const float dir[3], float radius, float hit_dist, BVHTree_RayCastCallback callback, void *userdata)
Definition: BLI_kdopbvh.c:2015
MINLINE double max_dd(double a, double b)
Math vector functions needed specifically for mesh intersect and boolean.
void closest_on_tri_to_point_v3(float r[3], const float p[3], const float v1[3], const float v2[3], const float v3[3])
Definition: math_geom.c:980
bool isect_ray_tri_epsilon_v3(const float ray_origin[3], const float ray_direction[3], const float v0[3], const float v1[3], const float v2[3], float *r_lambda, float r_uv[2], float epsilon)
Definition: math_geom.c:1735
MINLINE float len_squared_v3v3(const float a[3], const float b[3]) ATTR_WARN_UNUSED_RESULT
MINLINE void copy_v3fl_v3db(float r[3], const double a[3])
const char * BLI_getenv(const char *env) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT
Definition: path_util.c:1168
#define UNUSED_VARS_NDEBUG(...)
#define UNUSED(x)
#define ELEM(...)
typedef double(DMatrix)[4][4]
static uint8 component(Color32 c, uint i)
Definition: ColorBlock.cpp:108
_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 GLdouble r _GL_VOID_RET _GL_VOID GLfloat GLfloat r _GL_VOID_RET _GL_VOID GLint GLint r _GL_VOID_RET _GL_VOID GLshort GLshort r _GL_VOID_RET _GL_VOID GLdouble GLdouble r
_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
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
DBVT_INLINE bool Intersect(const btDbvtAabbMm &a, const btDbvtAabbMm &b)
Definition: btDbvt.h:621
void sort(btMatrix3x3 &U, btVector3 &sigma, btMatrix3x3 &V, int t)
Helper function of 3X3 SVD for sorting singular values.
SIMD_FORCE_INLINE btVector3 & operator[](int i)
Get a mutable reference to a row of the matrix as a vector.
Definition: btMatrix3x3.h:157
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
bool closest(btVector3 &v)
Patch()
Definition: subd/patch.h:14
SyclQueue * queue
void * tree
BLI_INLINE float fb(float length, float L)
IconTextureDrawCall normal
int count
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
ccl_device_inline float len_squared(const float3 a)
Definition: math_float3.h:423
struct Edge Edge
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
std::ostream & operator<<(std::ostream &stream, const AssetCatalogPath &path_to_append)
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, Size > normalize(const vec_base< T, Size > &v)
T abs(const T &a)
T length_squared(const vec_base< T, Size > &a)
std::ostream & operator<<(std::ostream &stream, const FatCo< T > &co)
Definition: delaunay_2d.cc:177
void parallel_for(IndexRange range, int64_t grain_size, const Function &function)
Definition: BLI_task.hh:51
Value parallel_reduce(IndexRange range, int64_t grain_size, const Value &identity, const Function &function, const Reduction &reduction)
Definition: BLI_task.hh:73
constexpr bool operator==(StringRef a, StringRef b)
vec_base< double, 3 > double3
int orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
uint64_t get_default_hash_2(const T1 &v1, const T2 &v2)
Definition: BLI_hash.hh:223
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
static const pxr::TfToken g("g", pxr::TfToken::Immortal)
static void copy(bNodeTree *dest_ntree, bNode *dest_node, const bNode *src_node)
#define hash
Definition: noise.c:153
unsigned __int64 uint64_t
Definition: stdint.h:90
float co[3]
Definition: bmesh_class.h:87
float origin[3]
Definition: BLI_kdopbvh.h:54
float direction[3]
Definition: BLI_kdopbvh.h:56
double PIL_check_seconds_timer(void)
Definition: time.c:64