Geogram  Version 1.9.1
A programming library of geometric algorithms
CDT_2d.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2000-2022 Inria
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright notice,
9  * this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright notice,
11  * this list of conditions and the following disclaimer in the documentation
12  * and/or other materials provided with the distribution.
13  * * Neither the name of the ALICE Project-Team nor the names of its
14  * contributors may be used to endorse or promote products derived from this
15  * software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  *
29  * Contact: Bruno Levy
30  *
31  * https://www.inria.fr/fr/bruno-levy
32  *
33  * Inria,
34  * Domaine de Voluceau,
35  * 78150 Le Chesnay - Rocquencourt
36  * FRANCE
37  *
38  */
39 
40 #ifndef GEOGRAM_DELAUNAY_CDT_2D
41 #define GEOGRAM_DELAUNAY_CDT_2D
42 
43 #include <geogram/basic/common.h>
44 #include <geogram/basic/geometry.h>
47 #include <geogram/mesh/index.h>
48 #include <functional>
49 
56 namespace GEO {
57 
62  struct CDT2d_ConstraintWalker;
63 
75  class GEOGRAM_API CDTBase2d {
76  public:
81 
85  virtual ~CDTBase2d();
86 
90  virtual void clear();
91 
98 
114  bool remove_internal_holes=false
115  );
116 
124  void set_delaunay(bool delaunay) {
125  delaunay_ = delaunay;
126  }
127 
131  index_t nT() const {
132  return T_.size()/3;
133  }
134 
138  index_t nv() const {
139  return nv_;
140  }
141 
145  index_t ncnstr() const {
146  return ncnstr_;
147  }
148 
155  index_t Tv(index_t t, index_t lv) const {
156  geo_debug_assert(t<nT());
157  geo_debug_assert(lv<3);
158  return T_[3*t+lv];
159  }
160 
168  geo_debug_assert(t<nT());
169  geo_debug_assert(v<nv());
170  return find_3(T_.data()+3*t, v);
171  }
172 
180  index_t Tadj(index_t t, index_t le) const {
181  geo_debug_assert(t<nT());
182  geo_debug_assert(le<3);
183  return Tadj_[3*t+le];
184  }
185 
194  geo_debug_assert(t1<nT());
195  geo_debug_assert(t2<nT());
196  return find_3(Tadj_.data()+3*t1, t2);
197  }
198 
205  index_t vT(index_t v) const {
206  geo_debug_assert(v < nv());
207  return v2T_[v];
208  }
209 
210 
238  geo_debug_assert(t < nT());
239  geo_debug_assert(le < 3);
240  return Tecnstr_first_[3*t+le];
241  }
242 
251  return ecnstr_next_[ecit];
252  }
253 
261  index_t edge_cnstr(index_t ecit) const {
262  return ecnstr_val_[ecit];
263  }
264 
277  index_t result = 0;
278  for(
279  index_t ecit = Tedge_cnstr_first(t,le);
280  ecit != index_t(-1);
281  ecit = edge_cnstr_next(ecit)
282  ) {
283  ++result;
284  }
285  return result;
286  }
287 
288 
293  virtual void save(const std::string& filename) const = 0;
294 
299  bool Tedge_is_Delaunay(index_t t, index_t le) const;
300 
301  protected:
302  virtual void begin_insert_transaction();
303  virtual void commit_insert_transaction();
304  virtual void rollback_insert_transaction();
305 
316 
326 
336  index_t v1, index_t v2, index_t v3, index_t v4
337  );
338 
344  void Tset_flag(index_t t, index_t flag) {
345  geo_debug_assert(t < nT());
346  geo_debug_assert(flag < 8);
347  Tflags_[t] |= Numeric::uint8(1u << flag);
348  }
349 
355  void Treset_flag(index_t t, index_t flag) {
356  geo_debug_assert(t < nT());
357  geo_debug_assert(flag < 8);
358  Tflags_[t] &= Numeric::uint8(~(1u << flag));
359  }
360 
368  bool Tflag_is_set(index_t t, index_t flag) {
369  geo_debug_assert(t < nT());
370  geo_debug_assert(flag < 8);
371  return ((Tflags_[t] & (1u << flag)) != 0);
372  }
373 
377  enum {
378  DLIST_S_ID=0,
379  DLIST_Q_ID=1,
380  DLIST_N_ID=2,
381  DLIST_NB=3
382  };
383 
387  enum {
388  T_MARKED_FLAG = DLIST_NB,
389  T_VISITED_FLAG = DLIST_NB+1
390  };
391 
398  bool Tis_in_list(index_t t) const {
399  return (
400  (Tflags_[t] &
401  Numeric::uint8((1 << DLIST_NB)-1)
402  ) != 0
403  );
404  }
405 
418  struct DList {
424  DList(CDTBase2d& cdt, index_t list_id) :
425  cdt_(cdt), list_id_(list_id),
426  back_(index_t(-1)), front_(index_t(-1)) {
427  geo_debug_assert(list_id < DLIST_NB);
428  }
429 
438  DList(CDTBase2d& cdt) :
439  cdt_(cdt), list_id_(index_t(-1)),
440  back_(index_t(-1)), front_(index_t(-1)) {
441  }
442 
447  void initialize(index_t list_id) {
448  geo_debug_assert(!initialized());
449  geo_debug_assert(list_id < DLIST_NB);
450  list_id_ = list_id;
451  }
452 
456  bool initialized() const {
457  return (list_id_ != index_t(-1));
458  }
459 
460  ~DList() {
461  if(initialized()) {
462  clear();
463  }
464  }
465 
466  bool empty() const {
467  geo_debug_assert(initialized());
469  (back_==index_t(-1))==(front_==index_t(-1))
470  );
471  return (back_==index_t(-1));
472  }
473 
474  bool contains(index_t t) const {
475  geo_debug_assert(initialized());
476  return cdt_.Tflag_is_set(t, list_id_);
477  }
478 
479  index_t front() const {
480  geo_debug_assert(initialized());
481  return front_;
482  }
483 
484  index_t back() const {
485  geo_debug_assert(initialized());
486  return back_;
487  }
488 
489  index_t next(index_t t) const {
490  geo_debug_assert(initialized());
491  geo_debug_assert(contains(t));
492  return cdt_.Tnext_[t];
493  }
494 
495  index_t prev(index_t t) const {
496  geo_debug_assert(initialized());
497  geo_debug_assert(contains(t));
498  return cdt_.Tprev_[t];
499  }
500 
501  void clear() {
502  for(index_t t=front_; t!=index_t(-1); t = cdt_.Tnext_[t]) {
503  cdt_.Treset_flag(t,list_id_);
504  }
505  back_ = index_t(-1);
506  front_ = index_t(-1);
507  }
508 
509  index_t size() const {
510  geo_debug_assert(initialized());
511  index_t result = 0;
512  for(index_t t=front(); t!=index_t(-1); t = next(t)) {
513  ++result;
514  }
515  return result;
516  }
517 
518  void push_back(index_t t) {
519  geo_debug_assert(initialized());
520  geo_debug_assert(!cdt_.Tis_in_list(t));
521  cdt_.Tset_flag(t,list_id_);
522  if(empty()) {
523  back_ = t;
524  front_ = t;
525  cdt_.Tnext_[t] = index_t(-1);
526  cdt_.Tprev_[t] = index_t(-1);
527  } else {
528  cdt_.Tnext_[t] = index_t(-1);
529  cdt_.Tnext_[back_] = t;
530  cdt_.Tprev_[t] = back_;
531  back_ = t;
532  }
533  }
534 
535  index_t pop_back() {
536  geo_debug_assert(initialized());
537  geo_debug_assert(!empty());
538  index_t t = back_;
539  back_ = cdt_.Tprev_[back_];
540  if(back_ == index_t(-1)) {
541  geo_debug_assert(front_ == t);
542  front_ = index_t(-1);
543  } else {
544  cdt_.Tnext_[back_] = index_t(-1);
545  }
546  geo_debug_assert(contains(t));
547  cdt_.Treset_flag(t,list_id_);
548  return t;
549  }
550 
551  void push_front(index_t t) {
552  geo_debug_assert(initialized());
553  geo_debug_assert(!cdt_.Tis_in_list(t));
554  cdt_.Tset_flag(t,list_id_);
555  if(empty()) {
556  back_ = t;
557  front_ = t;
558  cdt_.Tnext_[t] = index_t(-1);
559  cdt_.Tprev_[t] = index_t(-1);
560  } else {
561  cdt_.Tprev_[t] = index_t(-1);
562  cdt_.Tprev_[front_] = t;
563  cdt_.Tnext_[t] = front_;
564  front_ = t;
565  }
566  }
567 
568  index_t pop_front() {
569  geo_debug_assert(initialized());
570  geo_debug_assert(!empty());
571  index_t t = front_;
572  front_ = cdt_.Tnext_[front_];
573  if(front_ == index_t(-1)) {
574  geo_debug_assert(back_ == t);
575  back_ = index_t(-1);
576  } else {
577  cdt_.Tprev_[front_] = index_t(-1);
578  }
579  geo_debug_assert(contains(t));
580  cdt_.Treset_flag(t,list_id_);
581  return t;
582  }
583 
584  void remove(index_t t) {
585  geo_debug_assert(initialized());
586  if(t == front_) {
587  pop_front();
588  } else if(t == back_) {
589  pop_back();
590  } else {
591  geo_debug_assert(contains(t));
592  index_t t_prev = cdt_.Tprev_[t];
593  index_t t_next = cdt_.Tnext_[t];
594  cdt_.Tprev_[t_next] = t_prev;
595  cdt_.Tnext_[t_prev] = t_next;
596  cdt_.Treset_flag(t,list_id_);
597  }
598  }
599 
600  void show(std::ostream& out = std::cerr) const {
601  switch(list_id_) {
602  case DLIST_S_ID:
603  out << "S";
604  break;
605  case DLIST_Q_ID:
606  out << "Q";
607  break;
608  case DLIST_N_ID:
609  out << "N";
610  break;
611  case index_t(-1):
612  out << "<uninitialized list>";
613  break;
614  default:
615  out << "<unknown list id:" << list_id_ << ">";
616  break;
617  }
618  out << "=";
619  for(index_t t=front(); t!=index_t(-1); t = next(t)) {
620  out << t << ";";
621  }
622  out << std::endl;
623  }
624 
625  private:
626  CDTBase2d& cdt_;
627  index_t list_id_;
628  index_t back_;
629  index_t front_;
630  };
631 
640 
648  DList S(*this);
649  insert_vertex_in_edge(v,t,le,S);
650  }
651 
659 
675 
679  void walk_constraint_v(CDT2d_ConstraintWalker& W);
680 
684  void walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q);
685 
697 
707 
720 
728 
729 
740  void Tset(
741  index_t t,
742  index_t v1, index_t v2, index_t v3,
743  index_t adj1, index_t adj2, index_t adj3,
744  index_t e1cnstr = index_t(-1),
745  index_t e2cnstr = index_t(-1),
746  index_t e3cnstr = index_t(-1)
747  ) {
748  geo_debug_assert(t < nT());
749  geo_debug_assert(v1 < nv());
750  geo_debug_assert(v2 < nv());
751  geo_debug_assert(v3 < nv());
752  geo_debug_assert(adj1 < nT() || adj1 == index_t(-1));
753  geo_debug_assert(adj2 < nT() || adj2 == index_t(-1));
754  geo_debug_assert(adj3 < nT() || adj3 == index_t(-1));
755  geo_debug_assert(v1 != v2);
756  geo_debug_assert(v2 != v3);
757  geo_debug_assert(v3 != v1);
758  geo_debug_assert(adj1 != adj2 || adj1 == index_t(-1));
759  geo_debug_assert(adj2 != adj3 || adj2 == index_t(-1));
760  geo_debug_assert(adj3 != adj1 || adj3 == index_t(-1));
761  geo_debug_assert(orient2d(v1,v2,v3) != ZERO);
762  T_[3*t ] = v1;
763  T_[3*t+1] = v2;
764  T_[3*t+2] = v3;
765  Tadj_[3*t ] = adj1;
766  Tadj_[3*t+1] = adj2;
767  Tadj_[3*t+2] = adj3;
768  Tecnstr_first_[3*t] = e1cnstr;
769  Tecnstr_first_[3*t+1] = e2cnstr;
770  Tecnstr_first_[3*t+2] = e3cnstr;
771  v2T_[v1] = t;
772  v2T_[v2] = t;
773  v2T_[v3] = t;
774  }
775 
783  void Trot(index_t t, index_t lv) {
784  geo_debug_assert(t < nT());
785  geo_debug_assert(lv < 3);
786  if(lv != 0) {
787  index_t i = 3*t+lv;
788  index_t j = 3*t+((lv+1)%3);
789  index_t k = 3*t+((lv+2)%3);
790  Tset(
791  t,
792  T_[i], T_[j], T_[k],
793  Tadj_[i], Tadj_[j], Tadj_[k],
794  Tecnstr_first_[i], Tecnstr_first_[j], Tecnstr_first_[k]
795  );
796  }
797  }
798 
811  void swap_edge(index_t t1, bool swap_t1_t2=false);
812 
820  void Tadj_set(index_t t, index_t le, index_t adj) {
821  geo_debug_assert(t < nT());
822  geo_debug_assert(adj < nT());
823  geo_debug_assert(le < 3);
824  Tadj_[3*t+le] = adj;
825  }
826 
831  index_t Topp(index_t t, index_t e=0) const {
832  index_t t2 = Tadj(t,e);
833  if(t2 == index_t(-1)) {
834  return index_t(-1);
835  }
836  index_t e2 = Tadj_find(t2,t);
837  return Tv(t2,e2);
838  }
839 
853  index_t t1, index_t le1, index_t prev_t2_adj_e2
854  ) {
855  geo_debug_assert(t1 < nT());
856  geo_debug_assert(le1 < 3);
857  index_t t2 = Tadj(t1,le1);
858  if(t2 == index_t(-1)) {
859  return;
860  }
861  index_t le2 = Tadj_find(t2,prev_t2_adj_e2);
862  Tadj_set(t2,le2,t1);
863  Tset_edge_cnstr_first(t1,le1,Tedge_cnstr_first(t2,le2));
864  }
865 
871  index_t t = nT();
872  index_t nc = (t+1)*3; // new number of corners
873  T_.resize(nc, index_t(-1));
874  Tadj_.resize(nc, index_t(-1));
875  Tecnstr_first_.resize(nc, index_t(-1));
876  Tflags_.resize(t+1,0);
877  Tnext_.resize(t+1,index_t(-1));
878  Tprev_.resize(t+1,index_t(-1));
879  return t;
880  }
881 
890  index_t t, index_t le, index_t ecit
891  ) {
892  geo_debug_assert(t < nT());
893  geo_debug_assert(le < 3);
894  Tecnstr_first_[3*t+le] = ecit;
895  }
896 
904  index_t t, index_t le, index_t cnstr_id
905  ) {
906  geo_debug_assert(t < nT());
907  geo_debug_assert(le < 3);
908  // Check whether the edge is already constrained with the
909  // same constraint.
910  // TODO (if possible): understand how this can happen and
911  // remove this bloc of code that is not super elegant
912  // (it seems to be when we arrive at j and coming from a vertex
913  // traversed by the edge, both conditions make the constraint
914  // added to the traversed edge).
915  for(
916  index_t ecit = Tedge_cnstr_first(t,le);
917  ecit != index_t(-1);
918  ecit = edge_cnstr_next(ecit)
919  ) {
920  if(edge_cnstr(ecit) == cnstr_id) {
921  return;
922  }
923  }
924  ecnstr_val_.push_back(cnstr_id);
925  ecnstr_next_.push_back(Tedge_cnstr_first(t,le));
926  Tset_edge_cnstr_first(t,le, ecnstr_val_.size()-1);
927  }
928 
937  index_t t, index_t le, index_t cnstr_id
938  ) {
939  geo_debug_assert(t < nT());
940  geo_debug_assert(le < 3);
941 #ifdef GEO_DEBUG
942  index_t t_e_cnstr_first = Tedge_cnstr_first(t,le);
943 #endif
944  Tadd_edge_cnstr(t, le, cnstr_id);
945  index_t t2 = Tadj(t,le);
946  if(t2 != index_t(-1)) {
947  index_t le2 = Tadj_find(t2,t);
948  // Sanity check: make sure the two edges always share the
949  // same constraint list.
950  geo_debug_assert(Tedge_cnstr_first(t2,le2) == t_e_cnstr_first);
951  Tset_edge_cnstr_first(t2,le2,Tedge_cnstr_first(t,le));
952  }
953  }
954 
963  return (Tedge_cnstr_first(t,le) != index_t(-1));
964  }
965 
976  index_t v, std::function<bool(index_t t, index_t lv)> doit
977  ) {
978  index_t t = vT(v);
979  index_t lv = index_t(-1);
980  do {
981  lv = Tv_find(t,v);
982  if(doit(t,lv)) {
983  return;
984  }
985  t = Tadj(t, (lv+1)%3);
986  } while(t != vT(v) && t != index_t(-1));
987 
988  // We are done, this was an interior vertex
989  if(t != index_t(-1)) {
990  return;
991  }
992 
993  // It was a vertex on the border, so we need
994  // to traverse the triangle fan in the other
995  // direction until we reach the border again
996  t = vT(v);
997  lv = Tv_find(t,v);
998  t = Tadj(t, (lv+2)%3);
999  while(t != index_t(-1)) {
1000  lv = Tv_find(t,v);
1001  if(doit(t,lv)) {
1002  return;
1003  }
1004  t = Tadj(t, (lv+2)%3);
1005  }
1006  }
1007 
1008 
1020  index_t v, index_t hint = index_t(-1), Sign* orient = nullptr
1021  ) const;
1022 
1030  bool is_convex_quad(index_t t) const;
1031 
1037  virtual Sign orient2d(index_t i,index_t j,index_t k) const=0;
1038 
1050  virtual Sign incircle(index_t i,index_t j,index_t k,index_t l) const=0;
1051 
1072  index_t E1, index_t i, index_t j,
1073  index_t E2, index_t k, index_t l
1074  ) = 0;
1075 
1084  static inline index_t find_3(const index_t* T, index_t v) {
1085  // The following expression is 10% faster than using
1086  // if() statements. This uses the C++ norm, that
1087  // ensures that the 'true' boolean value converted to
1088  // an int is always 1. With most compilers, this avoids
1089  // generating branching instructions.
1090  // Thank to Laurent Alonso for this idea.
1091  index_t result = index_t( (T[1] == v) | ((T[2] == v) * 2) );
1092  // Sanity check, important if it was T[0], not explicitly
1093  // tested (detects input that does not meet the precondition).
1094  geo_debug_assert(T[result] == v);
1095  return result;
1096  }
1097 
1098  /*******************************************************************/
1099 
1109 
1110  /******************** Debugging ************************************/
1111 
1117  void Tcheck(index_t t) const {
1118  if(t == index_t(-1)) {
1119  return;
1120  }
1121  for(index_t e=0; e<3; ++e) {
1122  geo_assert(Tv(t,e) != Tv(t,(e+1)%3));
1123  if(Tadj(t,e) == index_t(-1)) {
1124  continue;
1125  }
1126  geo_assert(Tadj(t,e) != Tadj(t,(e+1)%3));
1127  index_t t2 = Tadj(t,e);
1128  index_t e2 = Tadj_find(t2,t);
1129  geo_assert(Tadj(t2,e2) == t);
1130  }
1131  }
1132 
1139  void debug_Tcheck(index_t t) const {
1140 #ifdef GEO_DEBUG
1141  Tcheck(t);
1142 #else
1143  geo_argused(t);
1144 #endif
1145  }
1146 
1151  void check_combinatorics() const {
1152  for(index_t t=0; t<nT(); ++t) {
1153  Tcheck(t);
1154  }
1155  }
1156 
1163 #ifdef GEO_DEBUG
1164  check_combinatorics();
1165 #endif
1166  }
1167 
1172  virtual void check_geometry() const;
1173 
1179  void debug_check_geometry() const {
1180 #ifdef GEO_DEBUG
1181  check_geometry();
1182 #endif
1183  }
1184 
1185 
1186  public:
1191  void check_consistency() const {
1192  check_combinatorics();
1193  check_geometry();
1194  }
1195 
1196  protected:
1203  debug_check_combinatorics();
1204  debug_check_geometry();
1205  }
1206 
1216  index_t u1, index_t u2, index_t v1, index_t v2
1217  ) const {
1218  if(orient2d(u1,u2,v1)*orient2d(u1,u2,v2) > 0) {
1219  return false;
1220  }
1221  return (orient2d(v1,v2,u1)*orient2d(v1,v2,u2) < 0);
1222  }
1223 
1235  index_t v1, index_t v2, index_t t, index_t le
1236  ) const {
1237  index_t u1 = Tv(t,(le + 1)%3);
1238  index_t u2 = Tv(t,(le + 2)%3);
1239  return segment_segment_intersect(u1,u2,v1,v2);
1240  }
1241 
1250  index_t v1, index_t v2, const DList& Q
1251  );
1252 
1253  typedef std::pair<index_t, index_t> Edge;
1254 
1260  index_t eT(Edge E) {
1261  index_t v1 = E.first;
1262  index_t v2 = E.second;
1263  index_t result = index_t(-1);
1264  for_each_T_around_v(
1265  v1, [&](index_t t, index_t lv)->bool {
1266  if(Tv(t, (lv+1)%3) == v2) {
1267  if(Tv(t, (lv+2)%3) != v1) {
1268  Trot(t, (lv+2)%3);
1269  }
1270  result = t;
1271  return true;
1272  } else if(Tv(t, (lv+1)%3) == v1) {
1273  if(Tv(t, (lv+2)%3) != v2) {
1274  Trot(t, (lv+2)%3);
1275  }
1276  result = t;
1277  return true;
1278  }
1279  return false;
1280  }
1281  );
1282  geo_debug_assert(result != index_t(-1));
1284  (Tv(result,1) == v1 && Tv(result,2) == v2) ||
1285  (Tv(result,1) == v2 && Tv(result,2) == v1)
1286  );
1287  return result;
1288  }
1289 
1295  index_t v, index_t hint = index_t(-1), Sign* orient = nullptr
1296  ) const;
1297 
1303  index_t i, index_t j, DList& Q, vector<Edge>& N
1304  );
1305 
1311 
1312  protected:
1313  index_t nv_;
1314  index_t ncnstr_;
1324  bool delaunay_;
1328  };
1329 
1330  /*****************************************************************/
1331 
1381  class GEOGRAM_API CDT2d: public CDTBase2d {
1382  public:
1383 
1384  CDT2d();
1385 
1386  ~CDT2d() override;
1387 
1391  void clear() override;
1392 
1400  const vec2& p1, const vec2& p2, const vec2& p3
1401  );
1402 
1411  const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4
1412  );
1413 
1414 
1422  double x1, double y1, double x2, double y2
1423  ) {
1424  create_enclosing_quad(
1425  vec2(x1,y1),
1426  vec2(x2,y1),
1427  vec2(x2,y2),
1428  vec2(x1,y2)
1429  );
1430  }
1431 
1440  index_t insert(const vec2& p, index_t hint = index_t(-1)) {
1441  debug_check_consistency();
1442  point_.push_back(p);
1443  index_t v = CDTBase2d::insert(point_.size()-1, hint);
1444  // If inserted point already existed in
1445  // triangulation, then nv() did not increase
1446  if(point_.size() > nv()) {
1447  point_.pop_back();
1448  }
1449  debug_check_consistency();
1450  return v;
1451  }
1452 
1477  void insert(
1478  index_t nb_points, const double* points,
1479  index_t* indices = nullptr,
1480  bool remove_unreferenced_vertices = false
1481  );
1482 
1486  void save(const std::string& filename) const override;
1487 
1493  const vec2 point(index_t v) const {
1494  geo_debug_assert(v < nv());
1495  return point_[v];
1496  }
1497 
1498  protected:
1502  Sign orient2d(index_t i, index_t j, index_t k) const override;
1503 
1507  Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
1508 
1513  index_t E1, index_t i, index_t j,
1514  index_t E2, index_t k, index_t l
1515  ) override;
1516 
1517  protected:
1518  vector<vec2> point_;
1519  };
1520 
1521  /*****************************************************************/
1522 
1537  class GEOGRAM_API ExactCDT2d : public CDTBase2d {
1538  public:
1539  typedef exact::vec2h ExactPoint;
1540 
1545 
1549  ~ExactCDT2d() override;
1550 
1554  void clear() override;
1555 
1568  const ExactPoint& p, index_t id=0, index_t hint = index_t(-1)
1569  );
1570 
1580  void insert_constraint(index_t v1, index_t v2, index_t operand_bits=0) {
1581  constraints_.push_back(bindex(v1,v2,bindex::KEEP_ORDER));
1582  cnstr_operand_bits_.push_back(operand_bits);
1584  }
1585 
1594  const ExactPoint& p1, const ExactPoint& p2,
1595  const ExactPoint& p3, const ExactPoint& p4
1596  );
1597 
1605  double x1, double y1, double x2, double y2
1606  ) {
1607  create_enclosing_quad(
1608  ExactPoint(vec2(x1,y1)),
1609  ExactPoint(vec2(x2,y1)),
1610  ExactPoint(vec2(x2,y2)),
1611  ExactPoint(vec2(x1,y2))
1612  );
1613  }
1614 
1620  const ExactPoint& point(index_t v) const {
1621  geo_debug_assert(v < nv());
1622  return point_[v];
1623  }
1624 
1631  geo_debug_assert(v < nv());
1632  return id_[v];
1633  }
1634 
1641  geo_debug_assert(v < nv());
1642  id_[v] = id;
1643  }
1644 
1660  const std::string& boolean_expression, bool mark_only=false
1661  );
1662 
1666  void save(const std::string& filename) const override;
1667 
1668  protected:
1669  void add_point(const ExactPoint& p, index_t id = index_t(-1));
1670  void begin_insert_transaction() override;
1671  void commit_insert_transaction() override;
1672  void rollback_insert_transaction() override;
1673 
1677  Sign orient2d(index_t i, index_t j, index_t k) const override;
1678 
1682  Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
1683 
1688  index_t E1, index_t i, index_t j,
1689  index_t E2, index_t k, index_t l
1690  ) override;
1691 
1692  protected:
1693  vector<ExactPoint> point_;
1694 #ifndef GEOGRAM_USE_EXACT_NT
1695  vector<double> length_;
1696 #endif
1697  vector<index_t> id_;
1698  vector<index_t> cnstr_operand_bits_;
1699  vector<index_t> facet_inclusion_bits_;
1700  mutable std::map<trindex, Sign> pred_cache_;
1701  bool use_pred_cache_insert_buffer_;
1702  mutable std::vector<std::pair<trindex, Sign>> pred_cache_insert_buffer_;
1703  vector<bindex> constraints_;
1704  };
1705 
1706  /*****************************************************************/
1707 
1708 }
1709 
1710 #endif
#define geo_assert(x)
Verifies that a condition is met.
Definition: assert.h:149
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
Constrained Delaunay triangulation.
Definition: CDT_2d.h:1381
index_t create_intersection(index_t E1, index_t i, index_t j, index_t E2, index_t k, index_t l) override
Given two segments that have an intersection, create the intersection.
void create_enclosing_quad(const vec2 &p1, const vec2 &p2, const vec2 &p3, const vec2 &p4)
Creates a first large enclosing quad.
index_t insert(const vec2 &p, index_t hint=index_t(-1))
Inserts a point.
Definition: CDT_2d.h:1440
Sign orient2d(index_t i, index_t j, index_t k) const override
Computes the orientation predicate in 2d.
Sign incircle(index_t i, index_t j, index_t k, index_t l) const override
Incircle predicate.
void insert(index_t nb_points, const double *points, index_t *indices=nullptr, bool remove_unreferenced_vertices=false)
Batch-inserts a set of point.
const vec2 point(index_t v) const
Gets a point by index.
Definition: CDT_2d.h:1493
void clear() override
Removes everything from this triangulation.
void create_enclosing_triangle(const vec2 &p1, const vec2 &p2, const vec2 &p3)
Creates a first large enclosing triangle.
void create_enclosing_rectangle(double x1, double y1, double x2, double y2)
Creates a first large enclosing rectangle.
Definition: CDT_2d.h:1421
void save(const std::string &filename) const override
Saves this CDT to a geogram mesh file.
Base class for constrained Delaunay triangulation.
Definition: CDT_2d.h:75
vector< index_t > ecnstr_val_
Definition: CDT_2d.h:1320
index_t insert(index_t v, index_t hint=index_t(-1))
Inserts a new point.
bool Tis_in_list(index_t t) const
Tests whether a triangle is in a DList.
Definition: CDT_2d.h:398
bool exact_incircle_
Definition: CDT_2d.h:1326
index_t locate(index_t v, index_t hint=index_t(-1), Sign *orient=nullptr) const
Locates a vertex.
index_t nv() const
Gets the number of vertices.
Definition: CDT_2d.h:138
void Trot(index_t t, index_t lv)
Rotates indices in triangle t in such a way that a given vertex becomes vertex 0.
Definition: CDT_2d.h:783
void debug_check_consistency() const
Checks both combinatorics and geometry in debug mode, ignored in release mode, aborts on unconsistenc...
Definition: CDT_2d.h:1202
CDTBase2d()
CDTBase2d constructor.
void Delaunayize_vertex_neighbors(index_t v, DList &S)
Restores Delaunay condition starting from the triangles incident to a given vertex.
void check_combinatorics() const
Consistency combinatorial check for all the triangles.
Definition: CDT_2d.h:1151
index_t Tedge_cnstr_first(index_t t, index_t le) const
Gets the constraint associated with an edge.
Definition: CDT_2d.h:237
static index_t find_3(const index_t *T, index_t v)
Finds the index of an integer in an array of three integers.
Definition: CDT_2d.h:1084
index_t find_intersected_edges(index_t i, index_t j, DList &Q)
Finds the edges intersected by a constraint.
void Delaunayize_new_edges(DList &N)
Restores Delaunay condition for a set of edges after inserting a constrained edge.
vector< uint8_t > Tflags_
Definition: CDT_2d.h:1318
vector< index_t > Tadj_
Definition: CDT_2d.h:1316
void insert_vertex_in_edge(index_t v, index_t t, index_t le, DList &S)
Inserts a vertex in an edge.
bool Tedge_is_constrained(index_t t, index_t le) const
Tests whether an edge is constrained.
Definition: CDT_2d.h:962
void insert_vertex_in_triangle(index_t v, index_t t, DList &S)
Inserts a vertex in a triangle.
virtual ~CDTBase2d()
CDTBase2d destructor.
vector< index_t > Tecnstr_first_
Definition: CDT_2d.h:1319
index_t Tv(index_t t, index_t lv) const
Gets a vertex of a triangle.
Definition: CDT_2d.h:155
index_t vT(index_t v) const
Gets a triangle incident to a given vertex.
Definition: CDT_2d.h:205
void create_enclosing_triangle(index_t v1, index_t v2, index_t v3)
Creates the combinatorics for a first large enclosing triangle.
index_t edge_cnstr_next(index_t ecit) const
Gets the successor of an edge constraint iterator.
Definition: CDT_2d.h:250
void Delaunayize_vertex_neighbors(index_t from_v)
Restores Delaunay condition starting from the triangles incident to a given vertex.
void remove_marked_triangles()
Removes all the triangles that have the flag T_MARKED_FLAG set.
bool segment_edge_intersect(index_t v1, index_t v2, index_t t, index_t le) const
Tests whether an edge triangle and a segment have a frank intersection.
Definition: CDT_2d.h:1234
virtual index_t create_intersection(index_t E1, index_t i, index_t j, index_t E2, index_t k, index_t l)=0
Given two segments that have an intersection, create the intersection.
index_t Topp(index_t t, index_t e=0) const
Gets the neighboring triangle vertex opposite to a given vertex.
Definition: CDT_2d.h:831
virtual Sign orient2d(index_t i, index_t j, index_t k) const =0
Orientation predicate.
bool Tflag_is_set(index_t t, index_t flag)
Tests a triangle flag.
Definition: CDT_2d.h:368
index_t Tnew()
Creates a new triangle.
Definition: CDT_2d.h:870
virtual Sign incircle(index_t i, index_t j, index_t k, index_t l) const =0
Incircle predicate.
void Tadd_edge_cnstr_with_neighbor(index_t t, index_t le, index_t cnstr_id)
Adds a constraint to a triangle edge and to the neighboring edge if it exists.
Definition: CDT_2d.h:936
void walk_constraint_t(CDT2d_ConstraintWalker &W, DList &Q)
Used by find_intersected_edges()
vector< index_t > Tnext_
Definition: CDT_2d.h:1322
void Tset_flag(index_t t, index_t flag)
Sets a triangle flag.
Definition: CDT_2d.h:344
bool exact_intersections_
Definition: CDT_2d.h:1327
void remove_external_triangles(bool remove_internal_holes=false)
Recursively removes all the triangles adjacent to the border, and keeps what's surrounded by constrai...
void debug_check_geometry() const
Consistency geometrical check for all the triangles in debug mode, ignored in release mode.
Definition: CDT_2d.h:1179
void check_consistency() const
Checks both combinatorics and geometry, aborts on unconsistency.
Definition: CDT_2d.h:1191
void Tadj_back_connect(index_t t1, index_t le1, index_t prev_t2_adj_e2)
After having changed connections from triangle to a neighbor, creates connections from neighbor to tr...
Definition: CDT_2d.h:852
virtual void save(const std::string &filename) const =0
Saves this CDT to a geogram mesh file.
void swap_edge(index_t t1, bool swap_t1_t2=false)
Swaps an edge.
void create_enclosing_quad(index_t v1, index_t v2, index_t v3, index_t v4)
Creates the combinatorics for a first large enclosing quad.
virtual void clear()
Removes everything from this triangulation.
index_t nT() const
Gets the number of triangles.
Definition: CDT_2d.h:131
void debug_Tcheck(index_t t) const
Consistency check for a triangle in debug mode, ignored in release mode.
Definition: CDT_2d.h:1139
void walk_constraint_v(CDT2d_ConstraintWalker &W)
Used by find_intersected_edges()
index_t eT(Edge E)
Gets a triangle incident a a given edge.
Definition: CDT_2d.h:1260
Sign orient_012_
Definition: CDT_2d.h:1325
void insert_vertex_in_edge(index_t v, index_t t, index_t le)
Inserts a vertex in an edge.
Definition: CDT_2d.h:647
void for_each_T_around_v(index_t v, std::function< bool(index_t t, index_t lv)> doit)
Calls a user-defined function for each triangle around a vertex.
Definition: CDT_2d.h:975
void Tadj_set(index_t t, index_t le, index_t adj)
Sets a triangle adjacency relation.
Definition: CDT_2d.h:820
virtual void check_geometry() const
Consistency geometrical check for all the triangles.
index_t edge_cnstr(index_t ecit) const
Gets an edge constraint from an edge constraint iterator.
Definition: CDT_2d.h:261
index_t Tadj_find(index_t t1, index_t t2) const
Finds the edge accross which a triangle is adjacent to another one.
Definition: CDT_2d.h:193
void debug_check_combinatorics() const
Consistency combinatorial check for all the triangles in debug mode, ignored in release mode.
Definition: CDT_2d.h:1162
void Tadd_edge_cnstr(index_t t, index_t le, index_t cnstr_id)
Adds a constraint to a triangle edge.
Definition: CDT_2d.h:903
void Tset_edge_cnstr_first(index_t t, index_t le, index_t ecit)
Sets the constraints list associated with an edge.
Definition: CDT_2d.h:889
index_t Tv_find(index_t t, index_t v) const
Finds the local index of a vertex in a triangle.
Definition: CDT_2d.h:167
vector< index_t > T_
Definition: CDT_2d.h:1315
vector< index_t > v2T_
Definition: CDT_2d.h:1317
void constrain_edges_naive(index_t i, index_t j, DList &Q, vector< Edge > &N)
Simpler version of constrain_edges() kept for reference.
void Tset(index_t t, index_t v1, index_t v2, index_t v3, index_t adj1, index_t adj2, index_t adj3, index_t e1cnstr=index_t(-1), index_t e2cnstr=index_t(-1), index_t e3cnstr=index_t(-1))
Sets all the combinatorial information of a triangle and edge flags.
Definition: CDT_2d.h:740
void constrain_edges(index_t i, index_t j, DList &Q, DList &N)
Constrains an edge by iteratively flipping the intersected edges.
vector< index_t > Tprev_
Definition: CDT_2d.h:1323
void Tcheck(index_t t) const
Consistency check for a triangle.
Definition: CDT_2d.h:1117
index_t ncnstr() const
Gets the number of constraints.
Definition: CDT_2d.h:145
void Treset_flag(index_t t, index_t flag)
Resets a triangle flag.
Definition: CDT_2d.h:355
index_t locate_naive(index_t v, index_t hint=index_t(-1), Sign *orient=nullptr) const
Simpler version of locate() kept for reference.
bool is_convex_quad(index_t t) const
Tests whether triange t and its neighbor accross edge 0 form a strictly convex quad.
void insert_constraint(index_t i, index_t j)
Inserts a constraint.
bool segment_segment_intersect(index_t u1, index_t u2, index_t v1, index_t v2) const
Tests whether two segments have a frank intersection.
Definition: CDT_2d.h:1215
void Delaunayize_new_edges_naive(vector< Edge > &N)
Simpler version of Delaunayize_new_edges() that uses a vector instead of a DList, kept for reference.
index_t Tedge_cnstr_nb(index_t t, index_t le) const
Gets the number of constraints associated with a triange edge.
Definition: CDT_2d.h:276
bool delaunay_
Definition: CDT_2d.h:1324
index_t Tadj(index_t t, index_t le) const
Gets a triangle adjacent to a triangle.
Definition: CDT_2d.h:180
void set_delaunay(bool delaunay)
Specifies whether a constrained Delaunay triangulation should be constructed, or just a plain constra...
Definition: CDT_2d.h:124
void check_edge_intersections(index_t v1, index_t v2, const DList &Q)
Checks that the edges stored in a DList exactly correspond to all edge intersections between a segmen...
bool Tedge_is_Delaunay(index_t t, index_t le) const
Tests whether a triangle edge is Delaunay.
vector< index_t > ecnstr_next_
Definition: CDT_2d.h:1321
Constrained Delaunay Triangulation with vertices that are exact points. Can be used to implement 2D C...
Definition: CDT_2d.h:1537
~ExactCDT2d() override
ExactCDT2d destructor.
ExactCDT2d()
ExactCDT2d constructor.
void create_enclosing_quad(const ExactPoint &p1, const ExactPoint &p2, const ExactPoint &p3, const ExactPoint &p4)
Creates a first large enclosing quad.
void classify_triangles(const std::string &boolean_expression, bool mark_only=false)
Used by 2D CSG operations, discards triangles according to a boolean operation.
Sign orient2d(index_t i, index_t j, index_t k) const override
Computes the orientation predicate in 2d.
void set_vertex_id(index_t v, index_t id)
Sets a vertex id by vertex index.
Definition: CDT_2d.h:1640
const ExactPoint & point(index_t v) const
Gets a point by vertex index.
Definition: CDT_2d.h:1620
index_t insert(const ExactPoint &p, index_t id=0, index_t hint=index_t(-1))
Inserts a point.
Sign incircle(index_t i, index_t j, index_t k, index_t l) const override
Incircle predicate.
index_t create_intersection(index_t E1, index_t i, index_t j, index_t E2, index_t k, index_t l) override
Given two segments that have an intersection, create the intersection.
index_t vertex_id(index_t v) const
Gets a vertex id by vertex index.
Definition: CDT_2d.h:1630
void clear() override
Removes everything from this triangulation.
void create_enclosing_rectangle(double x1, double y1, double x2, double y2)
Creates a first large enclosing rectangle.
Definition: CDT_2d.h:1604
void save(const std::string &filename) const override
void insert_constraint(index_t v1, index_t v2, index_t operand_bits=0)
Inserts a constraint.
Definition: CDT_2d.h:1580
2d vector with homogeneous coordinates
Definition: vechg.h:60
Vector with aligned memory allocation.
Definition: memory.h:635
Exact predicates and constructs.
Common include file, providing basic definitions. Should be included before anything else by all head...
Geometric functions in 2d and 3d.
Classes for managing tuples of indices.
basic_bindex< index_t > bindex
A basic_bindex made of 2 unsigned integers.
Definition: index.h:205
void clear(void *addr, size_t size)
Clears a memory block.
Definition: memory.h:116
uint8_t uint8
Definition: numeric.h:135
Global Vorpaline namespace.
Definition: basic.h:55
void geo_argused(const T &)
Suppresses compiler warnings about unused parameters.
Definition: argused.h:60
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
Sign
Integer constants that represent the sign of a value.
Definition: numeric.h:68
@ ZERO
Definition: numeric.h:72
vecng< 2, Numeric::float64 > vec2
Represents points and vectors in 2d.
Definition: geometry.h:59
Filtered exact predicates for restricted Voronoi diagrams.
Doubly connected triangle list.
Definition: CDT_2d.h:418
bool initialized() const
Tests whether a DList is initialized.
Definition: CDT_2d.h:456
void initialize(index_t list_id)
Initializes a list.
Definition: CDT_2d.h:447
DList(CDTBase2d &cdt, index_t list_id)
Constructs an empty DList.
Definition: CDT_2d.h:424
DList(CDTBase2d &cdt)
Creates an uninitialized DList.
Definition: CDT_2d.h:438