Fréchet View  1.6.0
A Tool for Exploring Fréchet Distance Algorithms
triangulation.cpp
Go to the documentation of this file.
1 
2 #include <triangulation.h>
5 
6 using namespace frechet;
7 using namespace poly;
8 
10  : tds(), vertices(), neighbors(), poly_n(-1)
11 {
12  if (capacity > 0)
13  vertices.reserve(capacity);
14  tds.set_dimension(2);
15  south_pole = tds.create_vertex(Vertex_base<TDS>(Point(NAN,NAN),SOUTH_POLE_INDEX));
16 }
17 
19 : tds(that.tds), vertices(), neighbors(), poly_n(that.poly_n)
20 {
21  // deep copy vertices
23  vertices.resize(that.vertices.size());
24  auto vi = tds.vertices_begin();
25  auto vend = tds.vertices_end();
26  for( ; vi != vend; ++vi) {
27  vi->prev=nullptr;
28  if (vi->pindex==SOUTH_POLE_INDEX)
29  south_pole = vi;
30  else if (vi->pindex >= 0)
31  vertices[vi->pindex] = vi;
32  }
33  Q_ASSERT((south_pole!=Vertex_handle()) && south_pole->pindex==SOUTH_POLE_INDEX);
34 
35  // DON'T deep copy neigbors (and edges)
36 
37  Q_ASSERT(assertEquals(that));
38 }
39 
41 {
42  vertices.clear();
43  neighbors.clear();
44  tds.clear();
45 }
46 
47 
49 {
50  return v->pindex==SOUTH_POLE_INDEX;
51 }
52 
54 {
55  return is_south_pole(f->vertex(0));
56 }
57 
59 {
60  return ! is_outer_face(f);
61 }
62 
64 {
65  return is_outer_face(e.first) && (e.second!=0);
66 }
67 
69 {
70  return is_outer_face(e.first)
71  != is_outer_face(e.first->neighbor(e.second));
72 }
73 
75 {
76  return is_inner_face(e.first)
77  && is_inner_face(e.first->neighbor(e.second));
78 }
79 
80 void Triangulation::mirror(Edge& e) const {
81  e = tds.mirror_edge(e);
82 }
83 
85  mirror(e2.first);
86  mirror(e2.second);
87 }
88 
90 {
91  for(const frechet::reach::Placement& p : ps)
92  addTargetVertex(&p);
93 }
94 
96 {
97  Vertex_handle v;
98  int i = floor(w);
99  if (i == w) { // existing vertex
100  Q_ASSERT(i>=0 && i < vertices.size());
101  return vertices[i];
102  }
103  else {
104  return insertVertexOnEdge(w);
105  }
106 }
107 
109 {
110  Vertex_handle v = addVertex(p->w);
111 
112  Q_ASSERT(p->fs.lower() >= 0);
113  Q_ASSERT(p->fs.upper() >= p->fs.lower());
114  v->target = p;
115  return v;
116 }
117 
119 {
120  this->poly_n = n;
121 }
122 
124 {
125  return vertices.size();
126 }
127 
129 {
130  Q_ASSERT(sz >= 0);
131  if (sz < size()) {
132  for (int i=vertices.size()-1; i >= sz; --i)
134  vertices.resize(sz);
135  }
136  else {
137  int i = size();
138  vertices.resize(sz);
139  for( ; i < sz; ++i)
140  vertices[i] = Vertex_handle();
141  }
142  //Q_ASSERT(tds.is_valid());
143 }
144 
145 
146 
148 {
149  Q_ASSERT(i >= 0);
150  Vertex_handle v = (*this)[i];
151  if (v == 0) {
152  if (i >= vertices.size())
153  resize(i+1);
154  vertices[i] = v = tds.create_vertex(Vertex_base<TDS>(p,i));
155  }
156  return v;
157 }
158 
189 {
190 
191  // (1) start at i = floor(x). We know that there is a polygon vertex at (i+1).
192  // find the adjacent south-pole face
193  int i=floor(x);
194  Vertex_handle v1 = (*this)[i];
195 
196  Edge e;
197  bool ok = tds.is_edge(south_pole,v1, e.first,e.second);
198  Q_ASSERT(ok);
199 
200  Face_handle f;
201  if (e.second==2) {
202  // wrong direction
203  f = e.first->neighbor(2);
204  }
205  else {
206  f = e.first;
207  }
208 
209  // (2) find the vertex with index < x < index+1
210  Vertex_handle v2 = f->vertex(1);
211  while((v2->offset != 0.0) && (v2->offset < x)) {
212  // advance
213  f = f->neighbor(2);
214  v2 = f->vertex(1);
215  }
216 
217  // (3) insert new vertex
218  Vertex_handle v = tds.insert_in_edge(f,0);
219  Point p = numeric::between(*v1,*v2,x-i);
220  addVertex(v, x, p);
221 
222  Q_ASSERT(tds.is_valid());
223  return v;
224 }
225 
227 {
228  Vertex_handle v1 = e.first->vertex((e.second+1)%3);
229  Vertex_handle v2 = e.first->vertex((e.second+2)%3);
230  return Segment(v1->pindex,v2->pindex);
231 }
232 
263 {
264 
265  Vertex_handle v1 = (*this)[index];
266  Face_circulator f (v1);
267  while(is_inner_face(f)) ++f;
268  Q_ASSERT(is_outer_face(f));
269  //std::cout << v1->index << "," << f << std::endl;
270 
271  //Face_handle f2 = f->neighbor(1);
272  Edge e (f,0);
274 
275  Vertex_handle v = tds.join_vertices(e);
276  // keep payload data from v2
277  v->setData(v2);
278  //v->set_face(f2); // right?
279 
280  vertices[v->pindex] = v;
281  vertices[index]=Vertex_handle();
282 
283  Q_ASSERT(v->is_valid());
284  Q_ASSERT(v->face()->is_valid());
285  Q_ASSERT(v->face()->has_vertex(v));
286  Q_ASSERT(tds.is_valid());
287 }
288 
289 
291 {
292  vertices.resize(size()+1);
293  vertices[size()-1] = v;
294  v->setData(p,size()-1,-1,x);
295 }
296 
298 {
299  Q_ASSERT(i >= 0);
300  if (i >= vertices.size())
301  return Vertex_handle();
302  else
303  return vertices[i];
304 }
305 
307 {
308  Q_ASSERT(i >= 0);
309  if (i >= vertices.size())
310  return Vertex_handle();
311  else
312  return vertices[i];
313 }
314 
316 {
317  Q_ASSERT(e.second!=SOUTH_POLE_INDEX);
318  Q_ASSERT(is_inner_face(e.first));
319 
320  if (!is_inner_face(e.first))
321  e = tds.mirror_edge(e);
322 
323  return EdgePair(Edge(e.first,e.first->ccw(e.second)),
324  Edge(e.first,e.first->cw(e.second)));
325 }
326 
335 {
336  Face_handle f = e.first;
337  Q_ASSERT(is_inner_face(f));
338  if (!is_inner_face(f))
339  f = f->neighbor(e.second);
340  Q_ASSERT(is_inner_face(f));
341 
342 #ifdef QT_DEBUG
343  // i1,i2 must be the incident vertices to e
345  if (vv.first->pindex==i2)
346  std::swap(vv.first,vv.second);
347  Q_ASSERT(vv.first->pindex==i1);
348  Q_ASSERT(vv.second->pindex==i2);
349 #endif
350  Vertex_handle v1 = vertices[i1];
351  Vertex_handle v2 = vertices[i2];
352  return EdgePair(Edge(f,f->index(v2)),
353  Edge(f,f->index(v1)));
354 }
355 
357 {
358  FacePair result;
359 
360  const Vertex_handle v1 = (*this)[s.first];
361  const Vertex_handle v2 = (*this)[s.second];
362 
363  /* TODO
364  * tds_is_edge() may need O(degree(v))
365  * later in the algorithm we want O(1)
366  *
367  * speed up acces by putting all edges into a map
368  */
369  int i;
370  if (tds.is_edge(v1,v2, result.first, i)) {
371  result.second = result.first->neighbor(i);
372  }
373  return result;
374 }
375 
377 {
378  return tds.incident_faces(v);
379 }
380 
382  e = tds.mirror_edge(e);
383  return e.first->vertex(e.second);
384 }
385 
386 int Triangulation::thirdVertex(int a, int b) const {
387  Edge e;
388  bool ok = tds.is_edge(vertices[a],vertices[b], e.first,e.second);
389  Q_ASSERT(ok);
390  Vertex_handle c = e.first->vertex(e.second);
391  if (is_south_pole(c)) {
392  e = tds.mirror_edge(e);
393  c = e.first->vertex(e.second);
394  }
395  Q_ASSERT(!is_south_pole(c));
396  return c->pindex;
397 }
398 
400 {
401  Vertex_handle vi = vertices[i];
402  Vertex_handle vj = vertices[j];
403  Edge e;
404  bool ok = tds.is_edge(vi,vj, e.first,e.second);
405  Q_ASSERT(ok);
406  if (e.second==SOUTH_POLE_INDEX)
407  e = tds.mirror_edge(e);
408  Q_ASSERT(e.second!=SOUTH_POLE_INDEX);
409  Q_ASSERT(is_inner_face(e.first));
410  return e;
411 }
412 
414 {
415  int i1 = e.first->ccw(e.second);
416  int i2 = e.first->cw(e.second);
417  return VertexPair(e.first->vertex(i1),
418  e.first->vertex(i2));
419 }
420 
422  auto v = tds.vertices_begin();
423  auto vend = tds.vertices_end();
424  for( ; v != vend; ++v)
425  v->prev = nullptr;
426 }
427 
429  auto v = tds.vertices_begin();
430  auto vend = tds.vertices_end();
431  for( ; v != vend; ++v)
432  v->target = nullptr;
433 }
434 
436 {
437  const Vertex_handle v1 = (*this)[s.first];
438  const Vertex_handle v2 = (*this)[s.second];
439 
440  TwoVertices result;
441  Edge e;
442  if (tds.is_edge(v1,v2, e.first,e.second))
443  {
444  result.first = e.first->vertex(e.second);
445  e = tds.mirror_edge(e);
446  result.second = e.first->vertex(e.second);
447  }
448 
449  return result;
450 }
451 
453 {
455 }
456 
457 void Triangulation::createFace(int a, int b, int c)
458 {
459  int count_before=countFaces();
460 
461  Vertex_handle va = (a==SOUTH_POLE_INDEX) ? south_pole : (*this)[a];
462  Vertex_handle vb = (*this)[b];
463  Vertex_handle vc = (*this)[c];
464 
465  TDS::Face_handle f = tds.create_face(va,vb,vc);
466 
467  //Q_ASSERT(tds.dimension()==2);
468  Q_ASSERT(f != TDS::Face_handle());
469 
470  Segment sa(a,b);
471  Segment sb(b,c);
472  Segment sc(c,a);
473 
474  linkNeighbors(sa, Edge(f,2));
475  linkNeighbors(sb, Edge(f,0));
476  linkNeighbors(sc, Edge(f,1));
477 
478  if (va->face()==Face_handle()) va->set_face(f);
479  if (vb->face()==Face_handle()) vb->set_face(f);
480  if (vc->face()==Face_handle()) vc->set_face(f);
481 
482  Q_ASSERT(countFaces()==count_before+1);
483 }
484 
485 
487 {
488  if (neighbors.empty()) return;
489 
490  // smallest vertex index
491  int v0 = neighbors.begin()->first.first;
492  // how often have we seen the first edge?
493  int first_edge=0;
494  while(! neighbors.empty()) {
495  auto n = neighbors.end();
496  --n;
497  Segment s = n->first;
498  Q_ASSERT(s.first>=0);
499  Q_ASSERT(s.second>=0);
500  int v1 = s.first;
501  int v2 = s.second;
502  // keep faces ALWAYS counter-clockwise
503  // special care needs to be taken when the first vertex is encountered
504  // sorting of the map ensures, that we encounter the [v0..]
505  // segment first
506  if (v1==v0) {
507  if (first_edge!=0) {
508  // wrap-around
509  std::swap(v1,v2);
510  }
511  ++first_edge;
512  Q_ASSERT(first_edge <= 2);
513  }
515  }
516  TDS::Vertex_iterator vit = tds.vertices_begin();
517  Q_ASSERT(vit==south_pole);
518  Q_ASSERT(south_pole->face() != Face_handle());
519  Q_ASSERT(tds.is_valid());
520 }
521 
523 {
524  // diagonals are visited twice.
525  // the first time, we record the face,
526  // the second time, we create the neighborhood relations
527  auto n = neighbors.find(s);
528  if (n==neighbors.end()) {
529  // First visit: store reference for later
530  neighbors.insert(std::make_pair(s,e));
531  }
532  else {
533  // Second visit: recover reference and link
534  Edge e2 = n->second;
535  tds.set_adjacency(e.first,e.second, e2.first,e2.second);
536  neighbors.erase(n);
537  }
538  // if all neighboring faces have been set, the 'neighbors' map will be empty again.
539  // @see isComplete()
540 }
541 
543 {
544  Q_ASSERT(v1!=Vertex_handle());
545  Q_ASSERT(v2!=Vertex_handle());
546 
547  Q_ASSERT(f->has_vertex(v1));
548  Q_ASSERT(f->has_vertex(v2));
549 
550  for(int i=0; i < 3; ++i) {
551  Vertex_handle v3 = f->vertex(i);
552  if (v3!=v1 && v3!=v2)
553  return v3;
554  }
555  Q_ASSERT(false);
556  return Vertex_handle();
557 }
558 
560 {
561  return third_vertex(f, (*this)[s.first], (*this)[s.second]);
562 }
563 
565 {
566  // inner faces: check embedding
567  return is_inner_face(f)
568  && PolygonShortestPaths::halfplaneTest<double>(f->vertex(0), f->vertex(1), f->vertex(2)) < 0
569  && PolygonShortestPaths::halfplaneTest<double>(f->vertex(1), f->vertex(2), f->vertex(0)) < 0
570  && PolygonShortestPaths::halfplaneTest<double>(f->vertex(2), f->vertex(0), f->vertex(1)) < 0;
571 }
572 /*
573 bool Triangulation::is_ccw_polygon_edge(Edge e)
574 {
575  // only applicable to polygon edges
576  Q_ASSERT(is_polygon_edge(e));
577 
578  int j = e.first->vertex((e.second+1)%3)->pindex;
579  int i = e.first->vertex((e.second+2)%3)->pindex;
580  return (i+1==j) || (j==0 || i==vertices.size());
581 }
582 */
583 template<class TDS>
585 
586 template<class TDS>
587 Vertex_base<TDS>::Vertex_base(const frechet::Point &p, int pi) : Vertex_base(p,pi,(double)pi) { }
588 
589 template<class TDS>
590 Vertex_base<TDS>::Vertex_base(const frechet::Point &p, int pi, double offs)
591  : Point(p),
592  CGAL::Triangulation_ds_vertex_base_2<TDS>(),
593  pindex(pi), dindex(-1), offset(offs),
594  prev(nullptr), target(nullptr),
595  L(),R(),B(), LR(),BR(),RR()
596 
597 { }
598 
599 template<class TDS>
601 
603  : TDS::Edge_iterator(i)
604 { }
605 
607  Face_handle f = e.first;
608  int i = e.second;
609  return f->vertex((i+1)%3);
610 }
611 
613  Face_handle f = e.first;
614  int i = e.second;
615  return f->vertex((i+2)%3);
616 }
617 
619  return Triangulation::vertex1(operator*());
620 }
621 
623  return Triangulation::vertex2(operator*());
624 }
625 
627 {
628  return Triangulation::is_outer_edge(operator*());
629 }
630 
632 {
633  return Triangulation::is_polygon_edge(operator*());
634 }
635 
637 {
638  return Triangulation::is_diagonal(operator*());
639 }
640 
642  return QLineF(*vertex1(),*vertex2());
643 }
644 
646  return Triangulation::segment(operator*());
647 }
648 
650 {
651  Q_ASSERT(tds.number_of_vertices()==that.tds.number_of_vertices());
652  Q_ASSERT(tds.number_of_faces()==that.tds.number_of_faces());
653  Q_ASSERT(tds.number_of_edges()==that.tds.number_of_edges());
654 
655  Q_ASSERT(vertices.size()==that.vertices.size());
656  for(int i=0; i < vertices.size(); ++i)
657  {
658  Vertex_handle v1 = vertices[i];
659  Vertex_handle v2 = that.vertices[i];
660 
661  Q_ASSERT((v1!=Vertex_handle()) == (v2!=Vertex_handle()));
662  Q_ASSERT((v1==Vertex_handle()) || v1->assertEquals(*v2));
663  }
664 
665  auto v1i = tds.vertices_begin();
666  auto v1end = tds.vertices_end();
667 /*
668  auto v2i = that.tds.vertices_begin();
669  auto v2end = that.tds.vertices_end();
670 
671  for( ; v1i!=v1end && v2i!=v2end; ++v1i,++v2i)
672  Q_ASSERT(v1i->assertEquals(*v2i));
673 */
674  for( ; v1i != v1end; ++v1i)
675  if (v1i->pindex >= 0)
676  Q_ASSERT(v1i == vertices[v1i->pindex]);
677 
678  Q_ASSERT(v1i==v1end);
679 // Q_ASSERT(v2i==v2end);
680  return true;
681 }
682 
683 template<typename TDS>
685 {
686  Q_ASSERT(this->pindex == that.pindex);
687 
688  if (this->pindex==SOUTH_POLE_INDEX) return true;
689 
690  Q_ASSERT((Point)*this == (Point)that);
691  Q_ASSERT(this->dindex == that.dindex);
692  Q_ASSERT(this->offset == that.offset);
693  return true;
694 }
695 
696 
697 std::ostream& frechet::poly::operator<< (std::ostream& out, const Triangulation& tri)
698 {
699  out << "{";
700  auto tbegin = tri.edges_begin();
701  auto tend = tri.edges_end();
702  for(auto t=tbegin; t!=tend; ++t) {
703  if (t!=tbegin) out << ",";
704  out << t;
705  }
706  return out << "}";
707 }
708 
709 std::ostream& frechet::poly::operator<< (std::ostream& out, const Triangulation::Face_handle& f)
710 {
711  return out << "["
712  <<f->vertex(0)->pindex<<"-"
713  <<f->vertex(1)->pindex<<"-"
714  <<f->vertex(2)->pindex<<"]";
715 }
716 
717 std::ostream& frechet::poly::operator<< (std::ostream& out, const Triangulation::Edge& e)
718 {
720  out << "P";
721  else if (Triangulation::is_outer_edge(e))
722  out << "O";
723  else if (Triangulation::is_diagonal(e))
724  out << "D";
725 
726  return out << Triangulation::vertex1(e)->pindex << "--" << Triangulation::vertex2(e)->pindex;
727 }
728 
729 std::ostream& frechet::poly::operator<< (std::ostream& out, const Triangulation::Edge_iterator& e)
730 {
731  return frechet::poly::operator<<(out,*e);
732 }
733 
Represents a polygon line segment from node i to j.
Definition: types.h:39
void complement()
after all faces have been added, create the complementary "outer" edges (connected to the south-pole)...
void swap(gpuword **A, gpuword **B)
NeighborMap neighbors
Keeps track of diagonals and their neighborhood during construction.
void setPolygonSize(int n)
set number of polygon vertexes
static Vertex_handle vertex1(Edge e)
Edge edge(int i, int j) const
construct an edge from two vertexes
TwoVertices adjacentVertices(Segment s) const
std::vector< frechet::reach::Placement > Placements
a list of Placement objects
Definition: graph_model.h:144
static Vertex_handle third_vertex(Face_handle f, Vertex_handle v1, Vertex_handle v2)
given two vertexes and a face, find the third vertex of the triangle
std::pair< Vertex_handle, Vertex_handle > VertexPair
a pair of vertexes
Vertex_handle addTargetVertex(const frechet::reach::Placement *p)
add a vertex for computing valid placements
Edge_iterator(TDS::Edge_iterator i)
default constructor
global definitions for all algorithms.
void mirror(Edge &e) const
flip the direction of an edge. Keep in mind that Edge object are desribed as a Face + opposite vertex...
static Segment segment(Edge e)
TDS::Edge Edge
an edge of the triangulation either an edge of the polygon, or a diagonal, or an articial "outer" edg...
Face_circulator incidentFaces(Vertex_handle v) const
a circular (never-ending) iterator that traverses all faces incident to a given vertex.
bool assertEquals(const Vertex_base &that)
debug assertion
static bool is_inner_face(Face_handle f)
EdgePair adjacentEdges(Edge e) const
TDS::Face_handle Face_handle
handle to a triangulation face (handles are managed by TDS base class)
static bool is_polygon_edge(Edge e)
int thirdVertex(int a, int b) const
given two vertexes, find the third vertex to complement a triangle
int pindex
Index of polygon vertex (if this object refers to a polygon vertex; invalid for extra points)
void clearTargets()
clear targets for shortest-path tree search
Vertex_handle createVertex(const Point &p, int i)
create a new triangulation vertex
void resize(int sz)
reset number of vertexes in triangulation
static Vertex_handle vertex2(Edge e)
Vertex_handle addVertex(double w)
add a triangulation vertex on the edge of the polygon (does not correspond to a polygon vertex)
double upper() const
Definition: interval.h:96
static bool is_diagonal(Edge e)
FacePair adjacentFaces(Segment s) const
void addTargetVertices(const frechet::reach::Placements &p)
add vertexes for computing valid placements
QPointF Point
a point in the plane; with double floating point precision. This type is heavily used throughout all ...
Definition: types.h:14
void createFace(int a, int b, int c)
insert a new face. Caller is responsible for keeping the triangulation in a consistent state....
static bool is_south_pole(Vertex_handle v)
static const int SOUTH_POLE_INDEX
artifical vertex used by Triangulation
Definition: types.h:21
base class for vertexes in a Triangulation Data Structure.
Definition: triangulation.h:40
Point between(Point p, Point q, double w)
compute a point on a line segment
Definition: numeric.h:295
bool is_ccw_face(Face_handle f)
std::pair< Face_handle, Face_handle > FacePair
a pair of faces
std::pair< Edge, Edge > EdgePair
a pair of edges
void linkNeighbors(Segment s, Edge e)
keep track of neighboring faces during construction of the triangulation.
Edge_iterator edges_end() const
double lower() const
Definition: interval.h:92
Edge_iterator edges_begin() const
bool assertEquals(const Triangulation &that) const
equality assertion for debugging
Vertex_handle insertVertexOnEdge(double x)
split an edge into two
~Vertex_base()
destructor (empty, actually)
std::pair< Vertex_handle, Vertex_handle > TwoVertices
a pair of vertex handles
placement of a diagonal point in free-space diagram When calculating valid placements in [buchin06]
Definition: graph_model.h:132
Vertex_base()
empty constructore; creates an invalid vertex
Vertex_handle oppositeVertex(Edge e) const
the second vertex opposite to an edge. Keep in mind that Edge object are desribed as a Face + opposit...
TDS tds
CGAL Triangulation Data Structure.
CGAL::Triangulation_data_structure_2< Vertex_base<>, CGAL::Triangulation_ds_face_base_2<> > TDS
Triangulation Data Structure. Instantiates CGAL::Triangulation_data_structure_2 with our customized V...
VertexPair incidentVertexes(Edge e) const
TDS::Face_circulator Face_circulator
a circulating iterator
TDS::Vertex_handle south_pole
a special vertex that completes the TDS structure
Interval fs
Free-Space Interval.
Definition: graph_model.h:135
Wrapper for CGAL::Triangulation_Data_Structure https://doc.cgal.org/latest/TDS_2/index....
TDS::Vertex_handle Vertex_handle
handle to a vertex (handles are managed by TDS base class)
static bool is_outer_face(Face_handle f)
static bool is_outer_edge(Edge e)
void clearPredecessors()
used during shortest-paths computation: reset shortest-path tree information on all vertexes
Triangulation(int capacity=-1)
default constructor
Vertex_handle operator[](int i)
edge iterator. Traverses all edges of the triangulation. Derives from TDS::Edge_iterator and supplies...
std::ostream & operator<<(std::ostream &out, const frechet::poly::CgalPolygon &poly)
print operator (for debugging)
Definition: poly_utils.cpp:338
int poly_n
polygon size, number of vertexes in polyong
void removeVertexFromEdge(int index)
remove a vertex