Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
delaunay_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_DELAUNAY_2D
41 #define GEOGRAM_DELAUNAY_DELAUNAY_2D
42 
43 #include <geogram/basic/common.h>
46 #include <geogram/basic/geometry.h>
47 
48 #include <stack>
49 
55 namespace GEO {
56 
106  class GEOGRAM_API Delaunay2d : public Delaunay {
107  public:
121  Delaunay2d(coord_index_t dimension = 2);
122 
127  index_t nb_vertices, const double* vertices
128  ) override;
129 
133  index_t nearest_vertex(const double* p) const override;
134 
144  bool has_empty_cells() const {
145  return has_empty_cells_;
146  }
147 
148 
155  void abort_if_empty_cell(bool x) {
156  abort_if_empty_cell_ = x;
157  }
158 
159  protected:
160 
167  static const index_t NO_TRIANGLE = index_t(-1);
168 
181  index_t& iv0, index_t& iv1, index_t& iv2
182  );
183 
204  const double* p, index_t hint = NO_TRIANGLE,
205  bool thread_safe = false,
206  Sign* orient = nullptr
207  ) const;
208 
227  const double* p, index_t hint, index_t max_iter
228  ) const;
229 
238  index_t insert(index_t v, index_t hint = NO_TRIANGLE);
239 
264  index_t v,
265  index_t t, const Sign* orient,
266  index_t& t_bndry, index_t& e_bndry,
267  index_t& first, index_t& last
268  );
269 
286  const double* p, index_t t,
287  index_t& t_bndry, index_t& e_bndry,
288  index_t& first, index_t& last
289  );
290 
291 
307  index_t v,
308  index_t t_bndry, index_t e_bndry
309  );
310 
311  // _________ Combinatorics - new and delete _________________________
312 
320  index_t max_t() const {
321  return cell_to_v_store_.size() / 3;
322  }
323 
324 
340  static const index_t NOT_IN_LIST = index_t(~0);
341 
357  static const index_t NOT_IN_LIST_BIT = index_t(1u << 31);
358 
364  static const index_t END_OF_LIST = ~(NOT_IN_LIST_BIT);
365 
366 
384  bool triangle_is_in_list(index_t t) const {
385  geo_debug_assert(t < max_t());
386  return (cell_next_[t] & NOT_IN_LIST_BIT) == 0;
387  }
388 
400  geo_debug_assert(t < max_t());
401  geo_debug_assert(triangle_is_in_list(t));
402  return cell_next_[t];
403  }
404 
415  void add_triangle_to_list(index_t t, index_t& first, index_t& last) {
416  geo_debug_assert(t < max_t());
417  geo_debug_assert(!triangle_is_in_list(t));
418  if(last == END_OF_LIST) {
419  geo_debug_assert(first == END_OF_LIST);
420  first = last = t;
421  cell_next_[t] = END_OF_LIST;
422  } else {
423  cell_next_[t] = first;
424  first = t;
425  }
426  }
427 
438  geo_debug_assert(t < max_t());
439  geo_debug_assert(triangle_is_in_list(t));
440  cell_next_[t] = NOT_IN_LIST;
441  }
442 
449  static const signed_index_t VERTEX_AT_INFINITY = -1;
450 
461  bool triangle_is_finite(index_t t) const {
462  return
463  cell_to_v_store_[3 * t] >= 0 &&
464  cell_to_v_store_[3 * t + 1] >= 0 &&
465  cell_to_v_store_[3 * t + 2] >= 0 ;
466  }
467 
479  bool triangle_is_real(index_t t) const {
480  return !triangle_is_free(t) && triangle_is_finite(t);
481  }
482 
492  bool triangle_is_virtual(index_t t) const {
493  return
494  !triangle_is_free(t) && (
495  cell_to_v_store_[3 * t] == VERTEX_AT_INFINITY ||
496  cell_to_v_store_[3 * t + 1] == VERTEX_AT_INFINITY ||
497  cell_to_v_store_[3 * t + 2] == VERTEX_AT_INFINITY
498  );
499  }
500 
511  bool triangle_is_free(index_t t) const {
512  return triangle_is_in_list(t);
513  }
514 
523  index_t result;
524  if(first_free_ == END_OF_LIST) {
525  cell_to_v_store_.resize(cell_to_v_store_.size() + 3, -1);
526  cell_to_cell_store_.resize(cell_to_cell_store_.size() + 3, -1);
527  // index_t(NOT_IN_LIST) is necessary, else with
528  // NOT_IN_LIST alone the compiler tries to generate a
529  // reference to NOT_IN_LIST resulting in a link error.
530  cell_next_.push_back(index_t(NOT_IN_LIST));
531  result = max_t() - 1;
532  } else {
533  result = first_free_;
534  first_free_ = triangle_next(first_free_);
535  remove_triangle_from_list(result);
536  }
537 
538  cell_to_cell_store_[3 * result] = -1;
539  cell_to_cell_store_[3 * result + 1] = -1;
540  cell_to_cell_store_[3 * result + 2] = -1;
541 
542  return result;
543  }
544 
558  signed_index_t v3
559  ) {
560  index_t result = new_triangle();
561  cell_to_v_store_[3 * result] = v1;
562  cell_to_v_store_[3 * result + 1] = v2;
563  cell_to_v_store_[3 * result + 2] = v3;
564  return result;
565  }
566 
575  cur_stamp_ = (stamp | NOT_IN_LIST_BIT);
576  }
577 
593  bool triangle_is_marked(index_t t) const {
594  return cell_next_[t] == cur_stamp_;
595  }
596 
611  cell_next_[t] = cur_stamp_;
612  }
613 
614  // _________ Combinatorics ___________________________________
615 
633  geo_debug_assert(e < 3);
634  geo_debug_assert(v < 2);
635  return index_t(triangle_edge_vertex_[e][v]);
636  }
637 
646  geo_debug_assert(t < max_t());
647  geo_debug_assert(lv < 3);
648  return cell_to_v_store_[3 * t + lv];
649  }
650 
659  geo_debug_assert(t < max_t());
660  // Find local index of v in triangle t vertices.
661  const signed_index_t* T = &(cell_to_v_store_[3 * t]);
662  return find_3(T,v);
663  }
664 
665 
674  geo_debug_assert(t < max_t());
675  geo_debug_assert(lv < 3);
676  geo_debug_assert(cell_to_v_store_[3 * t + lv] != -1);
677  return index_t(cell_to_v_store_[3 * t + lv]);
678  }
679 
687  geo_debug_assert(t < max_t());
688  geo_debug_assert(lv < 3);
689  cell_to_v_store_[3 * t + lv] = v;
690  }
691 
699  geo_debug_assert(t < max_t());
700  geo_debug_assert(le < 3);
701  signed_index_t result = cell_to_cell_store_[3 * t + le];
702  return result;
703  }
704 
713  geo_debug_assert(t1 < max_t());
714  geo_debug_assert(t2 < max_t());
715  geo_debug_assert(le1 < 3);
716  geo_debug_assert(t1 != t2);
717  cell_to_cell_store_[3 * t1 + le1] = signed_index_t(t2);
718  }
719 
729  index_t t1, index_t t2_in
730  ) const {
731  geo_debug_assert(t1 < max_t());
732  geo_debug_assert(t2_in < max_t());
733  geo_debug_assert(t1 != t2_in);
734 
735  signed_index_t t2 = signed_index_t(t2_in);
736 
737  // Find local index of t2 in triangle t1 adajcent tets.
738  const signed_index_t* T = &(cell_to_cell_store_[3 * t1]);
739  index_t result = find_3(T,t2);
740 
741  // Sanity check: make sure that t1 is adjacent to t2
742  // only once!
743  geo_debug_assert(triangle_adjacent(t1,(result+1)%3) != t2);
744  geo_debug_assert(triangle_adjacent(t1,(result+2)%3) != t2);
745  return result;
746  }
747 
748 
749 
761  void set_tet(
762  index_t t,
764  signed_index_t v2,
765  index_t a0, index_t a1, index_t a2
766  ) {
767  geo_debug_assert(t < max_t());
768  cell_to_v_store_[3 * t] = v0;
769  cell_to_v_store_[3 * t + 1] = v1;
770  cell_to_v_store_[3 * t + 2] = v2;
771  cell_to_cell_store_[3 * t] = signed_index_t(a0);
772  cell_to_cell_store_[3 * t + 1] = signed_index_t(a1);
773  cell_to_cell_store_[3 * t + 2] = signed_index_t(a2);
774  }
775 
776  // _________ Predicates _____________________________________________
777 
791  bool triangle_is_conflict(index_t t, const double* p) const {
792 
793  // Lookup triangle vertices
794  const double* pv[3];
795  for(index_t i=0; i<3; ++i) {
796  signed_index_t v = triangle_vertex(t,i);
797  pv[i] = (v == -1) ? nullptr : vertex_ptr(index_t(v));
798  }
799 
800  // Check for virtual triangles (then in_circle()
801  // is replaced with orient2d())
802  for(index_t le = 0; le < 3; ++le) {
803 
804  if(pv[le] == nullptr) {
805 
806  // Facet of a virtual triangle opposite to
807  // infinite vertex corresponds to
808  // the triangle on the convex hull of the points.
809  // Orientation is obtained by replacing vertex lf
810  // with p.
811  pv[le] = p;
812  Sign sign = PCK::orient_2d(pv[0],pv[1],pv[2]);
813 
814  if(sign > 0) {
815  return true;
816  }
817 
818  if(sign < 0) {
819  return false;
820  }
821 
822  // If sign is zero, we check the real triangle
823  // adjacent to the facet on the convex hull.
824  geo_debug_assert(triangle_adjacent(t, le) >= 0);
825  index_t t2 = index_t(triangle_adjacent(t, le));
826  geo_debug_assert(!triangle_is_virtual(t2));
827 
828  // If t2 is already chained in the conflict list,
829  // then it is conflict
830  if(triangle_is_in_list(t2)) {
831  return true;
832  }
833 
834  // If t2 is marked, then it is not in conflict.
835  if(triangle_is_marked(t2)) {
836  return false;
837  }
838 
839  return triangle_is_conflict(t2, p);
840  }
841  }
842 
843  // If the triangle is a finite one, it is in conflict
844  // if its circumscribed sphere contains the point (this is
845  // the standard case).
846 
847  if(weighted_) {
848  double h0 = heights_[finite_triangle_vertex(t, 0)];
849  double h1 = heights_[finite_triangle_vertex(t, 1)];
850  double h2 = heights_[finite_triangle_vertex(t, 2)];
851  index_t pindex = index_t(
852  (p - vertex_ptr(0)) / int(vertex_stride_)
853  );
854  double h = heights_[pindex];
855  return (PCK::orient_2dlifted_SOS(
856  pv[0],pv[1],pv[2],p,h0,h1,h2,h
857  ) > 0) ;
858  }
859 
860  return (PCK::in_circle_2d_SOS(pv[0], pv[1], pv[2], p) > 0);
861  }
862 
863  protected:
864 
873  static inline index_t find_3(
874  const signed_index_t* T, signed_index_t v
875  ) {
876  // The following expression is 10% faster than using
877  // if() statements. This uses the C++ norm, that
878  // ensures that the 'true' boolean value converted to
879  // an int is always 1. With most compilers, this avoids
880  // generating branching instructions.
881  // Thank to Laurent Alonso for this idea.
882  index_t result = index_t( (T[1] == v) | ((T[2] == v) * 2) );
883  // Sanity check, important if it was T[0], not explicitly
884  // tested (detects input that does not meet the precondition).
885  geo_debug_assert(T[result] == v);
886  return result;
887  }
888 
892  ~Delaunay2d() override;
893 
898  void show_triangle(index_t t) const;
899 
907 
913  void show_list(
914  index_t first, const std::string& list_name
915  ) const;
916 
920  void check_combinatorics(bool verbose = false) const;
921 
925  void check_geometry(bool verbose = false) const;
926 
927  private:
928  vector<signed_index_t> cell_to_v_store_;
929  vector<signed_index_t> cell_to_cell_store_;
930  vector<index_t> cell_next_;
931  vector<index_t> reorder_;
932  index_t cur_stamp_; // used for marking
933  index_t first_free_;
934  bool weighted_;
935  vector<double> heights_; // only used in weighted mode
936 
940  bool debug_mode_;
941 
945  bool verbose_debug_mode_;
946 
950  bool benchmark_mode_;
951 
960  static char triangle_edge_vertex_[3][2];
961 
965  std::stack<index_t> S_;
966 
970  bool has_empty_cells_;
971 
976  bool abort_if_empty_cell_;
977  };
978 
979  /************************************************************************/
980 
992  class GEOGRAM_API RegularWeightedDelaunay2d : public Delaunay2d {
993  public:
1004 
1005  protected:
1010  };
1011 }
1012 
1013 #endif
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
Implementation of Delaunay in 2d.
Definition: delaunay_2d.h:106
void find_conflict_zone_iterative(const double *p, index_t t, index_t &t_bndry, index_t &e_bndry, index_t &first, index_t &last)
This function is used to implement find_conflict_zone.
bool triangle_is_finite(index_t t) const
Tests whether a given triangle is a finite one.
Definition: delaunay_2d.h:461
static index_t triangle_edge_vertex(index_t e, index_t v)
Returns the local index of a vertex by edge and by local vertex index in the edge.
Definition: delaunay_2d.h:632
bool triangle_is_real(index_t t) const
Tests whether a triangle is a real one.
Definition: delaunay_2d.h:479
void set_triangle_vertex(index_t t, index_t lv, signed_index_t v)
Sets a triangle-to-vertex adjacency.
Definition: delaunay_2d.h:686
void set_vertices(index_t nb_vertices, const double *vertices) override
Sets the vertices of this Delaunay, and recomputes the cells.
bool triangle_is_in_list(index_t t) const
Tests whether a triangle belongs to a linked list.
Definition: delaunay_2d.h:384
index_t triangle_next(index_t t) const
Gets the index of a successor of a triangle.
Definition: delaunay_2d.h:399
signed_index_t triangle_vertex(index_t t, index_t lv) const
Gets the index of a vertex of a triangle.
Definition: delaunay_2d.h:645
void set_tet(index_t t, signed_index_t v0, signed_index_t v1, signed_index_t v2, index_t a0, index_t a1, index_t a2)
Sets the vertices and adjacent triangles of a triangle.
Definition: delaunay_2d.h:761
void set_triangle_mark_stamp(index_t stamp)
Generates a unique stamp for marking tets.
Definition: delaunay_2d.h:574
index_t new_triangle(signed_index_t v1, signed_index_t v2, signed_index_t v3)
Creates a new triangle.
Definition: delaunay_2d.h:556
signed_index_t triangle_adjacent(index_t t, index_t le) const
Gets the index of a triangle adjacent to another one.
Definition: delaunay_2d.h:698
void show_list(index_t first, const std::string &list_name) const
For debugging purposes, displays a triangle.
index_t locate_inexact(const double *p, index_t hint, index_t max_iter) const
Finds the triangle that (approximately) contains a point using inexact predicates.
index_t max_t() const
Maximum valid index for a triangle.
Definition: delaunay_2d.h:320
Delaunay2d(coord_index_t dimension=2)
Constructs a new Delaunay2d.
~Delaunay2d() override
Delaunay2d destructor.
void show_triangle(index_t t) const
For debugging purposes, displays a triangle.
void set_triangle_adjacent(index_t t1, index_t le1, index_t t2)
Sets a triangle-to-triangle adjacency.
Definition: delaunay_2d.h:712
void check_geometry(bool verbose=false) const
For debugging purposes, test some geometrical properties.
bool triangle_is_virtual(index_t t) const
Tests whether a triangle is a virtual one.
Definition: delaunay_2d.h:492
void abort_if_empty_cell(bool x)
Specifies behavior if an empty cell is detected.
Definition: delaunay_2d.h:155
index_t stellate_conflict_zone(index_t v, index_t t_bndry, index_t e_bndry)
Creates a star of triangles filling the conflict zone.
bool create_first_triangle(index_t &iv0, index_t &iv1, index_t &iv2)
Finds in the pointset a set of three non-colinear points.
index_t find_triangle_vertex(index_t t, signed_index_t v) const
Finds the index of the vertex in a triangle.
Definition: delaunay_2d.h:658
index_t insert(index_t v, index_t hint=NO_TRIANGLE)
Inserts a point in the triangulation.
void find_conflict_zone(index_t v, index_t t, const Sign *orient, index_t &t_bndry, index_t &e_bndry, index_t &first, index_t &last)
Determines the list of triangles in conflict with a given point.
index_t finite_triangle_vertex(index_t t, index_t lv) const
Gets the index of a vertex of a triangle.
Definition: delaunay_2d.h:673
index_t find_triangle_adjacent(index_t t1, index_t t2_in) const
Finds the index of the edge accros which t1 is adjacent to t2_in.
Definition: delaunay_2d.h:728
void remove_triangle_from_list(index_t t)
Removes a triangle from the linked list it belongs to.
Definition: delaunay_2d.h:437
void mark_triangle(index_t t)
Marks a triangle.
Definition: delaunay_2d.h:610
index_t locate(const double *p, index_t hint=NO_TRIANGLE, bool thread_safe=false, Sign *orient=nullptr) const
Finds the triangle that contains a point.
bool has_empty_cells() const
Tests whether the Laguerre diagram has empty cells.
Definition: delaunay_2d.h:144
bool triangle_is_marked(index_t t) const
Tests whether a triangle is marked.
Definition: delaunay_2d.h:593
void add_triangle_to_list(index_t t, index_t &first, index_t &last)
Adds a triangle to a linked list.
Definition: delaunay_2d.h:415
bool triangle_is_free(index_t t) const
Tests whether a triangle is in the free list.
Definition: delaunay_2d.h:511
static index_t find_3(const signed_index_t *T, signed_index_t v)
Finds the index of an integer in an array of three integers.
Definition: delaunay_2d.h:873
index_t nearest_vertex(const double *p) const override
Computes the nearest vertex from a query point.
void check_combinatorics(bool verbose=false) const
For debugging purposes, tests some combinatorial properties.
index_t new_triangle()
Creates a new triangle.
Definition: delaunay_2d.h:522
void show_triangle_adjacent(index_t t, index_t le) const
For debugging purposes, displays a triangle adjacency.
bool triangle_is_conflict(index_t t, const double *p) const
Tests whether a given triangle is in conflict with a given 3d point.
Definition: delaunay_2d.h:791
Abstract interface for Delaunay triangulation in Nd.
Definition: delaunay.h:71
Regular Delaunay triangulation of weighted points.
Definition: delaunay_2d.h:992
~RegularWeightedDelaunay2d() override
RegularWeightedDelaunay2d destructor.
RegularWeightedDelaunay2d(coord_index_t dimension=3)
Constructs a new Regular Delaunay2d triangulation.
Abstract interface for Delaunay.
Common include file, providing basic definitions. Should be included before anything else by all head...
Geometric functions in 2d and 3d.
Sign in_circle_2d_SOS(const double *p0, const double *p1, const double *p2, const double *p3)
Tests whether a 2d point is inside the circumscribed circle of a 3d triangle.
Sign orient_2dlifted_SOS(const double *p0, const double *p1, const double *p2, const double *p3, double h0, double h1, double h2, double h3)
Computes the 3d orientation test with lifted points.
Sign orient_2d(const vec2HE &p0, const vec2HE &p1, const vec2HE &p2)
Computes the orientation predicate in 2d.
Global Vorpaline namespace.
Definition: basic.h:55
geo_signed_index_t signed_index_t
The type for storing and manipulating indices differences.
Definition: numeric.h:343
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
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition: numeric.h:363
Filtered exact predicates for restricted Voronoi diagrams.