Blender  V3.3
delaunay_2d.cc
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 #include <algorithm>
8 #include <atomic>
9 #include <fstream>
10 #include <iostream>
11 #include <sstream>
12 
13 #include "BLI_array.hh"
14 #include "BLI_linklist.h"
15 #include "BLI_math_boolean.hh"
16 #include "BLI_math_mpq.hh"
18 #include "BLI_set.hh"
19 #include "BLI_task.hh"
20 #include "BLI_vector.hh"
21 
22 #include "BLI_delaunay_2d.h"
23 
24 namespace blender::meshintersect {
25 
26 using namespace blender::math;
27 
28 /* Throughout this file, template argument T will be an
29  * arithmetic-like type, like float, double, or mpq_class. */
30 
31 template<typename T> T math_abs(const T v)
32 {
33  return (v < 0) ? -v : v;
34 }
35 
36 #ifdef WITH_GMP
37 template<> mpq_class math_abs<mpq_class>(const mpq_class v)
38 {
39  return abs(v);
40 }
41 #endif
42 
43 template<> double math_abs<double>(const double v)
44 {
45  return fabs(v);
46 }
47 
48 template<typename T> double math_to_double(const T UNUSED(v))
49 {
50  BLI_assert(false); /* Need implementation for other type. */
51  return 0.0;
52 }
53 
54 #ifdef WITH_GMP
55 template<> double math_to_double<mpq_class>(const mpq_class v)
56 {
57  return v.get_d();
58 }
59 #endif
60 
61 template<> double math_to_double<double>(const double v)
62 {
63  return v;
64 }
65 
77 template<typename Arith_t> struct CDTVert;
78 template<typename Arith_t> struct CDTEdge;
79 template<typename Arith_t> struct CDTFace;
80 
81 template<typename Arith_t> struct SymEdge {
85  SymEdge<Arith_t> *rot{nullptr};
87  CDTVert<Arith_t> *vert{nullptr};
89  CDTEdge<Arith_t> *edge{nullptr};
91  CDTFace<Arith_t> *face{nullptr};
92 
93  SymEdge() = default;
94 };
95 
99 template<typename T> inline SymEdge<T> *sym(const SymEdge<T> *se)
100 {
101  return se->next->rot;
102 }
103 
105 template<typename T> inline SymEdge<T> *prev(const SymEdge<T> *se)
106 {
107  return se->rot->next->rot;
108 }
109 
112 template<typename T> struct FatCo {
113  vec2<T> exact;
114  vec2<double> approx;
115  vec2<double> abs_approx;
116 
117  FatCo();
118 #ifdef WITH_GMP
119  FatCo(const vec2<mpq_class> &v);
120 #endif
121  FatCo(const vec2<double> &v);
122 };
123 
124 #ifdef WITH_GMP
125 template<> struct FatCo<mpq_class> {
126  vec2<mpq_class> exact;
127  vec2<double> approx;
128  vec2<double> abs_approx;
129 
130  FatCo()
131  : exact(vec2<mpq_class>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
132  {
133  }
134 
135  FatCo(const vec2<mpq_class> &v)
136  {
137  exact = v;
138  approx = vec2<double>(v.x.get_d(), v.y.get_d());
139  abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
140  }
141 
142  FatCo(const vec2<double> &v)
143  {
144  exact = vec2<mpq_class>(v.x, v.y);
145  approx = v;
146  abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
147  }
148 };
149 #endif
150 
151 template<> struct FatCo<double> {
152  vec2<double> exact;
153  vec2<double> approx;
154  vec2<double> abs_approx;
155 
156  FatCo() : exact(vec2<double>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
157  {
158  }
159 
160 #ifdef WITH_GMP
161  FatCo(const vec2<mpq_class> &v)
162  {
163  exact = vec2<double>(v.x.get_d(), v.y.get_d());
164  approx = exact;
165  abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
166  }
167 #endif
168 
169  FatCo(const vec2<double> &v)
170  {
171  exact = v;
172  approx = v;
173  abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
174  }
175 };
176 
177 template<typename T> std::ostream &operator<<(std::ostream &stream, const FatCo<T> &co)
178 {
179  stream << "(" << co.approx.x << ", " << co.approx.y << ")";
180  return stream;
181 }
182 
183 template<typename T> struct CDTVert {
187  SymEdge<T> *symedge{nullptr};
191  int index{-1};
193  int merge_to_index{-1};
195  int visit_index{0};
196 
197  CDTVert() = default;
198  explicit CDTVert(const vec2<T> &pt);
199 };
200 
201 template<typename Arith_t> struct CDTEdge {
208 
209  CDTEdge() = default;
210 };
211 
212 template<typename Arith_t> struct CDTFace {
214  SymEdge<Arith_t> *symedge{nullptr};
220  int visit_index{0};
222  bool deleted{false};
224  bool hole{false};
225 
226  CDTFace() = default;
227 };
228 
229 template<typename Arith_t> struct CDTArrangement {
230  /* The arrangement owns the memory pointed to by the pointers in these vectors.
231  * They are pointers instead of actual structures because these vectors may be resized and
232  * other elements refer to the elements by pointer. */
233 
241  CDTFace<Arith_t> *outer_face{nullptr};
242 
243  CDTArrangement() = default;
245 
248  void reserve(int verts_num, int edges_num, int faces_num);
249 
254  CDTVert<Arith_t> *add_vert(const vec2<Arith_t> &pt);
255 
264  CDTFace<Arith_t> *fleft,
265  CDTFace<Arith_t> *fright);
266 
272 
275 
283 
289 
295 
302 
308  {
309  CDTVert<Arith_t> *v = this->verts[i];
310  if (v->merge_to_index != -1) {
311  v = this->verts[v->merge_to_index];
312  }
313  return v;
314  }
315 };
316 
317 template<typename T> class CDT_state {
318  public:
332  bool need_ids;
333 
334  explicit CDT_state(
335  int input_verts_num, int input_edges_num, int input_faces_num, T epsilon, bool need_ids);
336 };
337 
339 {
340  for (int i : this->verts.index_range()) {
341  CDTVert<T> *v = this->verts[i];
342  v->input_ids.clear();
343  delete v;
344  this->verts[i] = nullptr;
345  }
346  for (int i : this->edges.index_range()) {
347  CDTEdge<T> *e = this->edges[i];
348  e->input_ids.clear();
349  delete e;
350  this->edges[i] = nullptr;
351  }
352  for (int i : this->faces.index_range()) {
353  CDTFace<T> *f = this->faces[i];
354  f->input_ids.clear();
355  delete f;
356  this->faces[i] = nullptr;
357  }
358 }
359 
360 #define DEBUG_CDT
361 #ifdef DEBUG_CDT
362 /* Some functions to aid in debugging. */
363 template<typename T> std::string vertname(const CDTVert<T> *v)
364 {
365  std::stringstream ss;
366  ss << "[" << v->index << "]";
367  return ss.str();
368 }
369 
370 /* Abbreviated pointer value is easier to read in dumps. */
371 static std::string trunc_ptr(const void *p)
372 {
373  constexpr int TRUNC_PTR_MASK = 0xFFFF;
374  std::stringstream ss;
375  ss << std::hex << (POINTER_AS_INT(p) & TRUNC_PTR_MASK);
376  return ss.str();
377 }
378 
379 template<typename T> std::string sename(const SymEdge<T> *se)
380 {
381  std::stringstream ss;
382  ss << "{" << trunc_ptr(se) << "}";
383  return ss.str();
384 }
385 
386 template<typename T> std::ostream &operator<<(std::ostream &os, const SymEdge<T> &se)
387 {
388  if (se.next) {
389  os << vertname(se.vert) << "(" << se.vert->co << "->" << se.next->vert->co << ")"
390  << vertname(se.next->vert);
391  }
392  else {
393  os << vertname(se.vert) << "(" << se.vert->co << "->NULL)";
394  }
395  return os;
396 }
397 
398 template<typename T> std::ostream &operator<<(std::ostream &os, const SymEdge<T> *se)
399 {
400  os << *se;
401  return os;
402 }
403 
404 template<typename T> std::string short_se_dump(const SymEdge<T> *se)
405 {
406  if (se == nullptr) {
407  return std::string("NULL");
408  }
409  return vertname(se->vert) +
410  (se->next == nullptr ? std::string("[NULL]") : vertname(se->next->vert));
411 }
412 
413 template<typename T> std::ostream &operator<<(std::ostream &os, const CDT_state<T> &cdt_state)
414 {
415  os << "\nCDT\n\nVERTS\n";
416  for (const CDTVert<T> *v : cdt_state.cdt.verts) {
417  os << vertname(v) << " " << trunc_ptr(v) << ": " << v->co
418  << " symedge=" << trunc_ptr(v->symedge);
419  if (v->merge_to_index == -1) {
420  os << "\n";
421  }
422  else {
423  os << " merge to " << vertname(cdt_state.cdt.verts[v->merge_to_index]) << "\n";
424  }
425  const SymEdge<T> *se = v->symedge;
426  int count = 0;
427  constexpr int print_count_limit = 25;
428  if (se) {
429  os << " edges out:\n";
430  do {
431  if (se->next == nullptr) {
432  os << " [NULL] next/rot symedge, se=" << trunc_ptr(se) << "\n";
433  break;
434  }
435  if (se->next->next == nullptr) {
436  os << " [NULL] next-next/rot symedge, se=" << trunc_ptr(se) << "\n";
437  break;
438  }
439  const CDTVert<T> *vother = sym(se)->vert;
440  os << " " << vertname(vother) << "(e=" << trunc_ptr(se->edge)
441  << ", se=" << trunc_ptr(se) << ")\n";
442  se = se->rot;
443  count++;
444  } while (se != v->symedge && count < print_count_limit);
445  os << "\n";
446  }
447  }
448  os << "\nEDGES\n";
449  for (const CDTEdge<T> *e : cdt_state.cdt.edges) {
450  if (e->symedges[0].next == nullptr) {
451  continue;
452  }
453  os << trunc_ptr(&e) << ":\n";
454  for (int i = 0; i < 2; ++i) {
455  const SymEdge<T> *se = &e->symedges[i];
456  os << " se[" << i << "] @" << trunc_ptr(se) << " next=" << trunc_ptr(se->next)
457  << ", rot=" << trunc_ptr(se->rot) << ", vert=" << trunc_ptr(se->vert) << " "
458  << vertname(se->vert) << " " << se->vert->co << ", edge=" << trunc_ptr(se->edge)
459  << ", face=" << trunc_ptr(se->face) << "\n";
460  }
461  }
462  os << "\nFACES\n";
463  os << "outer_face=" << trunc_ptr(cdt_state.cdt.outer_face) << "\n";
464  /* Only after prepare_output do faces have non-null symedges. */
465  if (cdt_state.cdt.outer_face->symedge != nullptr) {
466  for (const CDTFace<T> *f : cdt_state.cdt.faces) {
467  if (!f->deleted) {
468  os << trunc_ptr(f) << " symedge=" << trunc_ptr(f->symedge) << "\n";
469  }
470  }
471  }
472  return os;
473 }
474 
475 template<typename T> void cdt_draw(const std::string &label, const CDTArrangement<T> &cdt)
476 {
477  static bool append = false; /* Will be set to true after first call. */
478 
479 /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library.
480  * This is just for developer debugging anyway, and should never be called in production Blender.
481  */
482 # ifdef _WIN32
483  const char *drawfile = "./cdt_debug_draw.html";
484 # else
485  const char *drawfile = "/tmp/cdt_debug_draw.html";
486 # endif
487  constexpr int max_draw_width = 1800;
488  constexpr int max_draw_height = 1600;
489  constexpr double margin_expand = 0.05;
490  constexpr int thin_line = 1;
491  constexpr int thick_line = 4;
492  constexpr int vert_radius = 3;
493  constexpr bool draw_vert_labels = true;
494  constexpr bool draw_edge_labels = false;
495  constexpr bool draw_face_labels = false;
496 
497  if (cdt.verts.size() == 0) {
498  return;
499  }
500  vec2<double> vmin(DBL_MAX, DBL_MAX);
501  vec2<double> vmax(-DBL_MAX, -DBL_MAX);
502  for (const CDTVert<T> *v : cdt.verts) {
503  for (int i = 0; i < 2; ++i) {
504  double dvi = v->co.approx[i];
505  if (dvi < vmin[i]) {
506  vmin[i] = dvi;
507  }
508  if (dvi > vmax[i]) {
509  vmax[i] = dvi;
510  }
511  }
512  }
513  double draw_margin = ((vmax.x - vmin.x) + (vmax.y - vmin.y)) * margin_expand;
514  double minx = vmin.x - draw_margin;
515  double maxx = vmax.x + draw_margin;
516  double miny = vmin.y - draw_margin;
517  double maxy = vmax.y + draw_margin;
518 
519  double width = maxx - minx;
520  double height = maxy - miny;
521  double aspect = height / width;
522  int view_width = max_draw_width;
523  int view_height = static_cast<int>(view_width * aspect);
524  if (view_height > max_draw_height) {
525  view_height = max_draw_height;
526  view_width = static_cast<int>(view_height / aspect);
527  }
528  double scale = view_width / width;
529 
530 # define SX(x) (((x)-minx) * scale)
531 # define SY(y) ((maxy - (y)) * scale)
532 
533  std::ofstream f;
534  if (append) {
535  f.open(drawfile, std::ios_base::app);
536  }
537  else {
538  f.open(drawfile);
539  }
540  if (!f) {
541  std::cout << "Could not open file " << drawfile << "\n";
542  return;
543  }
544 
545  f << "<div>" << label << "</div>\n<div>\n"
546  << "<svg version=\"1.1\" "
547  "xmlns=\"http://www.w3.org/2000/svg\" "
548  "xmlns:xlink=\"http://www.w3.org/1999/xlink\" "
549  "xml:space=\"preserve\"\n"
550  << "width=\"" << view_width << "\" height=\"" << view_height << "\">n";
551 
552  for (const CDTEdge<T> *e : cdt.edges) {
553  if (e->symedges[0].next == nullptr) {
554  continue;
555  }
556  const CDTVert<T> *u = e->symedges[0].vert;
557  const CDTVert<T> *v = e->symedges[1].vert;
558  const vec2<double> &uco = u->co.approx;
559  const vec2<double> &vco = v->co.approx;
560  int strokew = e->input_ids.size() == 0 ? thin_line : thick_line;
561  f << R"(<line fill="none" stroke="black" stroke-width=")" << strokew << "\" x1=\""
562  << SX(uco[0]) << "\" y1=\"" << SY(uco[1]) << "\" x2=\"" << SX(vco[0]) << "\" y2=\""
563  << SY(vco[1]) << "\">\n";
564  f << " <title>" << vertname(u) << vertname(v) << "</title>\n";
565  f << "</line>\n";
566  if (draw_edge_labels) {
567  f << "<text x=\"" << SX((uco[0] + vco[0]) / 2) << "\" y=\"" << SY((uco[1] + vco[1]) / 2)
568  << R"(" font-size="small">)";
569  f << vertname(u) << vertname(v) << sename(&e->symedges[0]) << sename(&e->symedges[1])
570  << "</text>\n";
571  }
572  }
573 
574  int i = 0;
575  for (const CDTVert<T> *v : cdt.verts) {
576  f << R"(<circle fill="black" cx=")" << SX(v->co.approx[0]) << "\" cy=\"" << SY(v->co.approx[1])
577  << "\" r=\"" << vert_radius << "\">\n";
578  f << " <title>[" << i << "]" << v->co.approx << "</title>\n";
579  f << "</circle>\n";
580  if (draw_vert_labels) {
581  f << "<text x=\"" << SX(v->co.approx[0]) + vert_radius << "\" y=\""
582  << SY(v->co.approx[1]) - vert_radius << R"(" font-size="small">[)" << i << "]</text>\n";
583  }
584  ++i;
585  }
586 
587  if (draw_face_labels) {
588  for (const CDTFace<T> *face : cdt.faces) {
589  if (!face->deleted) {
590  /* Since may not have prepared output yet, need a slow find of a SymEdge for this face. */
591  const SymEdge<T> *se_face_start = nullptr;
592  for (const CDTEdge<T> *e : cdt.edges) {
593  if (e->symedges[0].face == face) {
594  se_face_start = &e->symedges[0];
595  break;
596  }
597  if (e->symedges[1].face == face) {
598  se_face_start = &e->symedges[1];
599  }
600  }
601  if (se_face_start != nullptr) {
602  /* Find center of face. */
603  int face_nverts = 0;
604  vec2<double> cen(0.0, 0.0);
605  if (face == cdt.outer_face) {
606  cen.x = minx;
607  cen.y = miny;
608  }
609  else {
610  const SymEdge<T> *se = se_face_start;
611  do {
612  if (se->face == face) {
613  cen = cen + se->vert->co.approx;
614  face_nverts++;
615  }
616  } while ((se = se->next) != se_face_start);
617  if (face_nverts > 0) {
618  cen = cen / double(face_nverts);
619  }
620  }
621  f << "<text x=\"" << SX(cen[0]) << "\" y=\"" << SY(cen[1]) << "\""
622  << " font-size=\"small\">[";
623  f << trunc_ptr(face);
624  if (face->input_ids.size() > 0) {
625  for (int id : face->input_ids) {
626  f << " " << id;
627  }
628  }
629  f << "]</text>\n";
630  }
631  }
632  }
633  }
634 
635  append = true;
636 # undef SX
637 # undef SY
638 }
639 #endif
640 
645 template<typename T>
646 static int filtered_orient2d(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c);
647 
648 #ifdef WITH_GMP
649 template<>
650 int filtered_orient2d<mpq_class>(const FatCo<mpq_class> &a,
651  const FatCo<mpq_class> &b,
652  const FatCo<mpq_class> &c)
653 {
654  double det = (a.approx[0] - c.approx[0]) * (b.approx[1] - c.approx[1]) -
655  (a.approx[1] - c.approx[1]) * (b.approx[0] - c.approx[0]);
656  double supremum = (a.abs_approx[0] + c.abs_approx[0]) * (b.abs_approx[1] + c.abs_approx[1]) +
657  (a.abs_approx[1] + c.abs_approx[1]) * (b.abs_approx[0] + c.abs_approx[0]);
658  constexpr double index_orient2d = 6;
659  double err_bound = supremum * index_orient2d * DBL_EPSILON;
660  if (fabs(det) > err_bound) {
661  return det > 0 ? 1 : -1;
662  }
663  return orient2d(a.exact, b.exact, c.exact);
664 }
665 #endif
666 
667 template<>
668 int filtered_orient2d<double>(const FatCo<double> &a,
669  const FatCo<double> &b,
670  const FatCo<double> &c)
671 {
672  return orient2d(a.approx, b.approx, c.approx);
673 }
674 
678 template<typename T>
679 static int filtered_incircle(const FatCo<T> &a,
680  const FatCo<T> &b,
681  const FatCo<T> &c,
682  const FatCo<T> &d);
683 
684 #ifdef WITH_GMP
685 template<>
686 int filtered_incircle<mpq_class>(const FatCo<mpq_class> &a,
687  const FatCo<mpq_class> &b,
688  const FatCo<mpq_class> &c,
689  const FatCo<mpq_class> &d)
690 {
691  double adx = a.approx[0] - d.approx[0];
692  double bdx = b.approx[0] - d.approx[0];
693  double cdx = c.approx[0] - d.approx[0];
694  double ady = a.approx[1] - d.approx[1];
695  double bdy = b.approx[1] - d.approx[1];
696  double cdy = c.approx[1] - d.approx[1];
697  double bdxcdy = bdx * cdy;
698  double cdxbdy = cdx * bdy;
699  double alift = adx * adx + ady * ady;
700  double cdxady = cdx * ady;
701  double adxcdy = adx * cdy;
702  double blift = bdx * bdx + bdy * bdy;
703  double adxbdy = adx * bdy;
704  double bdxady = bdx * ady;
705  double clift = cdx * cdx + cdy * cdy;
706  double det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);
707 
708  double sup_adx = a.abs_approx[0] + d.abs_approx[0]; /* index 2. */
709  double sup_bdx = b.abs_approx[0] + d.abs_approx[0];
710  double sup_cdx = c.abs_approx[0] + d.abs_approx[0];
711  double sup_ady = a.abs_approx[1] + d.abs_approx[1];
712  double sup_bdy = b.abs_approx[1] + d.abs_approx[1];
713  double sup_cdy = c.abs_approx[1] + d.abs_approx[1];
714  double sup_bdxcdy = sup_bdx * sup_cdy; /* index 5. */
715  double sup_cdxbdy = sup_cdx * sup_bdy;
716  double sup_alift = sup_adx * sup_adx + sup_ady * sup_ady; /* index 6. */
717  double sup_cdxady = sup_cdx * sup_ady;
718  double sup_adxcdy = sup_adx * sup_cdy;
719  double sup_blift = sup_bdx * sup_bdx + sup_bdy * sup_bdy;
720  double sup_adxbdy = sup_adx * sup_bdy;
721  double sup_bdxady = sup_bdx * sup_ady;
722  double sup_clift = sup_cdx * sup_cdx + sup_cdy * sup_cdy;
723  double sup_det = sup_alift * (sup_bdxcdy + sup_cdxbdy) + sup_blift * (sup_cdxady + sup_adxcdy) +
724  sup_clift * (sup_adxbdy + sup_bdxady);
725  int index = 14;
726  double err_bound = sup_det * index * DBL_EPSILON;
727  if (fabs(det) > err_bound) {
728  return det < 0.0 ? -1 : (det > 0.0 ? 1 : 0);
729  }
730  return incircle(a.exact, b.exact, c.exact, d.exact);
731 }
732 #endif
733 
734 template<>
735 int filtered_incircle<double>(const FatCo<double> &a,
736  const FatCo<double> &b,
737  const FatCo<double> &c,
738  const FatCo<double> &d)
739 {
740  return incircle(a.approx, b.approx, c.approx, d.approx);
741 }
742 
749 template<typename T> static bool in_line(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c);
750 
751 #ifdef WITH_GMP
752 template<>
753 bool in_line<mpq_class>(const FatCo<mpq_class> &a,
754  const FatCo<mpq_class> &b,
755  const FatCo<mpq_class> &c)
756 {
757  vec2<double> ab = b.approx - a.approx;
758  vec2<double> bc = c.approx - b.approx;
759  vec2<double> ac = c.approx - a.approx;
760  vec2<double> supremum_ab = b.abs_approx + a.abs_approx;
761  vec2<double> supremum_bc = c.abs_approx + b.abs_approx;
762  vec2<double> supremum_ac = c.abs_approx + a.abs_approx;
763  double dot_ab_ac = ab.x * ac.x + ab.y * ac.y;
764  double supremum_dot_ab_ac = supremum_ab.x * supremum_ac.x + supremum_ab.y * supremum_ac.y;
765  constexpr double index = 6;
766  double err_bound = supremum_dot_ab_ac * index * DBL_EPSILON;
767  if (dot_ab_ac < -err_bound) {
768  return false;
769  }
770  double dot_bc_ac = bc.x * ac.x + bc.y * ac.y;
771  double supremum_dot_bc_ac = supremum_bc.x * supremum_ac.x + supremum_bc.y * supremum_ac.y;
772  err_bound = supremum_dot_bc_ac * index * DBL_EPSILON;
773  if (dot_bc_ac < -err_bound) {
774  return false;
775  }
776  vec2<mpq_class> exact_ab = b.exact - a.exact;
777  vec2<mpq_class> exact_ac = c.exact - a.exact;
778  if (dot(exact_ab, exact_ac) < 0) {
779  return false;
780  }
781  vec2<mpq_class> exact_bc = c.exact - b.exact;
782  return dot(exact_bc, exact_ac) >= 0;
783 }
784 #endif
785 
786 template<>
787 bool in_line<double>(const FatCo<double> &a, const FatCo<double> &b, const FatCo<double> &c)
788 {
789  vec2<double> ab = b.approx - a.approx;
790  vec2<double> ac = c.approx - a.approx;
791  if (dot(ab, ac) < 0) {
792  return false;
793  }
794  vec2<double> bc = c.approx - b.approx;
795  return dot(bc, ac) >= 0;
796 }
797 
798 template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
799 {
800  this->co.exact = pt;
801  this->co.approx = pt;
802  this->co.abs_approx = pt; /* Not used, so doesn't matter. */
803  this->symedge = nullptr;
804  this->index = -1;
805  this->merge_to_index = -1;
806  this->visit_index = 0;
807 }
808 
809 #ifdef WITH_GMP
810 template<> CDTVert<mpq_class>::CDTVert(const vec2<mpq_class> &pt)
811 {
812  this->co.exact = pt;
813  this->co.approx = double2(pt.x.get_d(), pt.y.get_d());
814  this->co.abs_approx = double2(fabs(this->co.approx.x), fabs(this->co.approx.y));
815  this->symedge = nullptr;
816  this->index = -1;
817  this->merge_to_index = -1;
818  this->visit_index = 0;
819 }
820 #endif
821 
822 template<typename T> CDTVert<T> *CDTArrangement<T>::add_vert(const vec2<T> &pt)
823 {
824  CDTVert<T> *v = new CDTVert<T>(pt);
825  int index = this->verts.append_and_get_index(v);
826  v->index = index;
827  return v;
828 }
829 
830 template<typename T>
831 CDTEdge<T> *CDTArrangement<T>::add_edge(CDTVert<T> *v1,
832  CDTVert<T> *v2,
833  CDTFace<T> *fleft,
834  CDTFace<T> *fright)
835 {
836  CDTEdge<T> *e = new CDTEdge<T>();
837  this->edges.append(e);
838  SymEdge<T> *se = &e->symedges[0];
839  SymEdge<T> *sesym = &e->symedges[1];
840  se->edge = sesym->edge = e;
841  se->face = fleft;
842  sesym->face = fright;
843  se->vert = v1;
844  if (v1->symedge == nullptr) {
845  v1->symedge = se;
846  }
847  sesym->vert = v2;
848  if (v2->symedge == nullptr) {
849  v2->symedge = sesym;
850  }
851  se->next = sesym->next = se->rot = sesym->rot = nullptr;
852  return e;
853 }
854 
855 template<typename T> CDTFace<T> *CDTArrangement<T>::add_face()
856 {
857  CDTFace<T> *f = new CDTFace<T>();
858  this->faces.append(f);
859  return f;
860 }
861 
862 template<typename T> void CDTArrangement<T>::reserve(int verts_num, int edges_num, int faces_num)
863 {
864  /* These reserves are just guesses; OK if they aren't exactly right since vectors will resize. */
865  this->verts.reserve(2 * verts_num);
866  this->edges.reserve(3 * verts_num + 2 * edges_num + 3 * 2 * faces_num);
867  this->faces.reserve(2 * verts_num + 2 * edges_num + 2 * faces_num);
868 }
869 
870 template<typename T>
872  int input_verts_num, int input_edges_num, int input_faces_num, T epsilon, bool need_ids)
873 {
874  this->input_vert_num = input_verts_num;
875  this->cdt.reserve(input_verts_num, input_edges_num, input_faces_num);
876  this->cdt.outer_face = this->cdt.add_face();
877  this->epsilon = epsilon;
878  this->need_ids = need_ids;
879  this->visit_count = 0;
880 }
881 
882 /* Is any id in (range_start, range_start+1, ... , range_end) in id_list? */
883 static bool id_range_in_list(const blender::Set<int> &id_list, int range_start, int range_end)
884 {
885  for (int id : id_list) {
886  if (id >= range_start && id <= range_end) {
887  return true;
888  }
889  }
890  return false;
891 }
892 
893 static void add_to_input_ids(blender::Set<int> &dst, int input_id)
894 {
895  dst.add(input_id);
896 }
897 
899 {
900  for (int value : src) {
901  dst.add(value);
902  }
903 }
904 
905 template<typename T> inline bool is_border_edge(const CDTEdge<T> *e, const CDT_state<T> *cdt)
906 {
907  return e->symedges[0].face == cdt->outer_face || e->symedges[1].face == cdt->outer_face;
908 }
909 
910 template<typename T> inline bool is_constrained_edge(const CDTEdge<T> *e)
911 {
912  return e->input_ids.size() > 0;
913 }
914 
915 template<typename T> inline bool is_deleted_edge(const CDTEdge<T> *e)
916 {
917  return e->symedges[0].next == nullptr;
918 }
919 
920 template<typename T> inline bool is_original_vert(const CDTVert<T> *v, CDT_state<T> *cdt)
921 {
922  return (v->index < cdt->input_vert_num);
923 }
924 
928 template<typename T>
929 SymEdge<T> *find_symedge_between_verts(const CDTVert<T> *v1, const CDTVert<T> *v2)
930 {
931  SymEdge<T> *t = v1->symedge;
932  SymEdge<T> *tstart = t;
933  do {
934  if (t->next->vert == v2) {
935  return t;
936  }
937  } while ((t = t->rot) != tstart);
938  return nullptr;
939 }
940 
944 template<typename T> SymEdge<T> *find_symedge_with_face(const CDTVert<T> *v, const CDTFace<T> *f)
945 {
946  SymEdge<T> *t = v->symedge;
947  SymEdge<T> *tstart = t;
948  do {
949  if (t->face == f) {
950  return t;
951  }
952  } while ((t = t->rot) != tstart);
953  return nullptr;
954 }
955 
959 template<typename T> inline bool exists_edge(const CDTVert<T> *v1, const CDTVert<T> *v2)
960 {
961  return find_symedge_between_verts(v1, v2) != nullptr;
962 }
963 
967 template<typename T> bool vert_touches_face(const CDTVert<T> *v, const CDTFace<T> *f)
968 {
969  SymEdge<T> *se = v->symedge;
970  do {
971  if (se->face == f) {
972  return true;
973  }
974  } while ((se = se->rot) != v->symedge);
975  return false;
976 }
977 
986 template<typename T> CDTEdge<T> *CDTArrangement<T>::add_diagonal(SymEdge<T> *s1, SymEdge<T> *s2)
987 {
988  CDTFace<T> *fold = s1->face;
989  CDTFace<T> *fnew = this->add_face();
990  SymEdge<T> *s1prev = prev(s1);
991  SymEdge<T> *s1prevsym = sym(s1prev);
992  SymEdge<T> *s2prev = prev(s2);
993  SymEdge<T> *s2prevsym = sym(s2prev);
994  CDTEdge<T> *ediag = this->add_edge(s1->vert, s2->vert, fnew, fold);
995  SymEdge<T> *sdiag = &ediag->symedges[0];
996  SymEdge<T> *sdiagsym = &ediag->symedges[1];
997  sdiag->next = s2;
998  sdiagsym->next = s1;
999  s2prev->next = sdiagsym;
1000  s1prev->next = sdiag;
1001  s1->rot = sdiag;
1002  sdiag->rot = s1prevsym;
1003  s2->rot = sdiagsym;
1004  sdiagsym->rot = s2prevsym;
1005  for (SymEdge<T> *se = s2; se != sdiag; se = se->next) {
1006  se->face = fnew;
1007  }
1008  add_list_to_input_ids(fnew->input_ids, fold->input_ids);
1009  return ediag;
1010 }
1011 
1012 template<typename T>
1013 CDTEdge<T> *CDTArrangement<T>::add_vert_to_symedge_edge(CDTVert<T> *v, SymEdge<T> *se)
1014 {
1015  SymEdge<T> *se_rot = se->rot;
1016  SymEdge<T> *se_rotsym = sym(se_rot);
1017  /* TODO: check: I think last arg in next should be sym(se)->face. */
1018  CDTEdge<T> *e = this->add_edge(v, se->vert, se->face, se->face);
1019  SymEdge<T> *new_se = &e->symedges[0];
1020  SymEdge<T> *new_se_sym = &e->symedges[1];
1021  new_se->next = se;
1022  new_se_sym->next = new_se;
1023  new_se->rot = new_se;
1024  new_se_sym->rot = se_rot;
1025  se->rot = new_se_sym;
1026  se_rotsym->next = new_se_sym;
1027  return e;
1028 }
1029 
1035 template<typename T>
1036 CDTEdge<T> *CDTArrangement<T>::connect_separate_parts(SymEdge<T> *se1, SymEdge<T> *se2)
1037 {
1038  BLI_assert(se1->face == this->outer_face && se2->face == this->outer_face);
1039  SymEdge<T> *se1_rot = se1->rot;
1040  SymEdge<T> *se1_rotsym = sym(se1_rot);
1041  SymEdge<T> *se2_rot = se2->rot;
1042  SymEdge<T> *se2_rotsym = sym(se2_rot);
1043  CDTEdge<T> *e = this->add_edge(se1->vert, se2->vert, this->outer_face, this->outer_face);
1044  SymEdge<T> *new_se = &e->symedges[0];
1045  SymEdge<T> *new_se_sym = &e->symedges[1];
1046  new_se->next = se2;
1047  new_se_sym->next = se1;
1048  new_se->rot = se1_rot;
1049  new_se_sym->rot = se2_rot;
1050  se1->rot = new_se;
1051  se2->rot = new_se_sym;
1052  se1_rotsym->next = new_se;
1053  se2_rotsym->next = new_se_sym;
1054  return e;
1055 }
1056 
1062 template<typename T> CDTEdge<T> *CDTArrangement<T>::split_edge(SymEdge<T> *se, T lambda)
1063 {
1064  /* Split e at lambda. */
1065  const vec2<T> *a = &se->vert->co.exact;
1066  const vec2<T> *b = &se->next->vert->co.exact;
1067  SymEdge<T> *sesym = sym(se);
1068  SymEdge<T> *sesymprev = prev(sesym);
1069  SymEdge<T> *sesymprevsym = sym(sesymprev);
1070  SymEdge<T> *senext = se->next;
1071  CDTVert<T> *v = this->add_vert(interpolate(*a, *b, lambda));
1072  CDTEdge<T> *e = this->add_edge(v, se->next->vert, se->face, sesym->face);
1073  sesym->vert = v;
1074  SymEdge<T> *newse = &e->symedges[0];
1075  SymEdge<T> *newsesym = &e->symedges[1];
1076  se->next = newse;
1077  newsesym->next = sesym;
1078  newse->next = senext;
1079  newse->rot = sesym;
1080  sesym->rot = newse;
1081  senext->rot = newsesym;
1082  newsesym->rot = sesymprevsym;
1083  sesymprev->next = newsesym;
1084  if (newsesym->vert->symedge == sesym) {
1085  newsesym->vert->symedge = newsesym;
1086  }
1087  add_list_to_input_ids(e->input_ids, se->edge->input_ids);
1088  return e;
1089 }
1090 
1112 template<typename T> void CDTArrangement<T>::delete_edge(SymEdge<T> *se)
1113 {
1114  SymEdge<T> *sesym = sym(se);
1115  CDTVert<T> *v1 = se->vert;
1116  CDTVert<T> *v2 = sesym->vert;
1117  CDTFace<T> *aface = se->face;
1118  CDTFace<T> *bface = sesym->face;
1119  SymEdge<T> *f = se->next;
1120  SymEdge<T> *h = prev(se);
1121  SymEdge<T> *i = sesym->next;
1122  SymEdge<T> *j = prev(sesym);
1123  SymEdge<T> *jsym = sym(j);
1124  SymEdge<T> *hsym = sym(h);
1125  bool v1_isolated = (i == se);
1126  bool v2_isolated = (f == sesym);
1127 
1128  if (!v1_isolated) {
1129  h->next = i;
1130  i->rot = hsym;
1131  }
1132  if (!v2_isolated) {
1133  j->next = f;
1134  f->rot = jsym;
1135  }
1136  if (!v1_isolated && !v2_isolated && aface != bface) {
1137  for (SymEdge<T> *k = i; k != f; k = k->next) {
1138  k->face = aface;
1139  }
1140  }
1141 
1142  /* If e was representative symedge for v1 or v2, fix that. */
1143  if (v1_isolated) {
1144  v1->symedge = nullptr;
1145  }
1146  else if (v1->symedge == se) {
1147  v1->symedge = i;
1148  }
1149  if (v2_isolated) {
1150  v2->symedge = nullptr;
1151  }
1152  else if (v2->symedge == sesym) {
1153  v2->symedge = f;
1154  }
1155 
1156  /* Mark SymEdge as deleted by setting all its pointers to NULL. */
1157  se->next = se->rot = nullptr;
1158  sesym->next = sesym->rot = nullptr;
1159  if (!v1_isolated && !v2_isolated && aface != bface) {
1160  bface->deleted = true;
1161  if (this->outer_face == bface) {
1162  this->outer_face = aface;
1163  }
1164  }
1165 }
1166 
1167 template<typename T> class SiteInfo {
1168  public:
1169  CDTVert<T> *v;
1171 };
1172 
1176 template<typename T> bool site_lexicographic_sort(const SiteInfo<T> &a, const SiteInfo<T> &b)
1177 {
1178  const vec2<T> &co_a = a.v->co.exact;
1179  const vec2<T> &co_b = b.v->co.exact;
1180  if (co_a[0] < co_b[0]) {
1181  return true;
1182  }
1183  if (co_a[0] > co_b[0]) {
1184  return false;
1185  }
1186  if (co_a[1] < co_b[1]) {
1187  return true;
1188  }
1189  if (co_a[1] > co_b[1]) {
1190  return false;
1191  }
1192  return a.orig_index < b.orig_index;
1193 }
1194 
1200 template<typename T> void find_site_merges(Array<SiteInfo<T>> &sites)
1201 {
1202  int n = sites.size();
1203  for (int i = 0; i < n - 1; ++i) {
1204  int j = i + 1;
1205  while (j < n && sites[j].v->co.exact == sites[i].v->co.exact) {
1206  sites[j].v->merge_to_index = sites[i].orig_index;
1207  ++j;
1208  }
1209  if (j - i > 1) {
1210  i = j - 1; /* j-1 because loop head will add another 1. */
1211  }
1212  }
1213 }
1214 
1215 template<typename T> inline bool vert_left_of_symedge(CDTVert<T> *v, SymEdge<T> *se)
1216 {
1217  return filtered_orient2d(v->co, se->vert->co, se->next->vert->co) > 0;
1218 }
1219 
1220 template<typename T> inline bool vert_right_of_symedge(CDTVert<T> *v, SymEdge<T> *se)
1221 {
1222  return filtered_orient2d(v->co, se->next->vert->co, se->vert->co) > 0;
1223 }
1224 
1225 /* Is se above basel? */
1226 template<typename T>
1227 inline bool dc_tri_valid(SymEdge<T> *se, SymEdge<T> *basel, SymEdge<T> *basel_sym)
1228 {
1229  return filtered_orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0;
1230 }
1231 
1238 template<typename T>
1239 void dc_tri(CDTArrangement<T> *cdt,
1240  Array<SiteInfo<T>> &sites,
1241  int start,
1242  int end,
1243  SymEdge<T> **r_le,
1244  SymEdge<T> **r_re)
1245 {
1246  constexpr int dbg_level = 0;
1247  if (dbg_level > 0) {
1248  std::cout << "DC_TRI start=" << start << " end=" << end << "\n";
1249  }
1250  int n = end - start;
1251  if (n <= 1) {
1252  *r_le = nullptr;
1253  *r_re = nullptr;
1254  return;
1255  }
1256 
1257  /* Base case: if n <= 3, triangulate directly. */
1258  if (n <= 3) {
1259  CDTVert<T> *v1 = sites[start].v;
1260  CDTVert<T> *v2 = sites[start + 1].v;
1261  CDTEdge<T> *ea = cdt->add_edge(v1, v2, cdt->outer_face, cdt->outer_face);
1262  ea->symedges[0].next = &ea->symedges[1];
1263  ea->symedges[1].next = &ea->symedges[0];
1264  ea->symedges[0].rot = &ea->symedges[0];
1265  ea->symedges[1].rot = &ea->symedges[1];
1266  if (n == 2) {
1267  *r_le = &ea->symedges[0];
1268  *r_re = &ea->symedges[1];
1269  return;
1270  }
1271  CDTVert<T> *v3 = sites[start + 2].v;
1272  CDTEdge<T> *eb = cdt->add_vert_to_symedge_edge(v3, &ea->symedges[1]);
1273  int orient = filtered_orient2d(v1->co, v2->co, v3->co);
1274  if (orient > 0) {
1275  cdt->add_diagonal(&eb->symedges[0], &ea->symedges[0]);
1276  *r_le = &ea->symedges[0];
1277  *r_re = &eb->symedges[0];
1278  }
1279  else if (orient < 0) {
1280  cdt->add_diagonal(&ea->symedges[0], &eb->symedges[0]);
1281  *r_le = ea->symedges[0].rot;
1282  *r_re = eb->symedges[0].rot;
1283  }
1284  else {
1285  /* Collinear points. Just return a line. */
1286  *r_le = &ea->symedges[0];
1287  *r_re = &eb->symedges[0];
1288  }
1289  return;
1290  }
1291  /* Recursive case. Do left (L) and right (R) halves separately, then join. */
1292  int n2 = n / 2;
1293  BLI_assert(n2 >= 2 && end - (start + n2) >= 2);
1294  SymEdge<T> *ldo;
1295  SymEdge<T> *ldi;
1296  SymEdge<T> *rdi;
1297  SymEdge<T> *rdo;
1298  dc_tri(cdt, sites, start, start + n2, &ldo, &ldi);
1299  dc_tri(cdt, sites, start + n2, end, &rdi, &rdo);
1300  if (dbg_level > 0) {
1301  std::cout << "\nDC_TRI merge step for start=" << start << ", end=" << end << "\n";
1302  std::cout << "ldo " << ldo << "\n"
1303  << "ldi " << ldi << "\n"
1304  << "rdi " << rdi << "\n"
1305  << "rdo " << rdo << "\n";
1306  if (dbg_level > 1) {
1307  std::string lab = "dc_tri (" + std::to_string(start) + "," + std::to_string(start + n2) +
1308  ")(" + std::to_string(start + n2) + "," + std::to_string(end) + ")";
1309  cdt_draw(lab, *cdt);
1310  }
1311  }
1312  /* Find lower common tangent of L and R. */
1313  for (;;) {
1314  if (vert_left_of_symedge(rdi->vert, ldi)) {
1315  ldi = ldi->next;
1316  }
1317  else if (vert_right_of_symedge(ldi->vert, rdi)) {
1318  rdi = sym(rdi)->rot; /* Previous edge to rdi with same right face. */
1319  }
1320  else {
1321  break;
1322  }
1323  }
1324  if (dbg_level > 0) {
1325  std::cout << "common lower tangent in between\n"
1326  << "rdi " << rdi << "\n"
1327  << "ldi" << ldi << "\n";
1328  }
1329 
1330  CDTEdge<T> *ebasel = cdt->connect_separate_parts(sym(rdi)->next, ldi);
1331  SymEdge<T> *basel = &ebasel->symedges[0];
1332  SymEdge<T> *basel_sym = &ebasel->symedges[1];
1333  if (dbg_level > 1) {
1334  std::cout << "basel " << basel;
1335  cdt_draw("after basel made", *cdt);
1336  }
1337  if (ldi->vert == ldo->vert) {
1338  ldo = basel_sym;
1339  }
1340  if (rdi->vert == rdo->vert) {
1341  rdo = basel;
1342  }
1343 
1344  /* Merge loop. */
1345  for (;;) {
1346  /* Locate the first point lcand->next->vert encountered by rising bubble,
1347  * and delete L edges out of basel->next->vert that fail the circle test. */
1348  SymEdge<T> *lcand = basel_sym->rot;
1349  SymEdge<T> *rcand = basel_sym->next;
1350  if (dbg_level > 1) {
1351  std::cout << "\ntop of merge loop\n";
1352  std::cout << "lcand " << lcand << "\n"
1353  << "rcand " << rcand << "\n"
1354  << "basel " << basel << "\n";
1355  }
1356  if (dc_tri_valid(lcand, basel, basel_sym)) {
1357  if (dbg_level > 1) {
1358  std::cout << "found valid lcand\n";
1359  std::cout << " lcand" << lcand << "\n";
1360  }
1361  while (filtered_incircle(basel_sym->vert->co,
1362  basel->vert->co,
1363  lcand->next->vert->co,
1364  lcand->rot->next->vert->co) > 0.0) {
1365  if (dbg_level > 1) {
1366  std::cout << "incircle says to remove lcand\n";
1367  std::cout << " lcand" << lcand << "\n";
1368  }
1369  SymEdge<T> *t = lcand->rot;
1370  cdt->delete_edge(sym(lcand));
1371  lcand = t;
1372  }
1373  }
1374  /* Symmetrically, locate first R point to be hit and delete R edges. */
1375  if (dc_tri_valid(rcand, basel, basel_sym)) {
1376  if (dbg_level > 1) {
1377  std::cout << "found valid rcand\n";
1378  std::cout << " rcand" << rcand << "\n";
1379  }
1380  while (filtered_incircle(basel_sym->vert->co,
1381  basel->vert->co,
1382  rcand->next->vert->co,
1383  sym(rcand)->next->next->vert->co) > 0.0) {
1384  if (dbg_level > 0) {
1385  std::cout << "incircle says to remove rcand\n";
1386  std::cout << " rcand" << rcand << "\n";
1387  }
1388  SymEdge<T> *t = sym(rcand)->next;
1389  cdt->delete_edge(rcand);
1390  rcand = t;
1391  }
1392  }
1393  /* If both lcand and rcand are invalid, then basel is the common upper tangent. */
1394  bool valid_lcand = dc_tri_valid(lcand, basel, basel_sym);
1395  bool valid_rcand = dc_tri_valid(rcand, basel, basel_sym);
1396  if (dbg_level > 0) {
1397  std::cout << "after bubbling up, valid_lcand=" << valid_lcand
1398  << ", valid_rand=" << valid_rcand << "\n";
1399  std::cout << "lcand" << lcand << "\n"
1400  << "rcand" << rcand << "\n";
1401  }
1402  if (!valid_lcand && !valid_rcand) {
1403  break;
1404  }
1405  /* The next cross edge to be connected is to either `lcand->next->vert` or `rcand->next->vert`;
1406  * if both are valid, choose the appropriate one using the #incircle test. */
1407  if (!valid_lcand || (valid_rcand && filtered_incircle(lcand->next->vert->co,
1408  lcand->vert->co,
1409  rcand->vert->co,
1410  rcand->next->vert->co) > 0)) {
1411  if (dbg_level > 0) {
1412  std::cout << "connecting rcand\n";
1413  std::cout << " se1=basel_sym" << basel_sym << "\n";
1414  std::cout << " se2=rcand->next" << rcand->next << "\n";
1415  }
1416  ebasel = cdt->add_diagonal(rcand->next, basel_sym);
1417  }
1418  else {
1419  if (dbg_level > 0) {
1420  std::cout << "connecting lcand\n";
1421  std::cout << " se1=sym(lcand)" << sym(lcand) << "\n";
1422  std::cout << " se2=basel_sym->next" << basel_sym->next << "\n";
1423  }
1424  ebasel = cdt->add_diagonal(basel_sym->next, sym(lcand));
1425  }
1426  basel = &ebasel->symedges[0];
1427  basel_sym = &ebasel->symedges[1];
1428  BLI_assert(basel_sym->face == cdt->outer_face);
1429  if (dbg_level > 2) {
1430  cdt_draw("after adding new crossedge", *cdt);
1431  }
1432  }
1433  *r_le = ldo;
1434  *r_re = rdo;
1435  BLI_assert(sym(ldo)->face == cdt->outer_face && rdo->face == cdt->outer_face);
1436 }
1437 
1438 /* Guibas-Stolfi Divide-and_Conquer algorithm. */
1439 template<typename T> void dc_triangulate(CDTArrangement<T> *cdt, Array<SiteInfo<T>> &sites)
1440 {
1441  /* Compress sites in place to eliminated verts that merge to others. */
1442  int i = 0;
1443  int j = 0;
1444  int nsites = sites.size();
1445  while (j < nsites) {
1446  /* Invariant: `sites[0..i-1]` have non-merged verts from `0..(j-1)` in them. */
1447  sites[i] = sites[j++];
1448  if (sites[i].v->merge_to_index < 0) {
1449  i++;
1450  }
1451  }
1452  int n = i;
1453  if (n == 0) {
1454  return;
1455  }
1456  SymEdge<T> *le;
1457  SymEdge<T> *re;
1458  dc_tri(cdt, sites, 0, n, &le, &re);
1459 }
1460 
1479 template<typename T> void initial_triangulation(CDTArrangement<T> *cdt)
1480 {
1481  int n = cdt->verts.size();
1482  if (n <= 1) {
1483  return;
1484  }
1485  Array<SiteInfo<T>> sites(n);
1486  for (int i = 0; i < n; ++i) {
1487  sites[i].v = cdt->verts[i];
1488  sites[i].orig_index = i;
1489  }
1490  std::sort(sites.begin(), sites.end(), site_lexicographic_sort<T>);
1491  find_site_merges(sites);
1492  dc_triangulate(cdt, sites);
1493 }
1494 
1502 template<typename T> static void re_delaunay_triangulate(CDTArrangement<T> *cdt, SymEdge<T> *se)
1503 {
1504  if (se->face == cdt->outer_face || sym(se)->face == cdt->outer_face) {
1505  return;
1506  }
1507  /* 'se' is a diagonal just added, and it is base of area to retriangulate (face on its left) */
1508  int count = 1;
1509  for (SymEdge<T> *ss = se->next; ss != se; ss = ss->next) {
1510  count++;
1511  }
1512  if (count <= 3) {
1513  return;
1514  }
1515  /* First and last are the SymEdges whose verts are first and last off of base,
1516  * continuing from 'se'. */
1517  SymEdge<T> *first = se->next->next;
1518  /* We want to make a triangle with 'se' as base and some other c as 3rd vertex. */
1519  CDTVert<T> *a = se->vert;
1520  CDTVert<T> *b = se->next->vert;
1521  CDTVert<T> *c = first->vert;
1522  SymEdge<T> *cse = first;
1523  for (SymEdge<T> *ss = first->next; ss != se; ss = ss->next) {
1524  CDTVert<T> *v = ss->vert;
1525  if (filtered_incircle(a->co, b->co, c->co, v->co) > 0) {
1526  c = v;
1527  cse = ss;
1528  }
1529  }
1530  /* Add diagonals necessary to make `abc` a triangle. */
1531  CDTEdge<T> *ebc = nullptr;
1532  CDTEdge<T> *eca = nullptr;
1533  if (!exists_edge(b, c)) {
1534  ebc = cdt->add_diagonal(se->next, cse);
1535  }
1536  if (!exists_edge(c, a)) {
1537  eca = cdt->add_diagonal(cse, se);
1538  }
1539  /* Now recurse. */
1540  if (ebc) {
1541  re_delaunay_triangulate(cdt, &ebc->symedges[1]);
1542  }
1543  if (eca) {
1544  re_delaunay_triangulate(cdt, &eca->symedges[1]);
1545  }
1546 }
1547 
1548 template<typename T> inline int tri_orient(const SymEdge<T> *t)
1549 {
1550  return filtered_orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co);
1551 }
1552 
1578 template<typename T> class CrossData {
1579  public:
1580  T lambda = T(0);
1581  CDTVert<T> *vert;
1582  SymEdge<T> *in;
1583  SymEdge<T> *out;
1584 
1585  CrossData() : lambda(T(0)), vert(nullptr), in(nullptr), out(nullptr)
1586  {
1587  }
1588  CrossData(T l, CDTVert<T> *v, SymEdge<T> *i, SymEdge<T> *o) : lambda(l), vert(v), in(i), out(o)
1589  {
1590  }
1591  CrossData(const CrossData &other)
1592  : lambda(other.lambda), vert(other.vert), in(other.in), out(other.out)
1593  {
1594  }
1595  CrossData(CrossData &&other) noexcept
1596  : lambda(std::move(other.lambda)),
1597  vert(std::move(other.vert)),
1598  in(std::move(other.in)),
1599  out(std::move(other.out))
1600  {
1601  }
1602  ~CrossData() = default;
1604  {
1605  if (this != &other) {
1606  lambda = other.lambda;
1607  vert = other.vert;
1608  in = other.in;
1609  out = other.out;
1610  }
1611  return *this;
1612  }
1613  CrossData &operator=(CrossData &&other) noexcept
1614  {
1615  lambda = std::move(other.lambda);
1616  vert = std::move(other.vert);
1617  in = std::move(other.in);
1618  out = std::move(other.out);
1619  return *this;
1620  }
1621 };
1622 
1623 template<typename T>
1624 bool get_next_crossing_from_vert(CDT_state<T> *cdt_state,
1625  CrossData<T> *cd,
1626  CrossData<T> *cd_next,
1627  const CDTVert<T> *v2);
1628 
1636 template<typename T>
1638  SymEdge<T> *cd_out,
1639  CrossData<T> *cd,
1640  CrossData<T> *cd_next)
1641 {
1642  SymEdge<T> *se;
1643 
1644  cd_next->lambda = T(0);
1645  cd_next->vert = v;
1646  cd_next->in = nullptr;
1647  cd_next->out = nullptr;
1648  if (cd->lambda == 0) {
1649  cd->out = cd_out;
1650  }
1651  else {
1652  /* One of the edges in the triangle with edge sym(cd->in) contains v. */
1653  se = sym(cd->in);
1654  if (se->vert != v) {
1655  se = se->next;
1656  if (se->vert != v) {
1657  se = se->next;
1658  }
1659  }
1660  BLI_assert(se->vert == v);
1661  cd_next->in = se;
1662  }
1663 }
1664 
1678 template<typename T>
1679 void fill_crossdata_for_intersect(const FatCo<T> &curco,
1680  const CDTVert<T> *v2,
1681  SymEdge<T> *t,
1682  CrossData<T> *cd,
1683  CrossData<T> *cd_next,
1684  const T epsilon)
1685 {
1686  CDTVert<T> *va = t->vert;
1687  CDTVert<T> *vb = t->next->vert;
1688  CDTVert<T> *vc = t->next->next->vert;
1689  SymEdge<T> *se_vcvb = sym(t->next);
1690  SymEdge<T> *se_vcva = t->next->next;
1691  BLI_assert(se_vcva->vert == vc && se_vcva->next->vert == va);
1692  BLI_assert(se_vcvb->vert == vc && se_vcvb->next->vert == vb);
1693  UNUSED_VARS_NDEBUG(vc);
1694  auto isect = isect_seg_seg(va->co.exact, vb->co.exact, curco.exact, v2->co.exact);
1695  T &lambda = isect.lambda;
1696  switch (isect.kind) {
1697  case isect_result<vec2<T>>::LINE_LINE_CROSS: {
1698 #ifdef WITH_GMP
1699  if (!std::is_same<T, mpq_class>::value) {
1700 #else
1701  if (true) {
1702 #endif
1703  double len_ab = distance(va->co.approx, vb->co.approx);
1704  if (lambda * len_ab <= epsilon) {
1705  fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
1706  }
1707  else if ((1 - lambda) * len_ab <= epsilon) {
1708  fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
1709  }
1710  else {
1711  *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
1712  if (cd->lambda == 0) {
1713  cd->out = se_vcva;
1714  }
1715  }
1716  }
1717  else {
1718  *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
1719  if (cd->lambda == 0) {
1720  cd->out = se_vcva;
1721  }
1722  }
1723  break;
1724  }
1725  case isect_result<vec2<T>>::LINE_LINE_EXACT: {
1726  if (lambda == 0) {
1727  fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
1728  }
1729  else if (lambda == 1) {
1730  fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
1731  }
1732  else {
1733  *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
1734  if (cd->lambda == 0) {
1735  cd->out = se_vcva;
1736  }
1737  }
1738  break;
1739  }
1740  case isect_result<vec2<T>>::LINE_LINE_NONE: {
1741 #ifdef WITH_GMP
1742  if (std::is_same<T, mpq_class>::value) {
1743  BLI_assert(false);
1744  }
1745 #endif
1746  /* It should be very near one end or other of segment. */
1747  const T middle_lambda = 0.5;
1748  if (lambda <= middle_lambda) {
1749  fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
1750  }
1751  else {
1752  fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
1753  }
1754  break;
1755  }
1756  case isect_result<vec2<T>>::LINE_LINE_COLINEAR: {
1757  if (distance_squared(va->co.approx, v2->co.approx) <=
1758  distance_squared(vb->co.approx, v2->co.approx)) {
1759  fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
1760  }
1761  else {
1762  fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
1763  }
1764  break;
1765  }
1766  }
1767 } // namespace blender::meshintersect
1768 
1776 template<typename T>
1777 bool get_next_crossing_from_vert(CDT_state<T> *cdt_state,
1778  CrossData<T> *cd,
1779  CrossData<T> *cd_next,
1780  const CDTVert<T> *v2)
1781 {
1782  SymEdge<T> *tstart = cd->vert->symedge;
1783  SymEdge<T> *t = tstart;
1784  CDTVert<T> *vcur = cd->vert;
1785  bool ok = false;
1786  do {
1787  /* The ray from `vcur` to v2 has to go either between two successive
1788  * edges around `vcur` or exactly along them. This time through the
1789  * loop, check to see if the ray goes along `vcur-va`
1790  * or between `vcur-va` and `vcur-vb`, where va is the end of t
1791  * and vb is the next vertex (on the next rot edge around vcur, but
1792  * should also be the next vert of triangle starting with `vcur-va`. */
1793  if (t->face != cdt_state->cdt.outer_face && tri_orient(t) < 0) {
1794  BLI_assert(false); /* Shouldn't happen. */
1795  }
1796  CDTVert<T> *va = t->next->vert;
1797  CDTVert<T> *vb = t->next->next->vert;
1798  int orient1 = filtered_orient2d(t->vert->co, va->co, v2->co);
1799  if (orient1 == 0 && in_line<T>(vcur->co, va->co, v2->co)) {
1800  fill_crossdata_for_through_vert(va, t, cd, cd_next);
1801  ok = true;
1802  break;
1803  }
1804  if (t->face != cdt_state->cdt.outer_face) {
1805  int orient2 = filtered_orient2d(vcur->co, vb->co, v2->co);
1806  /* Don't handle orient2 == 0 case here: next rotation will get it. */
1807  if (orient1 > 0 && orient2 < 0) {
1808  /* Segment intersection. */
1809  t = t->next;
1810  fill_crossdata_for_intersect(vcur->co, v2, t, cd, cd_next, cdt_state->epsilon);
1811  ok = true;
1812  break;
1813  }
1814  }
1815  } while ((t = t->rot) != tstart);
1816  return ok;
1817 }
1818 
1827 template<typename T>
1829  CrossData<T> *cd_next,
1830  const CDTVert<T> *v2,
1831  const T epsilon)
1832 {
1833  CDTVert<T> *va = cd->in->vert;
1834  CDTVert<T> *vb = cd->in->next->vert;
1835  vec2<T> curco = interpolate(va->co.exact, vb->co.exact, cd->lambda);
1836  FatCo<T> fat_curco(curco);
1837  SymEdge<T> *se_ac = sym(cd->in)->next;
1838  CDTVert<T> *vc = se_ac->next->vert;
1839  int orient = filtered_orient2d(fat_curco, v2->co, vc->co);
1840  if (orient < 0) {
1841  fill_crossdata_for_intersect<T>(fat_curco, v2, se_ac->next, cd, cd_next, epsilon);
1842  }
1843  else if (orient > 0.0) {
1844  fill_crossdata_for_intersect(fat_curco, v2, se_ac, cd, cd_next, epsilon);
1845  }
1846  else {
1847  *cd_next = CrossData<T>{0.0, vc, se_ac->next, nullptr};
1848  }
1849 }
1850 
1851 constexpr int inline_crossings_size = 128;
1852 template<typename T>
1853 void dump_crossings(const Vector<CrossData<T>, inline_crossings_size> &crossings)
1854 {
1855  std::cout << "CROSSINGS\n";
1856  for (int i = 0; i < crossings.size(); ++i) {
1857  std::cout << i << ": ";
1858  const CrossData<T> &cd = crossings[i];
1859  if (cd.lambda == 0) {
1860  std::cout << "v" << cd.vert->index;
1861  }
1862  else {
1863  std::cout << "lambda=" << cd.lambda;
1864  }
1865  if (cd.in != nullptr) {
1866  std::cout << " in=" << short_se_dump(cd.in);
1867  std::cout << " out=" << short_se_dump(cd.out);
1868  }
1869  std::cout << "\n";
1870  }
1871 }
1872 
1884 template<typename T>
1886  CDT_state<T> *cdt_state, CDTVert<T> *v1, CDTVert<T> *v2, int input_id, LinkNode **r_edges)
1887 {
1888  constexpr int dbg_level = 0;
1889  if (dbg_level > 0) {
1890  std::cout << "\nADD EDGE CONSTRAINT\n" << vertname(v1) << " " << vertname(v2) << "\n";
1891  }
1892  LinkNodePair edge_list = {nullptr, nullptr};
1893 
1894  if (r_edges) {
1895  *r_edges = nullptr;
1896  }
1897 
1898  /*
1899  * Handle two special cases first:
1900  * 1) The two end vertices are the same (can happen because of merging).
1901  * 2) There is already an edge between v1 and v2.
1902  */
1903  if (v1 == v2) {
1904  return;
1905  }
1906  SymEdge<T> *t = find_symedge_between_verts(v1, v2);
1907  if (t != nullptr) {
1908  /* Segment already there. */
1909  add_to_input_ids(t->edge->input_ids, input_id);
1910  if (r_edges != nullptr) {
1911  BLI_linklist_append(&edge_list, t->edge);
1912  *r_edges = edge_list.list;
1913  }
1914  return;
1915  }
1916 
1917  /*
1918  * Fill crossings array with CrossData points for intersection path from v1 to v2.
1919  *
1920  * At every point, the crossings array has the path so far, except that
1921  * the .out field of the last element of it may not be known yet -- if that
1922  * last element is a vertex, then we won't know the output edge until we
1923  * find the next crossing.
1924  *
1925  * To protect against infinite loops, we keep track of which vertices
1926  * we have visited by setting their visit_index to a new visit epoch.
1927  *
1928  * We check a special case first: where the segment is already there in
1929  * one hop. Saves a bunch of orient2d tests in that common case.
1930  */
1931  int visit = ++cdt_state->visit_count;
1932  Vector<CrossData<T>, inline_crossings_size> crossings;
1933  crossings.append(CrossData<T>(T(0), v1, nullptr, nullptr));
1934  int n;
1935  while (!((n = crossings.size()) > 0 && crossings[n - 1].vert == v2)) {
1936  crossings.append(CrossData<T>());
1937  CrossData<T> *cd = &crossings[n - 1];
1938  CrossData<T> *cd_next = &crossings[n];
1939  bool ok;
1940  if (crossings[n - 1].lambda == 0) {
1941  ok = get_next_crossing_from_vert(cdt_state, cd, cd_next, v2);
1942  }
1943  else {
1944  get_next_crossing_from_edge<T>(cd, cd_next, v2, cdt_state->epsilon);
1945  ok = true;
1946  }
1947  constexpr int unreasonably_large_crossings = 100000;
1948  if (!ok || crossings.size() == unreasonably_large_crossings) {
1949  /* Shouldn't happen but if does, just bail out. */
1950  BLI_assert(false);
1951  return;
1952  }
1953  if (crossings[n].lambda == 0) {
1954  if (crossings[n].vert->visit_index == visit) {
1955  /* Shouldn't happen but if it does, just bail out. */
1956  BLI_assert(false);
1957  return;
1958  }
1959  crossings[n].vert->visit_index = visit;
1960  }
1961  }
1962 
1963  if (dbg_level > 0) {
1964  dump_crossings(crossings);
1965  }
1966 
1967  /*
1968  * Post-process crossings.
1969  * Some crossings may have an intersection crossing followed
1970  * by a vertex crossing that is on the same edge that was just
1971  * intersected. We prefer to go directly from the previous
1972  * crossing directly to the vertex. This may chain backwards.
1973  *
1974  * This loop marks certain crossings as "deleted", by setting
1975  * their lambdas to -1.0.
1976  */
1977  int ncrossings = crossings.size();
1978  for (int i = 2; i < ncrossings; ++i) {
1979  CrossData<T> *cd = &crossings[i];
1980  if (cd->lambda == 0.0) {
1981  CDTVert<T> *v = cd->vert;
1982  int j;
1983  CrossData<T> *cd_prev;
1984  for (j = i - 1; j > 0; --j) {
1985  cd_prev = &crossings[j];
1986  if ((cd_prev->lambda == 0.0 && cd_prev->vert != v) ||
1987  (cd_prev->lambda != 0.0 && cd_prev->in->vert != v && cd_prev->in->next->vert != v)) {
1988  break;
1989  }
1990  cd_prev->lambda = -1.0; /* Mark cd_prev as 'deleted'. */
1991  }
1992  if (j < i - 1) {
1993  /* Some crossings were deleted. Fix the in and out edges across gap. */
1994  cd_prev = &crossings[j];
1995  SymEdge<T> *se;
1996  if (cd_prev->lambda == 0.0) {
1997  se = find_symedge_between_verts(cd_prev->vert, v);
1998  if (se == nullptr) {
1999  return;
2000  }
2001  cd_prev->out = se;
2002  cd->in = nullptr;
2003  }
2004  else {
2005  se = find_symedge_with_face(v, sym(cd_prev->in)->face);
2006  if (se == nullptr) {
2007  return;
2008  }
2009  cd->in = se;
2010  }
2011  }
2012  }
2013  }
2014 
2015  /*
2016  * Insert all intersection points on constrained edges.
2017  */
2018  for (int i = 0; i < ncrossings; ++i) {
2019  CrossData<T> *cd = &crossings[i];
2020  if (cd->lambda != 0.0 && cd->lambda != -1.0 && is_constrained_edge(cd->in->edge)) {
2021  CDTEdge<T> *edge = cdt_state->cdt.split_edge(cd->in, cd->lambda);
2022  cd->vert = edge->symedges[0].vert;
2023  }
2024  }
2025 
2026  /*
2027  * Remove any crossed, non-intersected edges.
2028  */
2029  for (int i = 0; i < ncrossings; ++i) {
2030  CrossData<T> *cd = &crossings[i];
2031  if (cd->lambda != 0.0 && cd->lambda != -1.0 && !is_constrained_edge(cd->in->edge)) {
2032  cdt_state->cdt.delete_edge(cd->in);
2033  }
2034  }
2035 
2036  /*
2037  * Insert segments for v1->v2.
2038  */
2039  SymEdge<T> *tstart = crossings[0].out;
2040  for (int i = 1; i < ncrossings; i++) {
2041  CrossData<T> *cd = &crossings[i];
2042  if (cd->lambda == -1.0) {
2043  continue; /* This crossing was deleted. */
2044  }
2045  t = nullptr;
2046  SymEdge<T> *tnext = t;
2047  CDTEdge<T> *edge;
2048  if (cd->lambda != 0.0) {
2049  if (is_constrained_edge(cd->in->edge)) {
2050  t = cd->vert->symedge;
2051  tnext = sym(t)->next;
2052  }
2053  }
2054  else if (cd->lambda == 0.0) {
2055  t = cd->in;
2056  tnext = cd->out;
2057  if (t == nullptr) {
2058  /* Previous non-deleted crossing must also have been a vert, and segment should exist. */
2059  int j;
2060  CrossData<T> *cd_prev;
2061  for (j = i - 1; j >= 0; j--) {
2062  cd_prev = &crossings[j];
2063  if (cd_prev->lambda != -1.0) {
2064  break;
2065  }
2066  }
2067  BLI_assert(cd_prev->lambda == 0.0);
2068  BLI_assert(cd_prev->out->next->vert == cd->vert);
2069  edge = cd_prev->out->edge;
2070  add_to_input_ids(edge->input_ids, input_id);
2071  if (r_edges != nullptr) {
2072  BLI_linklist_append(&edge_list, edge);
2073  }
2074  }
2075  }
2076  if (t != nullptr) {
2077  if (tstart->next->vert == t->vert) {
2078  edge = tstart->edge;
2079  }
2080  else {
2081  edge = cdt_state->cdt.add_diagonal(tstart, t);
2082  }
2083  add_to_input_ids(edge->input_ids, input_id);
2084  if (r_edges != nullptr) {
2085  BLI_linklist_append(&edge_list, edge);
2086  }
2087  /* Now retriangulate upper and lower gaps. */
2088  re_delaunay_triangulate(&cdt_state->cdt, &edge->symedges[0]);
2089  re_delaunay_triangulate(&cdt_state->cdt, &edge->symedges[1]);
2090  }
2091  if (i < ncrossings - 1) {
2092  if (tnext != nullptr) {
2093  tstart = tnext;
2094  }
2095  }
2096  }
2097 
2098  if (r_edges) {
2099  *r_edges = edge_list.list;
2100  }
2101 }
2102 
2109 template<typename T> void add_edge_constraints(CDT_state<T> *cdt_state, const CDT_input<T> &input)
2110 {
2111  int ne = input.edge.size();
2112  int nv = input.vert.size();
2113  for (int i = 0; i < ne; i++) {
2114  int iv1 = input.edge[i].first;
2115  int iv2 = input.edge[i].second;
2116  if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
2117  /* Ignore invalid indices in edges. */
2118  continue;
2119  }
2120  CDTVert<T> *v1 = cdt_state->cdt.get_vert_resolve_merge(iv1);
2121  CDTVert<T> *v2 = cdt_state->cdt.get_vert_resolve_merge(iv2);
2122  int id = cdt_state->need_ids ? i : 0;
2123  add_edge_constraint(cdt_state, v1, v2, id, nullptr);
2124  }
2125  cdt_state->face_edge_offset = ne;
2126 }
2127 
2146 template<typename T>
2148  CDT_state<T> *cdt_state, SymEdge<T> *face_symedge, int face_id, int fedge_start, int fedge_end)
2149 {
2150  /* Can't loop forever since eventually would visit every face. */
2151  cdt_state->visit_count++;
2152  int visit = cdt_state->visit_count;
2153  Vector<SymEdge<T> *> stack;
2154  stack.append(face_symedge);
2155  while (!stack.is_empty()) {
2156  SymEdge<T> *se = stack.pop_last();
2157  CDTFace<T> *face = se->face;
2158  if (face->visit_index == visit) {
2159  continue;
2160  }
2161  face->visit_index = visit;
2162  add_to_input_ids(face->input_ids, face_id);
2163  SymEdge<T> *se_start = se;
2164  for (se = se->next; se != se_start; se = se->next) {
2165  if (!id_range_in_list(se->edge->input_ids, fedge_start, fedge_end)) {
2166  SymEdge<T> *se_sym = sym(se);
2167  CDTFace<T> *face_other = se_sym->face;
2168  if (face_other->visit_index != visit) {
2169  stack.append(se_sym);
2170  }
2171  }
2172  }
2173  }
2174 }
2175 
2176 /* Return a power of 10 that is greater than or equal to x. */
2178 {
2179  if (x <= 0) {
2180  return 1;
2181  }
2182  int ans = 1;
2183  BLI_assert(x < INT_MAX / 10);
2184  while (ans < x) {
2185  ans *= 10;
2186  }
2187  return ans;
2188 }
2189 
2200 template<typename T>
2201 int add_face_constraints(CDT_state<T> *cdt_state,
2202  const CDT_input<T> &input,
2203  CDT_output_type output_type)
2204 {
2205  int nv = input.vert.size();
2206  int nf = input.face.size();
2207  SymEdge<T> *face_symedge0 = nullptr;
2208  CDTArrangement<T> *cdt = &cdt_state->cdt;
2209  int maxflen = 0;
2210  for (int f = 0; f < nf; f++) {
2211  maxflen = max_ii(maxflen, input.face[f].size());
2212  }
2213  /* For convenience in debugging, make face_edge_offset be a power of 10. */
2214  cdt_state->face_edge_offset = power_of_10_greater_equal_to(
2215  max_ii(maxflen, cdt_state->face_edge_offset));
2216  /* The original_edge encoding scheme doesn't work if the following is false.
2217  * If we really have that many faces and that large a max face length that when multiplied
2218  * together the are >= INT_MAX, then the Delaunay calculation will take unreasonably long anyway.
2219  */
2220  BLI_assert(INT_MAX / cdt_state->face_edge_offset > nf);
2221  int faces_added = 0;
2222  for (int f = 0; f < nf; f++) {
2223  int flen = input.face[f].size();
2224  if (flen <= 2) {
2225  /* Ignore faces with fewer than 3 vertices. */
2226  continue;
2227  }
2228  int fedge_start = (f + 1) * cdt_state->face_edge_offset;
2229  for (int i = 0; i < flen; i++) {
2230  int face_edge_id = fedge_start + i;
2231  int iv1 = input.face[f][i];
2232  int iv2 = input.face[f][(i + 1) % flen];
2233  if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
2234  /* Ignore face edges with invalid vertices. */
2235  continue;
2236  }
2237  ++faces_added;
2238  CDTVert<T> *v1 = cdt->get_vert_resolve_merge(iv1);
2239  CDTVert<T> *v2 = cdt->get_vert_resolve_merge(iv2);
2240  LinkNode *edge_list;
2241  int id = cdt_state->need_ids ? face_edge_id : 0;
2242  add_edge_constraint(cdt_state, v1, v2, id, &edge_list);
2243  /* Set a new face_symedge0 each time since earlier ones may not
2244  * survive later symedge splits. Really, just want the one when
2245  * i == flen -1, but this code guards against that one somehow
2246  * being null.
2247  */
2248  if (edge_list != nullptr) {
2249  CDTEdge<T> *face_edge = static_cast<CDTEdge<T> *>(edge_list->link);
2250  face_symedge0 = &face_edge->symedges[0];
2251  if (face_symedge0->vert != v1) {
2252  face_symedge0 = &face_edge->symedges[1];
2253  BLI_assert(face_symedge0->vert == v1);
2254  }
2255  }
2256  BLI_linklist_free(edge_list, nullptr);
2257  }
2258  int fedge_end = fedge_start + flen - 1;
2259  if (face_symedge0 != nullptr) {
2260  /* We need to propagate face ids to all faces that represent #f, if #need_ids.
2261  * Even if `need_ids == false`, we need to propagate at least the fact that
2262  * the face ids set would be non-empty if the output type is one of the ones
2263  * making valid BMesh faces. */
2264  int id = cdt_state->need_ids ? f : 0;
2265  add_face_ids(cdt_state, face_symedge0, id, fedge_start, fedge_end);
2266  if (cdt_state->need_ids ||
2268  add_face_ids(cdt_state, face_symedge0, f, fedge_start, fedge_end);
2269  }
2270  }
2271  }
2272  return faces_added;
2273 }
2274 
2275 /* Delete_edge but try not to mess up outer face.
2276  * Also faces have symedges now, so make sure not
2277  * to mess those up either. */
2278 template<typename T> void dissolve_symedge(CDT_state<T> *cdt_state, SymEdge<T> *se)
2279 {
2280  CDTArrangement<T> *cdt = &cdt_state->cdt;
2281  SymEdge<T> *symse = sym(se);
2282  if (symse->face == cdt->outer_face) {
2283  se = sym(se);
2284  symse = sym(se);
2285  }
2286  if (ELEM(cdt->outer_face->symedge, se, symse)) {
2287  /* Advancing by 2 to get past possible 'sym(se)'. */
2288  if (se->next->next == se) {
2289  cdt->outer_face->symedge = nullptr;
2290  }
2291  else {
2292  cdt->outer_face->symedge = se->next->next;
2293  }
2294  }
2295  else {
2296  if (se->face->symedge == se) {
2297  se->face->symedge = se->next;
2298  }
2299  if (symse->face->symedge == symse) {
2300  symse->face->symedge = symse->next;
2301  }
2302  }
2303  cdt->delete_edge(se);
2304 }
2305 
2309 template<typename T> void remove_non_constraint_edges(CDT_state<T> *cdt_state)
2310 {
2311  for (CDTEdge<T> *e : cdt_state->cdt.edges) {
2312  SymEdge<T> *se = &e->symedges[0];
2313  if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
2314  dissolve_symedge(cdt_state, se);
2315  }
2316  }
2317 }
2318 
2319 /*
2320  * Remove the non-constraint edges, but leave enough of them so that all of the
2321  * faces that would be #BMesh faces (that is, the faces that have some input representative)
2322  * are valid: they can't have holes, they can't have repeated vertices, and they can't have
2323  * repeated edges.
2324  *
2325  * Not essential, but to make the result look more aesthetically nice,
2326  * remove the edges in order of decreasing length, so that it is more likely that the
2327  * final remaining support edges are short, and therefore likely to make a fairly
2328  * direct path from an outer face to an inner hole face.
2329  */
2330 
2334 template<typename T> struct EdgeToSort {
2335  double len_squared = 0.0;
2336  CDTEdge<T> *e{nullptr};
2337 
2338  EdgeToSort() = default;
2339  EdgeToSort(const EdgeToSort &other) : len_squared(other.len_squared), e(other.e)
2340  {
2341  }
2342  EdgeToSort(EdgeToSort &&other) noexcept : len_squared(std::move(other.len_squared)), e(other.e)
2343  {
2344  }
2345  ~EdgeToSort() = default;
2347  {
2348  if (this != &other) {
2349  len_squared = other.len_squared;
2350  e = other.e;
2351  }
2352  return *this;
2353  }
2355  {
2356  len_squared = std::move(other.len_squared);
2357  e = other.e;
2358  return *this;
2359  }
2360 };
2361 
2362 template<typename T> void remove_non_constraint_edges_leave_valid_bmesh(CDT_state<T> *cdt_state)
2363 {
2364  CDTArrangement<T> *cdt = &cdt_state->cdt;
2365  size_t nedges = cdt->edges.size();
2366  if (nedges == 0) {
2367  return;
2368  }
2369  Vector<EdgeToSort<T>> dissolvable_edges;
2370  dissolvable_edges.reserve(cdt->edges.size());
2371  int i = 0;
2372  for (CDTEdge<T> *e : cdt->edges) {
2373  if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
2374  dissolvable_edges.append(EdgeToSort<T>());
2375  dissolvable_edges[i].e = e;
2376  const vec2<double> &co1 = e->symedges[0].vert->co.approx;
2377  const vec2<double> &co2 = e->symedges[1].vert->co.approx;
2378  dissolvable_edges[i].len_squared = distance_squared(co1, co2);
2379  i++;
2380  }
2381  }
2382  std::sort(dissolvable_edges.begin(),
2383  dissolvable_edges.end(),
2384  [](const EdgeToSort<T> &a, const EdgeToSort<T> &b) -> bool {
2385  return (a.len_squared < b.len_squared);
2386  });
2387  for (EdgeToSort<T> &ets : dissolvable_edges) {
2388  CDTEdge<T> *e = ets.e;
2389  SymEdge<T> *se = &e->symedges[0];
2390  bool dissolve = true;
2391  CDTFace<T> *fleft = se->face;
2392  CDTFace<T> *fright = sym(se)->face;
2393  if (fleft != cdt->outer_face && fright != cdt->outer_face &&
2394  (fleft->input_ids.size() > 0 || fright->input_ids.size() > 0)) {
2395  /* Is there another #SymEdge with same left and right faces?
2396  * Or is there a vertex not part of e touching the same left and right faces? */
2397  for (SymEdge<T> *se2 = se->next; dissolve && se2 != se; se2 = se2->next) {
2398  if (sym(se2)->face == fright ||
2399  (se2->vert != se->next->vert && vert_touches_face(se2->vert, fright))) {
2400  dissolve = false;
2401  }
2402  }
2403  }
2404 
2405  if (dissolve) {
2406  dissolve_symedge(cdt_state, se);
2407  }
2408  }
2409 }
2410 
2411 template<typename T> void remove_outer_edges_until_constraints(CDT_state<T> *cdt_state)
2412 {
2413  int visit = ++cdt_state->visit_count;
2414 
2415  cdt_state->cdt.outer_face->visit_index = visit;
2416  /* Walk around outer face, adding faces on other side of dissolvable edges to stack. */
2417  Vector<CDTFace<T> *> fstack;
2418  SymEdge<T> *se_start = cdt_state->cdt.outer_face->symedge;
2419  SymEdge<T> *se = se_start;
2420  do {
2421  if (!is_constrained_edge(se->edge)) {
2422  CDTFace<T> *fsym = sym(se)->face;
2423  if (fsym->visit_index != visit) {
2424  fstack.append(fsym);
2425  }
2426  }
2427  } while ((se = se->next) != se_start);
2428 
2429  while (!fstack.is_empty()) {
2430  LinkNode *to_dissolve = nullptr;
2431  bool dissolvable;
2432  CDTFace<T> *f = fstack.pop_last();
2433  if (f->visit_index == visit) {
2434  continue;
2435  }
2436  BLI_assert(f != cdt_state->cdt.outer_face);
2437  f->visit_index = visit;
2438  se_start = se = f->symedge;
2439  do {
2440  dissolvable = !is_constrained_edge(se->edge);
2441  if (dissolvable) {
2442  CDTFace<T> *fsym = sym(se)->face;
2443  if (fsym->visit_index != visit) {
2444  fstack.append(fsym);
2445  }
2446  else {
2447  BLI_linklist_prepend(&to_dissolve, se);
2448  }
2449  }
2450  se = se->next;
2451  } while (se != se_start);
2452  while (to_dissolve != nullptr) {
2453  se = static_cast<SymEdge<T> *>(BLI_linklist_pop(&to_dissolve));
2454  if (se->next != nullptr) {
2455  dissolve_symedge(cdt_state, se);
2456  }
2457  }
2458  }
2459 }
2460 
2461 template<typename T> void remove_faces_in_holes(CDT_state<T> *cdt_state)
2462 {
2463  CDTArrangement<T> *cdt = &cdt_state->cdt;
2464  for (int i : cdt->faces.index_range()) {
2465  CDTFace<T> *f = cdt->faces[i];
2466  if (!f->deleted && f->hole) {
2467  f->deleted = true;
2468  SymEdge<T> *se = f->symedge;
2469  SymEdge<T> *se_start = se;
2470  SymEdge<T> *se_next = nullptr;
2471  do {
2472  BLI_assert(se != nullptr);
2473  se_next = se->next; /* In case we delete this edge. */
2474  if (se->edge && !is_constrained_edge(se->edge)) {
2475  /* Invalidate one half of this edge. The other has will be or has been
2476  * handled with the adjacent triangle is processed: it should be part of the same hole.
2477  */
2478  se->next = nullptr;
2479  }
2480  se = se_next;
2481  } while (se != se_start);
2482  }
2483  }
2484 }
2485 
2496 template<typename T> void detect_holes(CDT_state<T> *cdt_state)
2497 {
2498  CDTArrangement<T> *cdt = &cdt_state->cdt;
2499 
2500  /* Make it so that each face with the same visit_index is connected through a path of
2501  * non-constraint edges. */
2502  Vector<CDTFace<T> *> fstack;
2503  Vector<CDTFace<T> *> region_rep_face;
2504  for (int i : cdt->faces.index_range()) {
2505  cdt->faces[i]->visit_index = -1;
2506  }
2507  int cur_region = -1;
2508  cdt->outer_face->visit_index = -2; /* Don't visit this one. */
2509  for (int i : cdt->faces.index_range()) {
2510  CDTFace<T> *f = cdt->faces[i];
2511  if (!f->deleted && f->symedge && f->visit_index == -1) {
2512  fstack.append(f);
2513  ++cur_region;
2514  region_rep_face.append(f);
2515  while (!fstack.is_empty()) {
2516  CDTFace<T> *f = fstack.pop_last();
2517  if (f->visit_index != -1) {
2518  continue;
2519  }
2520  f->visit_index = cur_region;
2521  SymEdge<T> *se_start = f->symedge;
2522  SymEdge<T> *se = se_start;
2523  do {
2524  if (se->edge && !is_constrained_edge(se->edge)) {
2525  CDTFace<T> *fsym = sym(se)->face;
2526  if (fsym && !fsym->deleted && fsym->visit_index == -1) {
2527  fstack.append(fsym);
2528  }
2529  }
2530  se = se->next;
2531  } while (se != se_start);
2532  }
2533  }
2534  }
2535  cdt_state->visit_count = ++cur_region; /* Good start for next use of visit_count. */
2536 
2537  /* Now get hole status for each region_rep_face. */
2538 
2539  /* Pick a ray end almost certain to be outside everything and in direction
2540  * that is unlikely to hit a vertex or overlap an edge exactly. */
2541  FatCo<T> ray_end;
2542  ray_end.exact = vec2<T>(123456, 654321);
2543  for (int i : region_rep_face.index_range()) {
2544  CDTFace<T> *f = region_rep_face[i];
2545  FatCo<T> mid;
2546  mid.exact[0] = (f->symedge->vert->co.exact[0] + f->symedge->next->vert->co.exact[0] +
2547  f->symedge->next->next->vert->co.exact[0]) /
2548  3;
2549  mid.exact[1] = (f->symedge->vert->co.exact[1] + f->symedge->next->vert->co.exact[1] +
2550  f->symedge->next->next->vert->co.exact[1]) /
2551  3;
2552  std::atomic<int> hits = 0;
2553  /* TODO: Use CDT data structure here to greatly reduce search for intersections! */
2554  threading::parallel_for(cdt->edges.index_range(), 256, [&](IndexRange range) {
2555  for (const int i : range) {
2556  const CDTEdge<T> *e = cdt->edges[i];
2557  if (!is_deleted_edge(e) && is_constrained_edge(e)) {
2558  if (e->symedges[0].face->visit_index == e->symedges[1].face->visit_index) {
2559  continue; /* Don't count hits on edges between faces in same region. */
2560  }
2561  auto isect = isect_seg_seg(ray_end.exact,
2562  mid.exact,
2563  e->symedges[0].vert->co.exact,
2564  e->symedges[1].vert->co.exact);
2565  switch (isect.kind) {
2566  case isect_result<vec2<T>>::LINE_LINE_CROSS: {
2567  hits++;
2568  break;
2569  }
2570  case isect_result<vec2<T>>::LINE_LINE_EXACT:
2571  case isect_result<vec2<T>>::LINE_LINE_NONE:
2572  case isect_result<vec2<T>>::LINE_LINE_COLINEAR:
2573  break;
2574  }
2575  }
2576  }
2577  });
2578  f->hole = (hits.load() % 2) == 0;
2579  }
2580 
2581  /* Finally, propagate hole status to all holes of a region. */
2582  for (int i : cdt->faces.index_range()) {
2583  CDTFace<T> *f = cdt->faces[i];
2584  int region = f->visit_index;
2585  if (region < 0) {
2586  continue;
2587  }
2588  CDTFace<T> *f_region_rep = region_rep_face[region];
2589  if (i >= 0) {
2590  f->hole = f_region_rep->hole;
2591  }
2592  }
2593 }
2594 
2599 template<typename T>
2600 void prepare_cdt_for_output(CDT_state<T> *cdt_state, const CDT_output_type output_type)
2601 {
2602  CDTArrangement<T> *cdt = &cdt_state->cdt;
2603  if (cdt->edges.is_empty()) {
2604  return;
2605  }
2606 
2607  /* Make sure all non-deleted faces have a symedge. */
2608  for (CDTEdge<T> *e : cdt->edges) {
2609  if (!is_deleted_edge(e)) {
2610  if (e->symedges[0].face->symedge == nullptr) {
2611  e->symedges[0].face->symedge = &e->symedges[0];
2612  }
2613  if (e->symedges[1].face->symedge == nullptr) {
2614  e->symedges[1].face->symedge = &e->symedges[1];
2615  }
2616  }
2617  }
2618 
2619  bool need_holes = output_type == CDT_INSIDE_WITH_HOLES ||
2621 
2622  if (need_holes) {
2623  detect_holes(cdt_state);
2624  }
2625 
2626  if (output_type == CDT_CONSTRAINTS) {
2627  remove_non_constraint_edges(cdt_state);
2628  }
2629  else if (output_type == CDT_CONSTRAINTS_VALID_BMESH) {
2631  }
2632  else if (output_type == CDT_INSIDE) {
2634  }
2635  else if (output_type == CDT_INSIDE_WITH_HOLES) {
2637  remove_faces_in_holes(cdt_state);
2638  }
2639  else if (output_type == CDT_CONSTRAINTS_VALID_BMESH_WITH_HOLES) {
2642  remove_faces_in_holes(cdt_state);
2643  }
2644 }
2645 
2646 template<typename T>
2647 CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
2648  const CDT_input<T> UNUSED(input),
2649  CDT_output_type output_type)
2650 {
2651  CDT_output_type oty = output_type;
2652  prepare_cdt_for_output(cdt_state, oty);
2654  CDTArrangement<T> *cdt = &cdt_state->cdt;
2655  result.face_edge_offset = cdt_state->face_edge_offset;
2656 
2657  /* All verts without a merge_to_index will be output.
2658  * vert_to_output_map[i] will hold the output vertex index
2659  * corresponding to the vert in position i in cdt->verts.
2660  * This first loop sets vert_to_output_map for un-merged verts. */
2661  int verts_size = cdt->verts.size();
2662  Array<int> vert_to_output_map(verts_size);
2663  int nv = 0;
2664  for (int i = 0; i < verts_size; ++i) {
2665  CDTVert<T> *v = cdt->verts[i];
2666  if (v->merge_to_index == -1) {
2667  vert_to_output_map[i] = nv;
2668  ++nv;
2669  }
2670  }
2671  if (nv <= 0) {
2672  return result;
2673  }
2674  /* Now we can set vert_to_output_map for merged verts,
2675  * and also add the input indices of merged verts to the input_ids
2676  * list of the merge target if they were an original input id. */
2677  if (nv < verts_size) {
2678  for (int i = 0; i < verts_size; ++i) {
2679  CDTVert<T> *v = cdt->verts[i];
2680  if (v->merge_to_index != -1) {
2681  if (cdt_state->need_ids) {
2682  if (i < cdt_state->input_vert_num) {
2683  add_to_input_ids(cdt->verts[v->merge_to_index]->input_ids, i);
2684  }
2685  }
2686  vert_to_output_map[i] = vert_to_output_map[v->merge_to_index];
2687  }
2688  }
2689  }
2690  result.vert = Array<vec2<T>>(nv);
2691  if (cdt_state->need_ids) {
2692  result.vert_orig = Array<Vector<int>>(nv);
2693  }
2694  int i_out = 0;
2695  for (int i = 0; i < verts_size; ++i) {
2696  CDTVert<T> *v = cdt->verts[i];
2697  if (v->merge_to_index == -1) {
2698  result.vert[i_out] = v->co.exact;
2699  if (cdt_state->need_ids) {
2700  if (i < cdt_state->input_vert_num) {
2701  result.vert_orig[i_out].append(i);
2702  }
2703  for (int vert : v->input_ids) {
2704  result.vert_orig[i_out].append(vert);
2705  }
2706  }
2707  ++i_out;
2708  }
2709  }
2710 
2711  /* All non-deleted edges will be output. */
2712  int ne = std::count_if(cdt->edges.begin(), cdt->edges.end(), [](const CDTEdge<T> *e) -> bool {
2713  return !is_deleted_edge(e);
2714  });
2715  result.edge = Array<std::pair<int, int>>(ne);
2716  if (cdt_state->need_ids) {
2717  result.edge_orig = Array<Vector<int>>(ne);
2718  }
2719  int e_out = 0;
2720  for (const CDTEdge<T> *e : cdt->edges) {
2721  if (!is_deleted_edge(e)) {
2722  int vo1 = vert_to_output_map[e->symedges[0].vert->index];
2723  int vo2 = vert_to_output_map[e->symedges[1].vert->index];
2724  result.edge[e_out] = std::pair<int, int>(vo1, vo2);
2725  if (cdt_state->need_ids) {
2726  for (int edge : e->input_ids) {
2727  result.edge_orig[e_out].append(edge);
2728  }
2729  }
2730  ++e_out;
2731  }
2732  }
2733 
2734  /* All non-deleted, non-outer faces will be output. */
2735  int nf = std::count_if(cdt->faces.begin(), cdt->faces.end(), [=](const CDTFace<T> *f) -> bool {
2736  return !f->deleted && f != cdt->outer_face;
2737  });
2738  result.face = Array<Vector<int>>(nf);
2739  if (cdt_state->need_ids) {
2740  result.face_orig = Array<Vector<int>>(nf);
2741  }
2742  int f_out = 0;
2743  for (const CDTFace<T> *f : cdt->faces) {
2744  if (!f->deleted && f != cdt->outer_face) {
2745  SymEdge<T> *se = f->symedge;
2746  BLI_assert(se != nullptr);
2747  SymEdge<T> *se_start = se;
2748  do {
2749  result.face[f_out].append(vert_to_output_map[se->vert->index]);
2750  se = se->next;
2751  } while (se != se_start);
2752  if (cdt_state->need_ids) {
2753  for (int face : f->input_ids) {
2754  result.face_orig[f_out].append(face);
2755  }
2756  }
2757  ++f_out;
2758  }
2759  }
2760  return result;
2761 }
2762 
2767 template<typename T> void add_input_verts(CDT_state<T> *cdt_state, const CDT_input<T> &input)
2768 {
2769  for (int i = 0; i < cdt_state->input_vert_num; ++i) {
2770  cdt_state->cdt.add_vert(input.vert[i]);
2771  }
2772 }
2773 
2774 template<typename T>
2776 {
2777  int nv = input.vert.size();
2778  int ne = input.edge.size();
2779  int nf = input.face.size();
2780  CDT_state<T> cdt_state(nv, ne, nf, input.epsilon, input.need_ids);
2781  add_input_verts(&cdt_state, input);
2782  initial_triangulation(&cdt_state.cdt);
2783  add_edge_constraints(&cdt_state, input);
2784  int actual_nf = add_face_constraints(&cdt_state, input, output_type);
2785  if (actual_nf == 0 && !ELEM(output_type, CDT_FULL, CDT_INSIDE, CDT_CONSTRAINTS)) {
2786  /* Can't look for faces or holes if there were no valid input faces. */
2787  output_type = CDT_INSIDE;
2788  }
2789  return get_cdt_output(&cdt_state, input, output_type);
2790 }
2791 
2792 blender::meshintersect::CDT_result<double> delaunay_2d_calc(const CDT_input<double> &input,
2793  CDT_output_type output_type)
2794 {
2795  return delaunay_calc(input, output_type);
2796 }
2797 
2798 #ifdef WITH_GMP
2799 blender::meshintersect::CDT_result<mpq_class> delaunay_2d_calc(const CDT_input<mpq_class> &input,
2800  CDT_output_type output_type)
2801 {
2802  return delaunay_calc(input, output_type);
2803 }
2804 #endif
2805 
2806 } /* namespace blender::meshintersect */
2807 
2808 /* C interface. */
2809 
2818  const CDT_output_type output_type)
2819 {
2820  blender::meshintersect::CDT_input<double> in;
2822  in.edge = blender::Array<std::pair<int, int>>(input->edges_len);
2823  in.face = blender::Array<blender::Vector<int>>(input->faces_len);
2824  for (int v = 0; v < input->verts_len; ++v) {
2825  double x = static_cast<double>(input->vert_coords[v][0]);
2826  double y = static_cast<double>(input->vert_coords[v][1]);
2827  in.vert[v] = blender::meshintersect::vec2<double>(x, y);
2828  }
2829  for (int e = 0; e < input->edges_len; ++e) {
2830  in.edge[e] = std::pair<int, int>(input->edges[e][0], input->edges[e][1]);
2831  }
2832  for (int f = 0; f < input->faces_len; ++f) {
2833  in.face[f] = blender::Vector<int>(input->faces_len_table[f]);
2834  int fstart = input->faces_start_table[f];
2835  for (int j = 0; j < input->faces_len_table[f]; ++j) {
2836  in.face[f][j] = input->faces[fstart + j];
2837  }
2838  }
2839  in.epsilon = static_cast<double>(input->epsilon);
2840  in.need_ids = input->need_ids;
2841 
2842  blender::meshintersect::CDT_result<double> res = blender::meshintersect::delaunay_2d_calc(
2843  in, output_type);
2844 
2845  ::CDT_result *output = static_cast<::CDT_result *>(MEM_mallocN(sizeof(*output), __func__));
2846  int nv = output->verts_len = res.vert.size();
2847  int ne = output->edges_len = res.edge.size();
2848  int nf = output->faces_len = res.face.size();
2849  int tot_v_orig = 0;
2850  int tot_e_orig = 0;
2851  int tot_f_orig = 0;
2852  int tot_f_lens = 0;
2853  if (input->need_ids) {
2854  for (int v = 0; v < nv; ++v) {
2855  tot_v_orig += res.vert_orig[v].size();
2856  }
2857  for (int e = 0; e < ne; ++e) {
2858  tot_e_orig += res.edge_orig[e].size();
2859  }
2860  }
2861  for (int f = 0; f < nf; ++f) {
2862  if (input->need_ids) {
2863  tot_f_orig += res.face_orig[f].size();
2864  }
2865  tot_f_lens += res.face[f].size();
2866  }
2867 
2868  output->vert_coords = static_cast<decltype(output->vert_coords)>(
2869  MEM_malloc_arrayN(nv, sizeof(output->vert_coords[0]), __func__));
2870  output->edges = static_cast<decltype(output->edges)>(
2871  MEM_malloc_arrayN(ne, sizeof(output->edges[0]), __func__));
2872  output->faces = static_cast<int *>(MEM_malloc_arrayN(tot_f_lens, sizeof(int), __func__));
2873  output->faces_start_table = static_cast<int *>(MEM_malloc_arrayN(nf, sizeof(int), __func__));
2874  output->faces_len_table = static_cast<int *>(MEM_malloc_arrayN(nf, sizeof(int), __func__));
2875  if (input->need_ids) {
2876  output->verts_orig = static_cast<int *>(MEM_malloc_arrayN(tot_v_orig, sizeof(int), __func__));
2877  output->verts_orig_start_table = static_cast<int *>(
2878  MEM_malloc_arrayN(nv, sizeof(int), __func__));
2879  output->verts_orig_len_table = static_cast<int *>(
2880  MEM_malloc_arrayN(nv, sizeof(int), __func__));
2881  output->edges_orig = static_cast<int *>(MEM_malloc_arrayN(tot_e_orig, sizeof(int), __func__));
2882  output->edges_orig_start_table = static_cast<int *>(
2883  MEM_malloc_arrayN(ne, sizeof(int), __func__));
2884  output->edges_orig_len_table = static_cast<int *>(
2885  MEM_malloc_arrayN(ne, sizeof(int), __func__));
2886  output->faces_orig = static_cast<int *>(MEM_malloc_arrayN(tot_f_orig, sizeof(int), __func__));
2887  output->faces_orig_start_table = static_cast<int *>(
2888  MEM_malloc_arrayN(nf, sizeof(int), __func__));
2889  output->faces_orig_len_table = static_cast<int *>(
2890  MEM_malloc_arrayN(nf, sizeof(int), __func__));
2891  }
2892  else {
2893  output->verts_orig = nullptr;
2894  output->verts_orig_start_table = nullptr;
2895  output->verts_orig_len_table = nullptr;
2896  output->edges_orig = nullptr;
2897  output->edges_orig_start_table = nullptr;
2898  output->edges_orig_len_table = nullptr;
2899  output->faces_orig = nullptr;
2900  output->faces_orig_start_table = nullptr;
2901  output->faces_orig_len_table = nullptr;
2902  }
2903 
2904  int v_orig_index = 0;
2905  for (int v = 0; v < nv; ++v) {
2906  output->vert_coords[v][0] = static_cast<float>(res.vert[v][0]);
2907  output->vert_coords[v][1] = static_cast<float>(res.vert[v][1]);
2908  if (input->need_ids) {
2909  int this_start = v_orig_index;
2910  output->verts_orig_start_table[v] = this_start;
2911  for (int j : res.vert_orig[v].index_range()) {
2912  output->verts_orig[v_orig_index++] = res.vert_orig[v][j];
2913  }
2914  output->verts_orig_len_table[v] = v_orig_index - this_start;
2915  }
2916  }
2917  int e_orig_index = 0;
2918  for (int e = 0; e < ne; ++e) {
2919  output->edges[e][0] = res.edge[e].first;
2920  output->edges[e][1] = res.edge[e].second;
2921  if (input->need_ids) {
2922  int this_start = e_orig_index;
2923  output->edges_orig_start_table[e] = this_start;
2924  for (int j : res.edge_orig[e].index_range()) {
2925  output->edges_orig[e_orig_index++] = res.edge_orig[e][j];
2926  }
2927  output->edges_orig_len_table[e] = e_orig_index - this_start;
2928  }
2929  }
2930  int f_orig_index = 0;
2931  int f_index = 0;
2932  for (int f = 0; f < nf; ++f) {
2933 
2934  output->faces_start_table[f] = f_index;
2935  int flen = res.face[f].size();
2936  output->faces_len_table[f] = flen;
2937  for (int j = 0; j < flen; ++j) {
2938  output->faces[f_index++] = res.face[f][j];
2939  }
2940  if (input->need_ids) {
2941  int this_start = f_orig_index;
2942  output->faces_orig_start_table[f] = this_start;
2943  for (int k : res.face_orig[f].index_range()) {
2944  output->faces_orig[f_orig_index++] = res.face_orig[f][k];
2945  }
2946  output->faces_orig_len_table[f] = f_orig_index - this_start;
2947  }
2948  }
2949  return output;
2950 }
2951 
2953 {
2954  MEM_freeN(result->vert_coords);
2955  MEM_freeN(result->edges);
2956  MEM_freeN(result->faces);
2957  MEM_freeN(result->faces_start_table);
2958  MEM_freeN(result->faces_len_table);
2959  if (result->verts_orig) {
2960  MEM_freeN(result->verts_orig);
2961  MEM_freeN(result->verts_orig_start_table);
2962  MEM_freeN(result->verts_orig_len_table);
2963  MEM_freeN(result->edges_orig);
2964  MEM_freeN(result->edges_orig_start_table);
2965  MEM_freeN(result->edges_orig_len_table);
2966  MEM_freeN(result->faces_orig);
2967  MEM_freeN(result->faces_orig_start_table);
2968  MEM_freeN(result->faces_orig_len_table);
2969  }
2970  MEM_freeN(result);
2971 }
#define BLI_assert(a)
Definition: BLI_assert.h:46
CDT_output_type
@ CDT_CONSTRAINTS_VALID_BMESH
@ CDT_CONSTRAINTS_VALID_BMESH_WITH_HOLES
@ CDT_CONSTRAINTS
@ CDT_INSIDE
@ CDT_FULL
@ CDT_INSIDE_WITH_HOLES
struct CDT_input CDT_input
MINLINE int max_ii(int a, int b)
Math vector functions needed specifically for mesh intersect and boolean.
#define UNUSED_VARS_NDEBUG(...)
#define UNUSED(x)
#define POINTER_AS_INT(i)
#define ELEM(...)
typedef double(DMatrix)[4][4]
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei height
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei width
_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
ATTR_WARN_UNUSED_RESULT const BMVert * v2
ATTR_WARN_UNUSED_RESULT const BMLoop * l
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
ATTR_WARN_UNUSED_RESULT const BMVert * v
void sort(btMatrix3x3 &U, btVector3 &sigma, btMatrix3x3 &V, int t)
Helper function of 3X3 SVD for sorting singular values.
SymEdge< T > * in
SymEdge< T > * out
CrossData & operator=(CrossData &&other) noexcept
~CrossData()=default
CrossData(const CrossData &other)
CrossData(CrossData &&other) noexcept
CrossData & operator=(const CrossData &other)
CDTVert< T > * vert
CrossData(T l, CDTVert< T > *v, SymEdge< T > *i, SymEdge< T > *o)
int orig_index
CDTVert< T > * v
int64_t size() const
Definition: BLI_set.hh:539
bool add(const Key &key)
Definition: BLI_set.hh:253
void clear()
Definition: BLI_set.hh:516
CDT_state(int input_verts_num, int input_edges_num, int input_faces_num, T epsilon, bool need_ids)
Definition: delaunay_2d.cc:871
const char * label
bool vert_touches_face(const CDTVert< T > *v, const CDTFace< T > *f)
Definition: delaunay_2d.cc:967
bool is_deleted_edge(const CDTEdge< T > *e)
Definition: delaunay_2d.cc:915
constexpr int inline_crossings_size
bool vert_left_of_symedge(CDTVert< T > *v, SymEdge< T > *se)
void dump_crossings(const Vector< CrossData< T >, inline_crossings_size > &crossings)
static void add_to_input_ids(blender::Set< int > &dst, int input_id)
Definition: delaunay_2d.cc:893
void add_face_ids(CDT_state< T > *cdt_state, SymEdge< T > *face_symedge, int face_id, int fedge_start, int fedge_end)
void fill_crossdata_for_intersect(const FatCo< T > &curco, const CDTVert< T > *v2, SymEdge< T > *t, CrossData< T > *cd, CrossData< T > *cd_next, const T epsilon)
void prepare_cdt_for_output(CDT_state< T > *cdt_state, const CDT_output_type output_type)
static int power_of_10_greater_equal_to(int x)
bool is_constrained_edge(const CDTEdge< T > *e)
Definition: delaunay_2d.cc:910
static int filtered_orient2d(const FatCo< T > &a, const FatCo< T > &b, const FatCo< T > &c)
bool get_next_crossing_from_vert(CDT_state< T > *cdt_state, CrossData< T > *cd, CrossData< T > *cd_next, const CDTVert< T > *v2)
void detect_holes(CDT_state< T > *cdt_state)
bool is_original_vert(const CDTVert< T > *v, CDT_state< T > *cdt)
Definition: delaunay_2d.cc:920
void remove_non_constraint_edges(CDT_state< T > *cdt_state)
int filtered_orient2d< double >(const FatCo< double > &a, const FatCo< double > &b, const FatCo< double > &c)
Definition: delaunay_2d.cc:668
bool exists_edge(const CDTVert< T > *v1, const CDTVert< T > *v2)
Definition: delaunay_2d.cc:959
static void re_delaunay_triangulate(CDTArrangement< T > *cdt, SymEdge< T > *se)
void dc_triangulate(CDTArrangement< T > *cdt, Array< SiteInfo< T >> &sites)
SymEdge< T > * find_symedge_with_face(const CDTVert< T > *v, const CDTFace< T > *f)
Definition: delaunay_2d.cc:944
#define SY(y)
void get_next_crossing_from_edge(CrossData< T > *cd, CrossData< T > *cd_next, const CDTVert< T > *v2, const T epsilon)
static bool in_line(const FatCo< T > &a, const FatCo< T > &b, const FatCo< T > &c)
void remove_non_constraint_edges_leave_valid_bmesh(CDT_state< T > *cdt_state)
void BLI_delaunay_2d_cdt_free(::CDT_result *result)
CDT_result< T > get_cdt_output(CDT_state< T > *cdt_state, const CDT_input< T > UNUSED(input), CDT_output_type output_type)
int tri_orient(const SymEdge< T > *t)
void add_edge_constraint(CDT_state< T > *cdt_state, CDTVert< T > *v1, CDTVert< T > *v2, int input_id, LinkNode **r_edges)
bool vert_right_of_symedge(CDTVert< T > *v, SymEdge< T > *se)
#define SX(x)
int filtered_incircle< double >(const FatCo< double > &a, const FatCo< double > &b, const FatCo< double > &c, const FatCo< double > &d)
Definition: delaunay_2d.cc:735
void add_edge_constraints(CDT_state< T > *cdt_state, const CDT_input< T > &input)
bool dc_tri_valid(SymEdge< T > *se, SymEdge< T > *basel, SymEdge< T > *basel_sym)
bool in_line< double >(const FatCo< double > &a, const FatCo< double > &b, const FatCo< double > &c)
Definition: delaunay_2d.cc:787
void fill_crossdata_for_through_vert(CDTVert< T > *v, SymEdge< T > *cd_out, CrossData< T > *cd, CrossData< T > *cd_next)
void remove_outer_edges_until_constraints(CDT_state< T > *cdt_state)
static bool id_range_in_list(const blender::Set< int > &id_list, int range_start, int range_end)
Definition: delaunay_2d.cc:883
void remove_faces_in_holes(CDT_state< T > *cdt_state)
void dissolve_symedge(CDT_state< T > *cdt_state, SymEdge< T > *se)
bool site_lexicographic_sort(const SiteInfo< T > &a, const SiteInfo< T > &b)
int add_face_constraints(CDT_state< T > *cdt_state, const CDT_input< T > &input, CDT_output_type output_type)
void find_site_merges(Array< SiteInfo< T >> &sites)
void initial_triangulation(CDTArrangement< T > *cdt)
void add_input_verts(CDT_state< T > *cdt_state, const CDT_input< T > &input)
SymEdge< T > * find_symedge_between_verts(const CDTVert< T > *v1, const CDTVert< T > *v2)
Definition: delaunay_2d.cc:929
blender::meshintersect::CDT_result< double > delaunay_2d_calc(const CDT_input< double > &input, CDT_output_type output_type)
CDT_result< T > delaunay_calc(const CDT_input< T > &input, CDT_output_type output_type)
static int filtered_incircle(const FatCo< T > &a, const FatCo< T > &b, const FatCo< T > &c, const FatCo< T > &d)
void dc_tri(CDTArrangement< T > *cdt, Array< SiteInfo< T >> &sites, int start, int end, SymEdge< T > **r_le, SymEdge< T > **r_re)
static void add_list_to_input_ids(blender::Set< int > &dst, const blender::Set< int > &src)
Definition: delaunay_2d.cc:898
bool is_border_edge(const CDTEdge< T > *e, const CDT_state< T > *cdt)
Definition: delaunay_2d.cc:905
::CDT_result * BLI_delaunay_2d_cdt_calc(const ::CDT_input *input, const CDT_output_type output_type)
SyclQueue void void * src
#define rot(x, k)
static float verts[][3]
int count
ccl_global KernelShaderEvalInput ccl_global float * output
ccl_global KernelShaderEvalInput * input
void *(* MEM_malloc_arrayN)(size_t len, size_t size, const char *str)
Definition: mallocn.c:34
void(* MEM_freeN)(void *vmemh)
Definition: mallocn.c:27
void *(* MEM_mallocN)(size_t len, const char *str)
Definition: mallocn.c:33
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
static ulong * next
#define T
static char faces[256]
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
void interpolate(const Span< T > src, const Span< int > indices, const Span< float > factors, MutableSpan< T > dst)
T distance_squared(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
T dot(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
T distance(const T &a, const T &b)
isect_result< vec_base< T, Size > > isect_seg_seg(const vec_base< T, Size > &v1, const vec_base< T, Size > &v2, const vec_base< T, Size > &v3, const vec_base< T, Size > &v4)
T abs(const T &a)
double math_to_double< double >(const double v)
Definition: delaunay_2d.cc:61
std::string sename(const SymEdge< T > *se)
Definition: delaunay_2d.cc:379
std::string vertname(const CDTVert< T > *v)
Definition: delaunay_2d.cc:363
double math_to_double(const T UNUSED(v))
Definition: delaunay_2d.cc:48
void cdt_draw(const std::string &label, const CDTArrangement< T > &cdt)
Definition: delaunay_2d.cc:475
SymEdge< T > * sym(const SymEdge< T > *se)
Definition: delaunay_2d.cc:99
std::string short_se_dump(const SymEdge< T > *se)
Definition: delaunay_2d.cc:404
static std::string trunc_ptr(const void *p)
Definition: delaunay_2d.cc:371
SymEdge< T > * prev(const SymEdge< T > *se)
Definition: delaunay_2d.cc:105
std::ostream & operator<<(std::ostream &stream, const FatCo< T > &co)
Definition: delaunay_2d.cc:177
double math_abs< double >(const double v)
Definition: delaunay_2d.cc:43
T math_abs(const T v)
Definition: delaunay_2d.cc:31
static void add_edge(const Mesh &mesh, const int old_edge_i, const int v1, const int v2, Vector< int > &new_to_old_edges_map, Vector< MEdge > &new_edges, Vector< int > &loop_edges)
static double epsilon
void parallel_for(IndexRange range, int64_t grain_size, const Function &function)
Definition: BLI_task.hh:51
vec_base< double, 2 > double2
int orient2d(const double2 &a, const double2 &b, const double2 &c)
int incircle(const double2 &a, const double2 &b, const double2 &c, const double2 &d)
std::string to_string(const T &n)
static const pxr::TfToken out("out", pxr::TfToken::Immortal)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
float co[3]
Definition: bmesh_class.h:87
struct BMEdge * e
Definition: bmesh_class.h:97
double len_squared
EdgeToSort(EdgeToSort &&other) noexcept
EdgeToSort & operator=(const EdgeToSort &other)
EdgeToSort()=default
~EdgeToSort()=default
EdgeToSort & operator=(EdgeToSort &&other)
EdgeToSort(const EdgeToSort &other)
CDTEdge< T > * e
LinkNode * list
Definition: BLI_linklist.h:34
void * link
Definition: BLI_linklist.h:24
CDTEdge< Arith_t > * add_diagonal(SymEdge< Arith_t > *s1, SymEdge< Arith_t > *s2)
Vector< CDTFace< Arith_t > * > faces
Definition: delaunay_2d.cc:239
CDTVert< Arith_t > * get_vert_resolve_merge(int i)
Definition: delaunay_2d.cc:307
void reserve(int verts_num, int edges_num, int faces_num)
CDTVert< Arith_t > * add_vert(const vec2< Arith_t > &pt)
void delete_edge(SymEdge< Arith_t > *se)
CDTEdge< Arith_t > * split_edge(SymEdge< Arith_t > *se, Arith_t lambda)
CDTEdge< Arith_t > * add_vert_to_symedge_edge(CDTVert< Arith_t > *v, SymEdge< Arith_t > *se)
Vector< CDTVert< Arith_t > * > verts
Definition: delaunay_2d.cc:235
Vector< CDTEdge< Arith_t > * > edges
Definition: delaunay_2d.cc:237
CDTEdge< Arith_t > * connect_separate_parts(SymEdge< Arith_t > *se1, SymEdge< Arith_t > *se2)
CDTEdge< Arith_t > * add_edge(CDTVert< Arith_t > *v1, CDTVert< Arith_t > *v2, CDTFace< Arith_t > *fleft, CDTFace< Arith_t > *fright)
blender::Set< int > input_ids
Definition: delaunay_2d.cc:205
blender::Set< int > input_ids
Definition: delaunay_2d.cc:218
SymEdge< Arith_t > * symedge
Definition: delaunay_2d.cc:214
blender::Set< int > input_ids
Definition: delaunay_2d.cc:189
CDTVert(const vec2< T > &pt)
FatCo(const vec2< double > &v)
SymEdge< Arith_t > * next
Definition: delaunay_2d.cc:83
CDTFace< Arith_t > * face
Definition: delaunay_2d.cc:91
CDTVert< Arith_t > * vert
Definition: delaunay_2d.cc:87
CDTEdge< Arith_t > * edge
Definition: delaunay_2d.cc:89
SymEdge< Arith_t > * rot
Definition: delaunay_2d.cc:85
static const char hex[17]
Definition: thumbs.c:156