1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2012 Excensus LLC.
7 * Copyright (C) 2019 Daniel Baston
8 *
9 * This is free software; you can redistribute and/or modify it under
10 * the terms of the GNU Lesser General Licence as published
11 * by the Free Software Foundation.
12 * See the COPYING file for more information.
13 *
14 **********************************************************************
15 *
16 * Last port: triangulate/quadedge/QuadEdgeSubdivision.java r524
17 *
18 **********************************************************************/
19 #include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
20
21 #include <algorithm>
22 #include <vector>
23 #include <set>
24 #include <iostream>
25
26 #include <geos/geom/Polygon.h>
27 #include <geos/geom/LineSegment.h>
28 #include <geos/geom/LineString.h>
29 #include <geos/geom/CoordinateSequence.h>
30 #include <geos/geom/CoordinateArraySequence.h>
31 #include <geos/geom/CoordinateSequenceFactory.h>
32 #include <geos/geom/CoordinateArraySequenceFactory.h>
33 #include <geos/geom/CoordinateList.h>
34 #include <geos/geom/GeometryCollection.h>
35 #include <geos/geom/GeometryFactory.h>
36 #include <geos/util.h>
37 #include <geos/triangulate/quadedge/QuadEdge.h>
38 #include <geos/triangulate/quadedge/QuadEdgeLocator.h>
39 #include <geos/triangulate/quadedge/LastFoundQuadEdgeLocator.h>
40 #include <geos/triangulate/quadedge/LocateFailureException.h>
41 #include <geos/triangulate/quadedge/TriangleVisitor.h>
42 #include <geos/geom/Triangle.h>
43
44
45 using namespace geos::geom;
46 using namespace std;
47
48 namespace geos {
49 namespace triangulate { //geos.triangulate
50 namespace quadedge { //geos.triangulate.quadedge
51
52 void
getTriangleEdges(const QuadEdge & startQE,const QuadEdge * triEdge[3])53 QuadEdgeSubdivision::getTriangleEdges(const QuadEdge& startQE,
54 const QuadEdge* triEdge[3])
55 {
56 triEdge[0] = &startQE;
57 triEdge[1] = &triEdge[0]->lNext();
58 triEdge[2] = &triEdge[1]->lNext();
59 if(&triEdge[2]->lNext() != triEdge[0]) {
60 throw util::IllegalArgumentException("Edges do not form a triangle");
61 }
62 }
63
QuadEdgeSubdivision(const geom::Envelope & env,double p_tolerance)64 QuadEdgeSubdivision::QuadEdgeSubdivision(const geom::Envelope& env, double p_tolerance) :
65 tolerance(p_tolerance),
66 locator(new LastFoundQuadEdgeLocator(this)),
67 visit_state_clean(true)
68 {
69 edgeCoincidenceTolerance = tolerance / EDGE_COINCIDENCE_TOL_FACTOR;
70 createFrame(env);
71 initSubdiv();
72 }
73
74 void
createFrame(const geom::Envelope & env)75 QuadEdgeSubdivision::createFrame(const geom::Envelope& env)
76 {
77 double deltaX = env.getWidth();
78 double deltaY = env.getHeight();
79 double offset = 0.0;
80 if(deltaX > deltaY) {
81 offset = deltaX * 10.0;
82 }
83 else {
84 offset = deltaY * 10.0;
85 }
86
87 frameVertex[0] = Vertex((env.getMaxX() + env.getMinX()) / 2.0, env
88 .getMaxY() + offset);
89 frameVertex[1] = Vertex(env.getMinX() - offset, env.getMinY() - offset);
90 frameVertex[2] = Vertex(env.getMaxX() + offset, env.getMinY() - offset);
91
92 frameEnv = Envelope(frameVertex[0].getCoordinate(), frameVertex[1]
93 .getCoordinate());
94 frameEnv.expandToInclude(frameVertex[2].getCoordinate());
95 }
96 void
initSubdiv()97 QuadEdgeSubdivision::initSubdiv()
98 {
99 assert(quadEdges.empty());
100
101 // build initial subdivision from frame
102 startingEdges[0] = QuadEdge::makeEdge(frameVertex[0], frameVertex[1], quadEdges);
103 startingEdges[1] = QuadEdge::makeEdge(frameVertex[1], frameVertex[2], quadEdges);
104 QuadEdge::splice(startingEdges[0]->sym(), *startingEdges[1]);
105
106 startingEdges[2] = QuadEdge::makeEdge(frameVertex[2], frameVertex[0], quadEdges);
107 QuadEdge::splice(startingEdges[1]->sym(), *startingEdges[2]);
108 QuadEdge::splice(startingEdges[2]->sym(), *startingEdges[0]);
109 }
110
111 QuadEdge&
makeEdge(const Vertex & o,const Vertex & d)112 QuadEdgeSubdivision::makeEdge(const Vertex& o, const Vertex& d)
113 {
114 QuadEdge* e = QuadEdge::makeEdge(o, d, quadEdges);
115 return *e;
116 }
117
118 QuadEdge&
connect(QuadEdge & a,QuadEdge & b)119 QuadEdgeSubdivision::connect(QuadEdge& a, QuadEdge& b)
120 {
121 QuadEdge* e = QuadEdge::connect(a, b, quadEdges);
122 return *e;
123 }
124
125 void
remove(QuadEdge & e)126 QuadEdgeSubdivision::remove(QuadEdge& e)
127 {
128 QuadEdge::splice(e, e.oPrev());
129 QuadEdge::splice(e.sym(), e.sym().oPrev());
130
131 // because QuadEdge pointers must be stable, do not remove edge from quadedges container
132 // This is fine since they are detached from the subdivision
133
134 e.remove();
135 }
136
137 QuadEdge*
locateFromEdge(const Vertex & v,const QuadEdge & startEdge) const138 QuadEdgeSubdivision::locateFromEdge(const Vertex& v,
139 const QuadEdge& startEdge) const
140 {
141 ::geos::ignore_unused_variable_warning(startEdge);
142
143 size_t iter = 0;
144 auto maxIter = quadEdges.size();
145
146 QuadEdge* e = startingEdges[0];
147
148 for(;;) {
149 ++iter;
150 /*
151 * So far it has always been the case that failure to locate indicates an
152 * invalid subdivision. So just fail completely. (An alternative would be
153 * to perform an exhaustive search for the containing triangle, but this
154 * would mask errors in the subdivision topology)
155 *
156 * This can also happen if two vertices are located very close together,
157 * since the orientation predicates may experience precision failures.
158 */
159 if(iter > maxIter) {
160 throw LocateFailureException("Could not locate vertex.");
161 }
162
163 if((v.equals(e->orig())) || (v.equals(e->dest()))) {
164 break;
165 }
166 else if(v.rightOf(*e)) {
167 e = &e->sym();
168 }
169 else if(!v.rightOf(e->oNext())) {
170 e = &e->oNext();
171 }
172 else if(!v.rightOf(e->dPrev())) {
173 e = &e->dPrev();
174 }
175 else {
176 // on edge or in triangle containing edge
177 break;
178 }
179 }
180 return e;
181 }
182
183 QuadEdge*
locate(const Coordinate & p0,const Coordinate & p1)184 QuadEdgeSubdivision::locate(const Coordinate& p0, const Coordinate& p1)
185 {
186 // find an edge containing one of the points
187 QuadEdge* e = locator->locate(Vertex(p0));
188 if(e == nullptr) {
189 return nullptr;
190 }
191
192 // normalize so that p0 is origin of base edge
193 QuadEdge* base = e;
194 if(e->dest().getCoordinate().equals2D(p0)) {
195 base = &e->sym();
196 }
197 // check all edges around origin of base edge
198 QuadEdge* locEdge = base;
199 do {
200 if(locEdge->dest().getCoordinate().equals2D(p1)) {
201 return locEdge;
202 }
203 locEdge = &locEdge->oNext();
204 }
205 while(locEdge != base);
206 return nullptr;
207 }
208
209 QuadEdge&
insertSite(const Vertex & v)210 QuadEdgeSubdivision::insertSite(const Vertex& v)
211 {
212 QuadEdge* e = locate(v);
213
214 if((v.equals(e->orig(), tolerance)) || (v.equals(e->dest(), tolerance))) {
215 return *e; // point already in subdivision.
216 }
217
218 // Connect the new point to the vertices of the containing
219 // triangle (or quadrilateral, if the new point fell on an
220 // existing edge.)
221 QuadEdge* base = &makeEdge(e->orig(), v);
222 QuadEdge::splice(*base, *e);
223 QuadEdge* startEdge = base;
224 do {
225 base = &connect(*e, base->sym());
226 e = &base->oPrev();
227 }
228 while(&e->lNext() != startEdge);
229
230 return *startEdge;
231 }
232
233 bool
isFrameEdge(const QuadEdge & e) const234 QuadEdgeSubdivision::isFrameEdge(const QuadEdge& e) const
235 {
236 if(isFrameVertex(e.orig()) || isFrameVertex(e.dest())) {
237 return true;
238 }
239 return false;
240 }
241
242 bool
isFrameBorderEdge(const QuadEdge & e) const243 QuadEdgeSubdivision::isFrameBorderEdge(const QuadEdge& e) const
244 {
245 // check other vertex of triangle to left of edge
246 Vertex vLeftTriOther = e.lNext().dest();
247 if(isFrameVertex(vLeftTriOther)) {
248 return true;
249 }
250 // check other vertex of triangle to right of edge
251 Vertex vRightTriOther = e.sym().lNext().dest();
252 if(isFrameVertex(vRightTriOther)) {
253 return true;
254 }
255
256 return false;
257 }
258
259 bool
isFrameVertex(const Vertex & v) const260 QuadEdgeSubdivision::isFrameVertex(const Vertex& v) const
261 {
262 if(v.equals(frameVertex[0])) {
263 return true;
264 }
265 if(v.equals(frameVertex[1])) {
266 return true;
267 }
268 if(v.equals(frameVertex[2])) {
269 return true;
270 }
271 return false;
272 }
273
274 bool
isOnEdge(const QuadEdge & e,const Coordinate & p) const275 QuadEdgeSubdivision::isOnEdge(const QuadEdge& e, const Coordinate& p) const
276 {
277 geom::LineSegment seg;
278 seg.setCoordinates(e.orig().getCoordinate(), e.dest().getCoordinate());
279 double dist = seg.distance(p);
280 // heuristic (hack?)
281 return dist < edgeCoincidenceTolerance;
282 }
283
284 bool
isVertexOfEdge(const QuadEdge & e,const Vertex & v) const285 QuadEdgeSubdivision::isVertexOfEdge(const QuadEdge& e, const Vertex& v) const
286 {
287 if((v.equals(e.orig(), tolerance)) || (v.equals(e.dest(), tolerance))) {
288 return true;
289 }
290 return false;
291 }
292
293 std::unique_ptr<QuadEdgeSubdivision::QuadEdgeList>
getPrimaryEdges(bool includeFrame)294 QuadEdgeSubdivision::getPrimaryEdges(bool includeFrame)
295 {
296 QuadEdgeList* edges = new QuadEdgeList();
297 QuadEdgeStack edgeStack;
298
299 edgeStack.push(startingEdges[0]);
300
301 prepareVisit();
302
303 while(!edgeStack.empty()) {
304 QuadEdge* edge = edgeStack.top();
305 edgeStack.pop();
306 if(!edge->isVisited()) {
307 QuadEdge* priQE = (QuadEdge*)&edge->getPrimary();
308
309 if(includeFrame || ! isFrameEdge(*priQE)) {
310 edges->push_back(priQE);
311 }
312
313 edgeStack.push(&edge->oNext());
314 edgeStack.push(&edge->sym().oNext());
315
316 edge->setVisited(true);
317 edge->sym().setVisited(true);
318 }
319 }
320 return std::unique_ptr<QuadEdgeList>(edges);
321 }
322
323 QuadEdge**
fetchTriangleToVisit(QuadEdge * edge,QuadEdgeStack & edgeStack,bool includeFrame)324 QuadEdgeSubdivision::fetchTriangleToVisit(QuadEdge* edge,
325 QuadEdgeStack& edgeStack, bool includeFrame)
326 {
327 QuadEdge* curr = edge;
328 int edgeCount = 0;
329 bool isFrame = false;
330 do {
331 triEdges[edgeCount] = curr;
332
333 if(!includeFrame && isFrameEdge(*curr)) {
334 isFrame = true;
335 }
336
337 // push sym edges to visit next
338 QuadEdge* sym = &curr->sym();
339 if (!sym->isVisited()) {
340 edgeStack.push(sym);
341 }
342
343 // mark this edge as visited
344 curr->setVisited(true);
345
346 edgeCount++;
347 curr = &curr->lNext();
348 }
349 while(curr != edge);
350
351 if(!includeFrame && isFrame) {
352 return nullptr;
353 }
354 return triEdges;
355 }
356
357 class
358 QuadEdgeSubdivision::TriangleCoordinatesVisitor : public TriangleVisitor {
359 private:
360 QuadEdgeSubdivision::TriList* triCoords;
361 CoordinateArraySequenceFactory coordSeqFact;
362
363 public:
TriangleCoordinatesVisitor(QuadEdgeSubdivision::TriList * p_triCoords)364 TriangleCoordinatesVisitor(QuadEdgeSubdivision::TriList* p_triCoords): triCoords(p_triCoords)
365 {
366 }
367
368 void
visit(QuadEdge * triEdges[3])369 visit(QuadEdge* triEdges[3]) override
370 {
371 auto coordSeq = coordSeqFact.create(4, 0);
372 for(size_t i = 0; i < 3; i++) {
373 Vertex v = triEdges[i]->orig();
374 coordSeq->setAt(v.getCoordinate(), i);
375 }
376 coordSeq->setAt(triEdges[0]->orig().getCoordinate(), 3);
377 triCoords->push_back(std::move(coordSeq));
378 }
379 };
380
381
382 class
383 QuadEdgeSubdivision::TriangleCircumcentreVisitor : public TriangleVisitor {
384 public:
385 void
visit(QuadEdge * triEdges[3])386 visit(QuadEdge* triEdges[3]) override
387 {
388 Triangle triangle(triEdges[0]->orig().getCoordinate(),
389 triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
390 Coordinate cc;
391
392 //TODO: identify heuristic to allow calling faster circumcentre() when possible
393 triangle.circumcentreDD(cc);
394
395 Vertex ccVertex(cc);
396
397 for(int i = 0 ; i < 3 ; i++) {
398 triEdges[i]->rot().setOrig(ccVertex);
399 }
400 }
401 };
402
403
404 void
getTriangleCoordinates(QuadEdgeSubdivision::TriList * triList,bool includeFrame)405 QuadEdgeSubdivision::getTriangleCoordinates(QuadEdgeSubdivision::TriList* triList, bool includeFrame)
406 {
407 TriangleCoordinatesVisitor visitor(triList);
408 visitTriangles(&visitor, includeFrame);
409 }
410
411 void
prepareVisit()412 QuadEdgeSubdivision::prepareVisit() {
413 if (!visit_state_clean) {
414 for (auto& qe : quadEdges) {
415 qe.setVisited(false);
416 }
417 }
418
419 visit_state_clean = false;
420 }
421
422 void
visitTriangles(TriangleVisitor * triVisitor,bool includeFrame)423 QuadEdgeSubdivision::visitTriangles(TriangleVisitor* triVisitor, bool includeFrame)
424 {
425 QuadEdgeStack edgeStack;
426 edgeStack.push(startingEdges[0]);
427
428 prepareVisit();
429
430 while(!edgeStack.empty()) {
431 QuadEdge* edge = edgeStack.top();
432 edgeStack.pop();
433 if(!edge->isVisited()) {
434 QuadEdge** p_triEdges = fetchTriangleToVisit(edge, edgeStack, includeFrame);
435 if(p_triEdges != nullptr) {
436 triVisitor->visit(p_triEdges);
437 }
438 }
439 }
440 }
441
442 std::unique_ptr<geom::MultiLineString>
getEdges(const geom::GeometryFactory & geomFact)443 QuadEdgeSubdivision::getEdges(const geom::GeometryFactory& geomFact)
444 {
445 std::unique_ptr<QuadEdgeList> p_quadEdges(getPrimaryEdges(false));
446 std::vector<std::unique_ptr<Geometry>> edges;
447 const CoordinateSequenceFactory* coordSeqFact = geomFact.getCoordinateSequenceFactory();
448
449 edges.reserve(p_quadEdges->size());
450 for(const QuadEdge* qe : *p_quadEdges) {
451 auto coordSeq = coordSeqFact->create(2);
452
453 coordSeq->setAt(qe->orig().getCoordinate(), 0);
454 coordSeq->setAt(qe->dest().getCoordinate(), 1);
455
456 edges.emplace_back(geomFact.createLineString(coordSeq.release()));
457 }
458
459 return geomFact.createMultiLineString(std::move(edges));
460 }
461
462 std::unique_ptr<GeometryCollection>
getTriangles(const GeometryFactory & geomFact)463 QuadEdgeSubdivision::getTriangles(const GeometryFactory& geomFact)
464 {
465 TriList triPtsList;
466 getTriangleCoordinates(&triPtsList, false);
467 std::vector<std::unique_ptr<Geometry>> tris;
468 tris.reserve(triPtsList.size());
469
470 for(auto& coordSeq : triPtsList) {
471 tris.push_back(
472 geomFact.createPolygon(geomFact.createLinearRing(std::move(coordSeq))));
473 }
474
475 return geomFact.createGeometryCollection(std::move(tris));
476 }
477
478
479 //Methods for VoronoiDiagram
480 std::unique_ptr<geom::GeometryCollection>
getVoronoiDiagram(const geom::GeometryFactory & geomFact)481 QuadEdgeSubdivision::getVoronoiDiagram(const geom::GeometryFactory& geomFact)
482 {
483 return geomFact.createGeometryCollection(getVoronoiCellPolygons(geomFact));
484 }
485
486 std::unique_ptr<geom::MultiLineString>
getVoronoiDiagramEdges(const geom::GeometryFactory & geomFact)487 QuadEdgeSubdivision::getVoronoiDiagramEdges(const geom::GeometryFactory& geomFact)
488 {
489 return geomFact.createMultiLineString(getVoronoiCellEdges(geomFact));
490 }
491
492 std::vector<std::unique_ptr<geom::Geometry>>
getVoronoiCellPolygons(const geom::GeometryFactory & geomFact)493 QuadEdgeSubdivision::getVoronoiCellPolygons(const geom::GeometryFactory& geomFact)
494 {
495 std::vector<std::unique_ptr<geom::Geometry>> cells;
496 TriangleCircumcentreVisitor tricircumVisitor;
497
498 visitTriangles(&tricircumVisitor, true);
499
500 std::unique_ptr<QuadEdgeSubdivision::QuadEdgeList> edges = getVertexUniqueEdges(false);
501
502 cells.reserve(edges->size());
503 for(const QuadEdge* qe : *edges) {
504 cells.push_back(getVoronoiCellPolygon(qe, geomFact));
505 }
506
507 return cells;
508 }
509
510 std::vector<std::unique_ptr<geom::Geometry>>
getVoronoiCellEdges(const geom::GeometryFactory & geomFact)511 QuadEdgeSubdivision::getVoronoiCellEdges(const geom::GeometryFactory& geomFact)
512 {
513 std::vector<std::unique_ptr<geom::Geometry>> cells;
514 TriangleCircumcentreVisitor tricircumVisitor;
515
516 visitTriangles((TriangleVisitor*) &tricircumVisitor, true);
517
518 std::unique_ptr<QuadEdgeSubdivision::QuadEdgeList> edges = getVertexUniqueEdges(false);
519 cells.reserve(edges->size());
520
521 for(const QuadEdge* qe : *edges) {
522 cells.push_back(getVoronoiCellEdge(qe, geomFact));
523 }
524
525 return cells;
526 }
527
528 std::unique_ptr<geom::Geometry>
getVoronoiCellPolygon(const QuadEdge * qe,const geom::GeometryFactory & geomFact)529 QuadEdgeSubdivision::getVoronoiCellPolygon(const QuadEdge* qe, const geom::GeometryFactory& geomFact)
530 {
531 std::vector<Coordinate> cellPts;
532
533 const QuadEdge* startQE = qe;
534 do {
535 const Coordinate& cc = qe->rot().orig().getCoordinate();
536 if(cellPts.empty() || cellPts.back() != cc) { // no duplicates
537 cellPts.push_back(cc);
538 }
539 qe = &qe->oPrev();
540
541 }
542 while(qe != startQE);
543
544 // Close the ring
545 if (cellPts.front() != cellPts.back()) {
546 cellPts.push_back(cellPts.front());
547 }
548 if (cellPts.size() < 4) {
549 cellPts.push_back(cellPts.back());
550 }
551
552 auto seq = geomFact.getCoordinateSequenceFactory()->create(std::move(cellPts));
553 std::unique_ptr<Geometry> cellPoly = geomFact.createPolygon(geomFact.createLinearRing(std::move(seq)));
554
555 // FIXME why is this returning a pointer to a local variable?
556 Vertex v = startQE->orig();
557 Coordinate c(0, 0);
558 c = v.getCoordinate();
559 cellPoly->setUserData(reinterpret_cast<void*>(&c));
560 return cellPoly;
561 }
562
563 std::unique_ptr<geom::Geometry>
getVoronoiCellEdge(const QuadEdge * qe,const geom::GeometryFactory & geomFact)564 QuadEdgeSubdivision::getVoronoiCellEdge(const QuadEdge* qe, const geom::GeometryFactory& geomFact)
565 {
566 std::vector<Coordinate> cellPts;
567
568 const QuadEdge* startQE = qe;
569 do {
570 const Coordinate& cc = qe->rot().orig().getCoordinate();
571 if(cellPts.empty() || cellPts.back() != cc) { // no duplicates
572 cellPts.push_back(cc);
573 }
574 qe = &qe->oPrev();
575
576 }
577 while(qe != startQE);
578
579 // Close the ring
580 if (cellPts.front() != cellPts.back()) {
581 cellPts.push_back(cellPts.front());
582 }
583
584 std::unique_ptr<geom::Geometry> cellEdge(
585 geomFact.createLineString(new geom::CoordinateArraySequence(std::move(cellPts))));
586
587 // FIXME why is this returning a pointer to a local variable?
588 Vertex v = startQE->orig();
589 Coordinate c(0, 0);
590 c = v.getCoordinate();
591 cellEdge->setUserData(reinterpret_cast<void*>(&c));
592 return cellEdge;
593 }
594
595 std::unique_ptr<QuadEdgeSubdivision::QuadEdgeList>
getVertexUniqueEdges(bool includeFrame)596 QuadEdgeSubdivision::getVertexUniqueEdges(bool includeFrame)
597 {
598 auto edges = detail::make_unique<QuadEdgeList>();
599 std::set<Vertex> visitedVertices; // TODO unordered_set of Vertex* ?
600
601 for(auto& quartet : quadEdges) {
602 QuadEdge* qe = &quartet.base();
603 const Vertex& v = qe->orig();
604
605 if(visitedVertices.find(v) == visitedVertices.end()) { //if v not found
606 visitedVertices.insert(v);
607
608 if(includeFrame || ! QuadEdgeSubdivision::isFrameVertex(v)) {
609 edges->push_back(qe);
610 }
611 }
612 QuadEdge* qd = &(qe->sym());
613 const Vertex& vd = qd->orig();
614
615 if(visitedVertices.find(vd) == visitedVertices.end()) {
616 visitedVertices.insert(vd);
617 if(includeFrame || ! QuadEdgeSubdivision::isFrameVertex(vd)) {
618 edges->push_back(qd);
619 }
620 }
621 }
622 return edges;
623 }
624
625 } //namespace geos.triangulate.quadedge
626 } //namespace geos.triangulate
627 } //namespace goes
628