Geogram Version 1.9.2
A programming library of geometric algorithms
Loading...
Searching...
No Matches
generic_RVD_vertex.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_VORONOI_GENERIC_RVD_VERTEX
41#define GEOGRAM_VORONOI_GENERIC_RVD_VERTEX
42
44#include <geogram/mesh/mesh.h>
49
60namespace GEOGen {
61
62 using GEO::Delaunay;
63 using GEO::index_t;
65 using GEO::coord_index_t;
66 using GEO::Sign;
68 using GEO::Mesh;
69
78 template <class T, index_t DIM>
79 class small_set {
80
83
84 public:
86 typedef T* iterator;
87
89 typedef const T* const_iterator;
90
92 typedef T& reference;
93
95 typedef T value_type;
96
101 size_(0) {
102 }
103
107 index_t size() const {
108 return size_;
109 }
110
115 index_t capacity() const {
116 return (index_t) DIM;
117 }
118
123 return data_;
124 }
125
130 return data_ + size_;
131 }
132
138 return data_ + DIM;
139 }
140
145 return data_;
146 }
147
152 return data_ + size_;
153 }
154
160 return data_ + DIM;
161 }
162
169 iterator insert(const T& x) {
170 return insert(x, find_i(x));
171 }
172
181 iterator insert(const T& x, iterator where) {
182 if(where == end()) {
183 *where = x;
184 grow();
185 return where;
186 }
187 if(*where == x) {
188 return where;
189 }
190 grow();
191 if(where == end() - 1) {
192 *where = x;
193 return where;
194 }
195 for(iterator i = end() - 1; i != where; i--) {
196 geo_debug_assert(i != begin());
197 *i = *(i - 1);
198 }
199 *where = x;
200#ifdef GEO_DEBUG
201 for(iterator i = begin(); i != end() - 1; ++i) {
202 geo_debug_assert(*i < *(i + 1));
203 }
204#endif
205 return where;
206 }
207
211 void clear() {
212 size_ = 0;
213 }
214
220 iterator find(const T& x) {
221 iterator result = find_i(x);
222 if(*result != x) {
223 result = end();
224 }
225 return result;
226 }
227
233 const_iterator find(const T& x) const {
234 const_iterator result = find_i(x);
235 if(*result != x) {
236 result = end();
237 }
238 return result;
239 }
240
246 void push_back(const T& x) {
247#ifdef GEO_DEBUG
248 for(iterator i = begin(); i != end(); ++i) {
249 geo_debug_assert(*i < x);
250 }
251#endif
252 *end() = x;
253 grow();
254 }
255
259 void print(std::ostream& out) const {
260 out << "[ ";
261 for(const_iterator it = begin(); it != end(); ++it) {
262 out << *it << " ";
263 }
264 out << "]";
265 }
266
272 T& operator[] (signed_index_t i) {
273 geo_debug_assert(i >= 0);
274 geo_debug_assert(begin() + i < end());
275 return begin()[i];
276 }
277
283 const T& operator[] (signed_index_t i) const {
284 geo_debug_assert(i >= 0);
285 geo_debug_assert(begin() + i < end());
286 return begin()[i];
287 }
288
289 protected:
294 void grow() {
296 size_++;
297 }
298
299 // Note: maybe we should start from end() instead of begin()
300 // since negative indices are inserted first.
301
309 iterator find_i(const T& x) {
310 iterator result = begin();
311 while(result != end() && *result < x) {
312 result++;
313 }
314 return result;
315 }
316
323 const_iterator find_i(const T& x) const {
324 const_iterator result = begin();
325 while(result != end() && *result < x) {
326 result++;
327 }
328 return result;
329 }
330
331 protected:
332 T data_[DIM];
333 index_t size_;
334 };
335
339 template <class T, index_t DIM>
340 inline std::ostream& operator<< (
341 std::ostream& out,
342 const small_set<T, DIM>& S) {
343 S.print(out);
344 return out;
345 }
346
353 template <class T, index_t DIM1, index_t DIM2, index_t DIM3>
354 inline void sets_intersect(
355 const small_set<T, DIM1>& S1,
356 const small_set<T, DIM2>& S2,
358 ) {
359 I.clear();
360 auto i1 = S1.begin();
361 auto i2 = S2.begin();
362 while(i1 < S1.end() && i2 < S2.end()) {
363 if(*i1 < *i2) {
364 ++i1;
365 }
366 else if(*i2 < *i1) {
367 ++i2;
368 }
369 else {
370 I.push_back(*i1);
371 ++i1;
372 ++i2;
373 }
374 }
375 }
376
408 class SymbolicVertex : public small_set<GEO::signed_index_t, 3> {
409
412
415
416 public:
421 v1_(0),
422 v2_(0) {
423 }
424
428 void add_bisector(index_t i) {
429 baseclass::insert(signed_index_t(i) + 1);
430 }
431
435 void add_boundary_facet(index_t i) {
436 baseclass::insert(-signed_index_t(i) - 1);
437 }
438
443 index_t nb_boundary_facets() const {
444 index_t result = 0;
445 for(auto it = baseclass::begin();
446 it != baseclass::end() && *it < 0; ++it) {
447 result++;
448 }
449 return result;
450 }
451
455 index_t nb_bisectors() const {
456 index_t result = 0;
457 for(auto it = baseclass::end() - 1;
458 it != baseclass::begin() - 1 && *it > 0; --it) {
459 result++;
460 }
461 return result;
462 }
463
469 static index_t to_unsigned_int(signed_index_t x) {
470 geo_debug_assert(x >= 0);
471 return (index_t) (x);
472 }
473
481 index_t bisector(signed_index_t i) const {
482 geo_debug_assert(i < (signed_index_t) nb_bisectors());
483 return to_unsigned_int((baseclass::end()[-1 - i]) - 1);
484 }
485
492 index_t boundary_facet(signed_index_t i) const {
493 geo_debug_assert(i < (signed_index_t) nb_boundary_facets());
494 return to_unsigned_int(-(baseclass::begin()[i]) - 1);
495 }
496
502 bool has_bisector(index_t i) const {
503 return baseclass::find(signed_index_t(i) + 1) != baseclass::end();
504 }
505
511 bool has_boundary_facet(index_t i) const {
512 return baseclass::find(-signed_index_t(i) - 1) != baseclass::end();
513 }
514
520 index_t get_boundary_vertex() const {
522 geo_debug_assert(v1_ != 0);
523 return v1_ - 1;
524 }
525
533 void get_boundary_edge(index_t& v1, index_t& v2) const {
535 geo_debug_assert(v1_ != 0);
536 geo_debug_assert(v2_ != 0);
537 v1 = v1_ - 1;
538 v2 = v2_ - 1;
539 }
540
545 void set_boundary_vertex(index_t v) {
546 v1_ = v + 1;
547 v2_ = 0;
548 }
549
555 void set_boundary_edge(index_t v1, index_t v2) {
556 v1_ = v1 + 1;
557 v2_ = v2 + 1;
558 }
559
566 geo_debug_assert(rhs.nb_bisectors() == 1);
567 geo_debug_assert(rhs.v1_ > 0);
568 geo_debug_assert(rhs.v2_ > 0);
569 v1_ = rhs.v1_;
570 v2_ = rhs.v2_;
571 }
572
583 const thisclass& v1,
584 const thisclass& v2,
585 index_t E
586 ) {
587
588 // Compute the symbolic representation as the intersection
589 // of three planes.
590 sets_intersect(v1, v2, *this);
591 // this computes the set of planes that contain
592 // the edge [v1,v2]
593
594 add_bisector(E); // the intersection is on E.
595
596 // Compute the symbolic representation as intersection between
597 // bisector and boundary edge
598 // (it's redundant and less elegant than the representation
599 // as planes interactions,
600 // but we need this to handle degenerate configurations properly,
601 // and to use exact predicates with original boundary vertices
602 // coordinates).
603
604 if(nb_boundary_facets() == 2) {
605 // If *this is on the intersection of two boundary facets,
606 // then *this is on
607 // a boundary edge, and we need to retrieve the indices of the
608 // two extremities of this boundary edge.
609
610 index_t nb1 = v1.nb_boundary_facets();
611 index_t nb2 = v2.nb_boundary_facets();
612 if(nb1 == 3 && nb2 == 3) {
613 // If v1 and v2 are boundary vertices,
614 // then I is on the boundary
615 // edge that connects v1 and v2
619 );
620 } else if(nb1 == 2) {
622 // If v1 is on a boundary edge,
623 // then I is on the same boundary edge as v1
625 } else if(nb2 == 2) {
627 // If v2 is on a boundary edge,
628 // then I is on the same boundary edge as v2
630 }
631 }
632
633 // Sanity check: problem detected here, we
634 // notify the caller that will use a workaround
635 // (see clip_by_plane())
636 if(baseclass::size() != 3) {
637 return false;
638 }
639 return true;
640 }
641
642 private:
643 index_t v1_;
644 index_t v2_;
645 };
646
647
669 public:
674 PointAllocator(coord_index_t dim) :
675 size_(0),
676 capacity_(0),
677 dimension_(dim) {
678 }
679
685 double* new_item() {
686 if(size_ == capacity_) {
687 grow();
688 }
689 size_++;
690 return item(size_ - 1);
691 }
692
696 void clear() {
697 size_ = 0;
698 }
699
705 for(index_t c = 0; c < chunks_.size(); c++) {
706 GEO::Memory::aligned_free(chunks_[c]);
707 }
708 }
709
714 coord_index_t dimension() const {
715 return dimension_;
716 }
717
718 protected:
722 enum {
723 CHUNK_SHIFT = 8,
724 CHUNK_SIZE = 1 << CHUNK_SHIFT,
725 CHUNK_MASK = CHUNK_SIZE - 1
726 };
727
731 void grow() {
732 chunks_.push_back(
733 reinterpret_cast<double*>(
735 index_t(CHUNK_SIZE) * dimension_ * sizeof(double)
736 )
737 )
738 );
739 capacity_ += CHUNK_SIZE;
740 }
741
747 double* item(index_t i) {
748 geo_debug_assert(i < size_);
749 return &(chunks_[i >> CHUNK_SHIFT][(i & CHUNK_MASK) * dimension_]);
750 }
751
752 private:
753 index_t size_;
754 index_t capacity_;
755 coord_index_t dimension_;
756 std::vector<double*> chunks_;
757 };
758
762 enum {
764 INTERSECT = 2
765 };
766
771 typedef index_t EdgeFlags;
772
776 typedef index_t EdgeFlag;
777
787 class Vertex {
788
790 typedef Vertex thisclass;
791
792 public:
801 const double* p, double w, signed_index_t f,
802 const SymbolicVertex& sym
803 ) :
804 point_(p),
805 weight_(w),
806 f_(f),
807 seed_(-1),
808 sym_(sym),
809 flags_(ORIGINAL) {
810 }
811
818 Vertex(const double* p, double w, signed_index_t f) :
819 point_(p),
820 weight_(w),
821 f_(f),
822 seed_(-1),
823 flags_(ORIGINAL) {
824 }
825
830 point_(nullptr),
831 weight_(1.0),
832 f_(-1),
833 seed_(-1),
834 flags_(0) {
835 }
836
841 const double* point() const {
842 return point_;
843 }
844
849 void set_point(const double* p) {
850 point_ = p;
851 }
852
858 double weight() const {
859 return weight_;
860 }
861
867 void set_weight(double w) {
868 weight_ = w;
869 }
870
875 signed_index_t adjacent_seed() const {
876 return seed_;
877 }
878
883 void set_adjacent_seed(signed_index_t s) {
884 seed_ = s;
885 }
886
892 const SymbolicVertex& sym() const {
893 return sym_;
894 }
895
900 return sym_;
901 }
902
907 signed_index_t adjacent_facet() const {
908 return f_;
909 }
910
915 void set_adjacent_facet(signed_index_t f) {
916 f_ = f;
917 }
918
926 operator const double* () const {
927 return point_;
928 }
929
933 void clear() {
934 flags_ = 0;
935 f_ = -1;
936 }
937
942 flags_ |= f;
943 }
944
949 flags_ &= ~f;
950 }
951
955 bool check_flag(EdgeFlag f) const {
956 return (flags_ & f) != 0;
957 }
958
962 void copy_edge_from(const Vertex& rhs) {
964 flags_ = rhs.flags_;
965 }
966
976 template <index_t DIM>
978 PointAllocator& target_intersections,
979 const Vertex& vq1, const Vertex& vq2,
980 const double* p1, const double* p2
981 ) {
982 const double* q1 = vq1.point();
983 const double* q2 = vq2.point();
984 double* Ipoint = target_intersections.new_item();
985 set_point(Ipoint);
986 double d = 0.0, l1 = 0.0, l2 = 0.0;
987 for(coord_index_t c = 0; c < DIM; ++c) {
988 double n = p1[c] - p2[c];
989 d -= n * (p2[c] + p1[c]);
990 l1 += q2[c] * n;
991 l2 += q1[c] * n;
992 }
993 d = 0.5 * d;
994 l1 = ::fabs(l1 + d);
995 l2 = ::fabs(l2 + d);
996 double l12 = l1 + l2;
997 if(l12 > 1e-30) {
998 l1 /= l12;
999 l2 /= l12;
1000 } else {
1001 l1 = 0.5;
1002 l2 = 0.5;
1003 }
1004 for(coord_index_t c = 0; c < DIM; ++c) {
1005 Ipoint[c] = l1 * q1[c] + l2 * q2[c];
1006 }
1007 set_weight(l1 * vq1.weight() + l2 * vq2.weight());
1008 }
1009
1020 template <index_t DIM>
1022 const double* p1, const double* p2
1023 ) const {
1024 double r = 0.0;
1025 for(index_t c = 0; c < DIM; ++c) {
1026 r += GEO::geo_sqr(p2[c] - point()[c]);
1027 r -= GEO::geo_sqr(p1[c] - point()[c]);
1028 }
1029 return GEO::geo_sgn(r);
1030 }
1031
1032 private:
1033 const double* point_;
1034 double weight_;
1035
1040 signed_index_t f_;
1041
1047 signed_index_t seed_;
1048
1050 SymbolicVertex sym_;
1051
1056 EdgeFlags flags_;
1057 };
1058}
1059
1060/*
1061 namespace GEO {
1062 template <> struct can_be_used_as_attribute<GEOGen::SymbolicVertex> {
1063 static constexpr auto value = std::integral_constant<bool,true>();
1064 };
1065 }
1066*/
1067
1068
1069#endif
Assertion checking mechanism.
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition assert.h:196
Generic mechanism for attributes.
Common include file, providing basic definitions. Should be included before anything else by all head...
An allocator for points that are created from intersections in GenericVoronoiDiagram.
void clear()
Clears this PointAllocator.
~PointAllocator()
PointAllocator destructor.
double * new_item()
Allocates a new point.
void grow()
Allocates a new chunk of memory.
coord_index_t dimension() const
Gets the dimension of the points stored in this PointAllocator.
double * item(index_t i)
Gets a pointer to one of the allocated points from its index.
PointAllocator(coord_index_t dim)
Creates a new empty PointAllocator.
A set of three integers that encodes the equation of a vertex in GenericVoronoiDiagram.
bool has_bisector(index_t i) const
Tests whether a bisector is present in the symbolic representation this vertex.
void copy_boundary_edge_from(const thisclass &rhs)
Copies a boundary edge from the symbolic representation of another vertex.
void set_boundary_edge(index_t v1, index_t v2)
Sets the boundary edge on which this vertex is located.
index_t boundary_facet(signed_index_t i) const
Gets a boundary facet.
void add_boundary_facet(index_t i)
Adds a boundary facet to the symbolic representation.
void get_boundary_edge(index_t &v1, index_t &v2) const
Gets the global indices of the boundary vertices that define the boundary edge on which this vertex i...
index_t nb_boundary_facets() const
Gets the number of boundary facets in the symbolic representation.
index_t get_boundary_vertex() const
Gets the global index of the boundary vertex that corresponds to this vertex.
index_t nb_bisectors() const
Gets the number of bisectors in the symbolic representation.
void set_boundary_vertex(index_t v)
Sets the boundary vertex on which this vertex is located.
bool has_boundary_facet(index_t i) const
Tests whether a boundary facet is present in the symbolic representation of this vertex.
SymbolicVertex()
Creates an uninitialized SymbolicVertex.
static index_t to_unsigned_int(signed_index_t x)
Casts a signed_index_t into an (unsigned) index_t.
index_t bisector(signed_index_t i) const
Gets a bisector.
void add_bisector(index_t i)
Adds a bisector to the symbolic representation.
bool intersect_symbolic(const thisclass &v1, const thisclass &v2, index_t E)
Computes the symbolic representation of the intersection between a segment and a bisector.
Internal representation of vertices in GenericVoronoiDiagram.
void clear()
Clears this Vertex.
Vertex(const double *p, double w, signed_index_t f, const SymbolicVertex &sym)
Creates a new Vertex.
void set_adjacent_seed(signed_index_t s)
Sets the adjacent seed.
void set_flag(EdgeFlag f)
Sets an EdgeFlag in this Vertex.
Vertex()
Creates an uninitialized Vertex.
void intersect_geom(PointAllocator &target_intersections, const Vertex &vq1, const Vertex &vq2, const double *p1, const double *p2)
Computes the intersection between a segment and a bisector.
Sign side_fast(const double *p1, const double *p2) const
Computes the side of this vertex relative to a bisector.
SymbolicVertex & sym()
Gets the symbolic representation.
void set_adjacent_facet(signed_index_t f)
Sets the adjacent facet.
const double * point() const
Gets the geometric location at this Vertex.
void set_point(const double *p)
Sets the geometric location at this vertex.
double weight() const
Gets Vertex weight.
void copy_edge_from(const Vertex &rhs)
Copies adjacent facet and edge flags from another Vertex.
void unset_flag(EdgeFlag f)
Resets an EdgeFlag in this Vertex.
const SymbolicVertex & sym() const
Gets the symbolic representation.
void set_weight(double w)
Sets the vertex weight.
signed_index_t adjacent_seed() const
Gets the adjacent seed.
Vertex(const double *p, double w, signed_index_t f)
Creates a new Vertex.
signed_index_t adjacent_facet() const
Gets the adjacent facet.
bool check_flag(EdgeFlag f) const
Tests an EdgeFlag in this Vertex.
Small_set is similar to std::set, but with fixed maximum size (and no dynamic memory allocation).
const_iterator begin() const
Gets a const iterator to the first element..
iterator find(const T &x)
Finds an element by value.
void grow()
Increases the size of this small_set.
iterator end()
Gets an iterator one position past the last element.
T value_type
Type of the elements.
index_t capacity() const
Gets the maximum number of elements that can be stored in this small_set.
const T * const_iterator
A random access iterator to const elements.
small_set()
Constructs an empty small_set.
iterator insert(const T &x, iterator where)
Inserts a new element at a specified location..
const_iterator find(const T &x) const
Finds an element by value.
const_iterator find_i(const T &x) const
Finds where an element should be located from its value.
T * iterator
A random access iterator to elements
const_iterator end_of_storage() const
Gets a const iterator one position past the last element that can be stored.
iterator begin()
Gets an iterator to the first element.
void clear()
Clears this small_set.
T & reference
Reference to element.
iterator find_i(const T &x)
Finds where an element is or where it should be inserted from its value.
void print(std::ostream &out) const
Displays the stored elements.
T & operator[](signed_index_t i)
Direct access to an element.
iterator insert(const T &x)
Insert a new element.
void push_back(const T &x)
Appends an element to the end of the list.
index_t size() const
Gets the number of element in this small_set.
const_iterator end() const
Gets a const iterator one position past the last element.
iterator end_of_storage()
Gets an iterator one position past the last element that can be stored.
Abstract interface for Delaunay triangulation in Nd.
Definition delaunay.h:71
Represents a mesh.
Definition mesh.h:2701
Implementation of Delaunay using nearest neighbors.
index_t EdgeFlags
A set of EdgeFlags.
void sets_intersect(const small_set< T, DIM1 > &S1, const small_set< T, DIM2 > &S2, small_set< T, DIM3 > &I)
Computes the intersection between two small_sets.
index_t EdgeFlag
An individual edge flag.
The class that represents a mesh.
void * aligned_malloc(size_t size, size_t alignment=GEO_MEMORY_ALIGNMENT)
Allocates aligned memory.
Definition memory.h:281
void aligned_free(void *p)
Deallocates aligned memory.
Definition memory.h:309
T geo_sqr(T x)
Gets the square value of a value.
Definition numeric.h:301
geo_signed_index_t signed_index_t
The type for storing and manipulating indices differences.
Definition numeric.h:343
Sign geo_sgn(const T &x)
Gets the sign of a value.
Definition numeric.h:90
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
Function and classes for process manipulation.