CDT  1.4.2
C++ library for constrained Delaunay triangulation
 
Loading...
Searching...
No Matches
Triangulation.hpp
Go to the documentation of this file.
1/* This Source Code Form is subject to the terms of the Mozilla Public
2 * License, v. 2.0. If a copy of the MPL was not distributed with this
3 * file, You can obtain one at https://mozilla.org/MPL/2.0/. */
4
9
10#include "Triangulation.h"
11#include "portable_nth_element.hpp"
12
13#include <algorithm>
14#include <cassert>
15#include <cmath>
16#include <deque>
17#include <stdexcept>
18
19namespace CDT
20{
21
22typedef std::deque<TriInd> TriDeque;
23
24namespace detail
25{
26
27namespace defaults
28{
29
30const std::size_t nTargetVerts = 0;
32const VertexInsertionOrder::Enum vertexInsertionOrder =
34const IntersectingConstraintEdges::Enum intersectingEdgesStrategy =
36const float minDistToConstraintEdge(0);
37
38} // namespace defaults
39
40} // namespace detail
41
42template <typename T, typename TNearPointLocator>
44 : m_nTargetVerts(detail::defaults::nTargetVerts)
45 , m_superGeomType(detail::defaults::superGeomType)
46 , m_vertexInsertionOrder(detail::defaults::vertexInsertionOrder)
47 , m_intersectingEdgesStrategy(detail::defaults::intersectingEdgesStrategy)
48 , m_minDistToConstraintEdge(detail::defaults::minDistToConstraintEdge)
49#ifdef CDT_ENABLE_CALLBACK_HANDLER
50 , m_callbackHandler(NULL)
51#endif
52{}
53
54template <typename T, typename TNearPointLocator>
56 const VertexInsertionOrder::Enum vertexInsertionOrder)
57 : m_nTargetVerts(detail::defaults::nTargetVerts)
58 , m_superGeomType(detail::defaults::superGeomType)
59 , m_vertexInsertionOrder(vertexInsertionOrder)
60 , m_intersectingEdgesStrategy(detail::defaults::intersectingEdgesStrategy)
61 , m_minDistToConstraintEdge(detail::defaults::minDistToConstraintEdge)
62#ifdef CDT_ENABLE_CALLBACK_HANDLER
63 , m_callbackHandler(NULL)
64#endif
65{}
66
67template <typename T, typename TNearPointLocator>
69 const VertexInsertionOrder::Enum vertexInsertionOrder,
70 const IntersectingConstraintEdges::Enum intersectingEdgesStrategy,
71 const T minDistToConstraintEdge)
72 : m_nTargetVerts(detail::defaults::nTargetVerts)
73 , m_superGeomType(detail::defaults::superGeomType)
74 , m_vertexInsertionOrder(vertexInsertionOrder)
75 , m_intersectingEdgesStrategy(intersectingEdgesStrategy)
76 , m_minDistToConstraintEdge(minDistToConstraintEdge)
77#ifdef CDT_ENABLE_CALLBACK_HANDLER
78 , m_callbackHandler(NULL)
79#endif
80{}
81
82template <typename T, typename TNearPointLocator>
84 const VertexInsertionOrder::Enum vertexInsertionOrder,
85 const TNearPointLocator& nearPtLocator,
86 const IntersectingConstraintEdges::Enum intersectingEdgesStrategy,
87 const T minDistToConstraintEdge)
88 : m_nearPtLocator(nearPtLocator)
89 , m_nTargetVerts(detail::defaults::nTargetVerts)
90 , m_superGeomType(detail::defaults::superGeomType)
91 , m_vertexInsertionOrder(vertexInsertionOrder)
92 , m_intersectingEdgesStrategy(intersectingEdgesStrategy)
93 , m_minDistToConstraintEdge(minDistToConstraintEdge)
94#ifdef CDT_ENABLE_CALLBACK_HANDLER
95 , m_callbackHandler(NULL)
96#endif
97{}
98
99template <typename T, typename TNearPointLocator>
101{
102 if(m_superGeomType != SuperGeometryType::SuperTriangle)
103 return;
104 // find triangles adjacent to super-triangle's vertices
105 TriIndUSet toErase;
106 for(TriInd iT(0); iT < TriInd(triangles.size()); ++iT)
107 {
108 if(touchesSuperTriangle(triangles[iT]))
109 toErase.insert(iT);
110 }
111 finalizeTriangulation(toErase);
112}
113
114template <typename T, typename TNearPointLocator>
116{
117 assert(m_vertTris[0] != noNeighbor);
118 const std::stack<TriInd> seed(std::deque<TriInd>(1, m_vertTris[0]));
119 const TriIndUSet toErase = growToBoundary(seed);
120 finalizeTriangulation(toErase);
121}
122
123template <typename T, typename TNearPointLocator>
125{
126 const std::vector<LayerDepth> triDepths = calculateTriangleDepths();
127 TriIndUSet toErase;
128 toErase.reserve(triangles.size());
129 for(std::size_t iT = 0; iT != triangles.size(); ++iT)
130 {
131 if(triDepths[iT] % 2 == 0)
132 toErase.insert(static_cast<TriInd>(iT));
133 }
134 finalizeTriangulation(toErase);
135}
136
139{
140 return Edge(VertInd(e.v1() - 3), VertInd(e.v2() - 3));
141}
142
143template <typename T, typename TNearPointLocator>
145 const TriIndUSet& removedTriangles)
146{
147 if(removedTriangles.empty())
148 return;
149 // remove triangles and calculate triangle index mapping
150 TriIndUMap triIndMap;
151 for(TriInd iT(0), iTnew(0); iT < TriInd(triangles.size()); ++iT)
152 {
153 if(removedTriangles.count(iT))
154 continue;
155 triIndMap[iT] = iTnew;
156 triangles[iTnew] = triangles[iT];
157 iTnew++;
158 }
159 triangles.erase(triangles.end() - removedTriangles.size(), triangles.end());
160 // adjust triangles' neighbors
161 for(TriInd iT(0); iT < triangles.size(); ++iT)
162 {
163 Triangle& t = triangles[iT];
164 // update neighbors to account for removed triangles
165 NeighborsArr3& nn = t.neighbors;
166 for(NeighborsArr3::iterator n = nn.begin(); n != nn.end(); ++n)
167 {
168 if(removedTriangles.count(*n))
169 {
170 *n = noNeighbor;
171 }
172 else if(*n != noNeighbor)
173 {
174 *n = triIndMap[*n];
175 }
176 }
177 }
178}
179
180template <typename T, typename TNearPointLocator>
185
186template <typename T, typename TNearPointLocator>
188{
189 return m_vertTris;
190}
191
192template <typename T, typename TNearPointLocator>
193void Triangulation<T, TNearPointLocator>::finalizeTriangulation(
194 const TriIndUSet& removedTriangles)
195{
196 m_vertTris = TriIndVec();
197 // remove super-triangle
198 if(m_superGeomType == SuperGeometryType::SuperTriangle)
199 {
200 vertices.erase(
201 vertices.begin(), vertices.begin() + nSuperTriangleVertices);
202 // Edge re-mapping
203 { // fixed edges
204 EdgeUSet updatedFixedEdges;
205 typedef CDT::EdgeUSet::const_iterator It;
206 for(It e = fixedEdges.begin(); e != fixedEdges.end(); ++e)
207 {
208 updatedFixedEdges.insert(RemapNoSuperTriangle(*e));
209 }
210 fixedEdges = updatedFixedEdges;
211 }
212 { // overlap count
213 unordered_map<Edge, BoundaryOverlapCount> updatedOverlapCount;
214 typedef unordered_map<Edge, BoundaryOverlapCount>::const_iterator
215 It;
216 for(It it = overlapCount.begin(); it != overlapCount.end(); ++it)
217 {
218 updatedOverlapCount.insert(
219 std::make_pair(
220 RemapNoSuperTriangle(it->first), it->second));
221 }
222 overlapCount = updatedOverlapCount;
223 }
224 { // split edges mapping
225 unordered_map<Edge, EdgeVec> updatedPieceToOriginals;
226 typedef unordered_map<Edge, EdgeVec>::const_iterator It;
227 for(It it = pieceToOriginals.begin(); it != pieceToOriginals.end();
228 ++it)
229 {
230 EdgeVec ee = it->second;
231 for(EdgeVec::iterator eeIt = ee.begin(); eeIt != ee.end();
232 ++eeIt)
233 {
234 *eeIt = RemapNoSuperTriangle(*eeIt);
235 }
236 updatedPieceToOriginals.insert(
237 std::make_pair(RemapNoSuperTriangle(it->first), ee));
238 }
239 pieceToOriginals = updatedPieceToOriginals;
240 }
241 }
242 // remove other triangles
243 removeTriangles(removedTriangles);
244 // adjust triangle vertices: account for removed super-triangle
245 if(m_superGeomType == SuperGeometryType::SuperTriangle)
246 {
247 for(TriangleVec::iterator t = triangles.begin(); t != triangles.end();
248 ++t)
249 {
250 VerticesArr3& vv = t->vertices;
251 for(VerticesArr3::iterator v = vv.begin(); v != vv.end(); ++v)
252 {
253 *v -= nSuperTriangleVertices;
254 }
255 }
256 }
257}
258
259template <typename T, typename TNearPointLocator>
261{
262 m_nearPtLocator.initialize(vertices);
263 m_nTargetVerts = static_cast<IndexSizeType>(vertices.size());
264 m_superGeomType = SuperGeometryType::Custom;
265}
266
267template <typename T, typename TNearPointLocator>
268TriIndUSet Triangulation<T, TNearPointLocator>::growToBoundary(
269 std::stack<TriInd> seeds) const
270{
271 TriIndUSet traversed;
272 while(!seeds.empty())
273 {
274 const TriInd iT = seeds.top();
275 seeds.pop();
276 traversed.insert(iT);
277 const Triangle& t = triangles[iT];
278 for(Index i(0); i < Index(3); ++i)
279 {
280 const Edge opEdge(t.vertices[ccw(i)], t.vertices[cw(i)]);
281 if(fixedEdges.count(opEdge))
282 continue;
283 const TriInd iN = t.neighbors[opoNbr(i)];
284 if(iN != noNeighbor && traversed.count(iN) == 0)
285 seeds.push(iN);
286 }
287 }
288 return traversed;
289}
290
291template <typename T, typename TNearPointLocator>
292TriInd Triangulation<T, TNearPointLocator>::addTriangle(const Triangle& t)
293{
294 const TriInd iT(triangles.size());
295 triangles.push_back(t);
296 return iT;
297}
298
299template <typename T, typename TNearPointLocator>
300TriInd Triangulation<T, TNearPointLocator>::addTriangle()
301{
302 return addTriangle(Triangle());
303}
304
305template <typename T, typename TNearPointLocator>
307 const std::vector<Edge>& edges)
308{
309 insertEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2);
310}
311
312template <typename T, typename TNearPointLocator>
314 const std::vector<Edge>& edges)
315{
316 conformToEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2);
317}
318
319template <typename T, typename TNearPointLocator>
320void Triangulation<T, TNearPointLocator>::fixEdge(const Edge& edge)
321{
322 if(!fixedEdges.insert(edge).second)
323 {
324 ++overlapCount[edge]; // if edge is already fixed increment the counter
325 }
326}
327
328namespace detail
329{
330
331// add element to 'to' if not already in 'to'
332template <typename T, typename Allocator1>
333void insert_unique(std::vector<T, Allocator1>& to, const T& elem)
334{
335 if(std::find(to.begin(), to.end(), elem) == to.end())
336 {
337 to.push_back(elem);
338 }
339}
340
341// add elements of 'from' that are not present in 'to' to 'to'
342template <typename T, typename Allocator1, typename Allocator2>
343void insert_unique(
344 std::vector<T, Allocator1>& to,
345 const std::vector<T, Allocator2>& from)
346{
347 typedef typename std::vector<T, Allocator2>::const_iterator Cit;
348 to.reserve(to.size() + from.size());
349 for(Cit cit = from.begin(); cit != from.end(); ++cit)
350 {
351 insert_unique(to, *cit);
352 }
353}
354
355} // namespace detail
356
357template <typename T, typename TNearPointLocator>
358void Triangulation<T, TNearPointLocator>::splitFixedEdge(
359 const Edge& edge,
360 const VertInd iSplitVert)
361{
362 // split constraint (fixed) edge that already exists in triangulation
363 const Edge half1(edge.v1(), iSplitVert);
364 const Edge half2(iSplitVert, edge.v2());
365 // remove the edge that and add its halves
366 fixedEdges.erase(edge);
367 fixEdge(half1);
368 fixEdge(half2);
369 // maintain overlaps
370 typedef unordered_map<Edge, BoundaryOverlapCount>::const_iterator It;
371 const It overlapIt = overlapCount.find(edge);
372 if(overlapIt != overlapCount.end())
373 {
374 overlapCount[half1] += overlapIt->second;
375 overlapCount[half2] += overlapIt->second;
376 overlapCount.erase(overlapIt);
377 }
378 // maintain piece-to-original mapping
379 EdgeVec newOriginals(1, edge);
380 const unordered_map<Edge, EdgeVec>::const_iterator originalsIt =
381 pieceToOriginals.find(edge);
382 if(originalsIt != pieceToOriginals.end())
383 { // edge being split was split before: pass-through originals
384 newOriginals = originalsIt->second;
385 pieceToOriginals.erase(originalsIt);
386 }
387 detail::insert_unique(pieceToOriginals[half1], newOriginals);
388 detail::insert_unique(pieceToOriginals[half2], newOriginals);
389}
390
391template <typename T, typename TNearPointLocator>
392VertInd Triangulation<T, TNearPointLocator>::addSplitEdgeVertex(
393 const V2d<T>& splitVert,
394 const TriInd iT,
395 const TriInd iTopo)
396{
397 // add a new point on the edge that splits an edge in two
398 const VertInd iSplitVert = static_cast<VertInd>(vertices.size());
399 addNewVertex(splitVert, noNeighbor);
400
401#ifdef CDT_ENABLE_CALLBACK_HANDLER
402 if(m_callbackHandler)
403 {
404 m_callbackHandler->onAddVertexStart(
406 }
407#endif
408
409 std::stack<TriInd> triStack = insertVertexOnEdge(iSplitVert, iT, iTopo);
410 tryAddVertexToLocator(iSplitVert);
411 ensureDelaunayByEdgeFlips(iSplitVert, triStack);
412 return iSplitVert;
413}
414
415template <typename T, typename TNearPointLocator>
416VertInd Triangulation<T, TNearPointLocator>::splitFixedEdgeAt(
417 const Edge& edge,
418 const V2d<T>& splitVert,
419 const TriInd iT,
420 const TriInd iTopo)
421{
422 const VertInd iSplitVert = addSplitEdgeVertex(splitVert, iT, iTopo);
423 splitFixedEdge(edge, iSplitVert);
424 return iSplitVert;
425}
426
427template <typename T, typename TNearPointLocator>
428void Triangulation<T, TNearPointLocator>::fixEdge(
429 const Edge& edge,
430 const Edge& originalEdge)
431{
432 fixEdge(edge);
433 if(edge != originalEdge)
434 detail::insert_unique(pieceToOriginals[edge], originalEdge);
435}
436
437namespace detail
438{
439
440template <typename T>
441T lerp(const T& a, const T& b, const T t)
442{
443 return (T(1) - t) * a + t * b;
444}
445
446// Precondition: ab and cd intersect normally
447template <typename T>
448V2d<T> intersectionPosition(
449 const V2d<T>& a,
450 const V2d<T>& b,
451 const V2d<T>& c,
452 const V2d<T>& d)
453{
454 using namespace predicates::adaptive;
455
456 // note: for better accuracy we interpolate x and y separately
457 // on a segment with the shortest x/y-projection correspondingly
458 const T a_cd = orient2d(c.x, c.y, d.x, d.y, a.x, a.y);
459 const T b_cd = orient2d(c.x, c.y, d.x, d.y, b.x, b.y);
460 const T t_ab = a_cd / (a_cd - b_cd);
461
462 const T c_ab = orient2d(a.x, a.y, b.x, b.y, c.x, c.y);
463 const T d_ab = orient2d(a.x, a.y, b.x, b.y, d.x, d.y);
464 const T t_cd = c_ab / (c_ab - d_ab);
465
466 return V2d<T>(
467 std::fabs(a.x - b.x) < std::fabs(c.x - d.x) ? lerp(a.x, b.x, t_ab)
468 : lerp(c.x, d.x, t_cd),
469 std::fabs(a.y - b.y) < std::fabs(c.y - d.y) ? lerp(a.y, b.y, t_ab)
470 : lerp(c.y, d.y, t_cd));
471}
472
473} // namespace detail
474
475template <typename T, typename TNearPointLocator>
476void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
477 const Edge edge,
478 const Edge originalEdge,
479 EdgeVec& remaining,
480 std::vector<TriangulatePseudoPolygonTask>& tppIterations)
481{
482 const VertInd iA = edge.v1();
483 VertInd iB = edge.v2();
484 if(iA == iB) // edge connects a vertex to itself
485 return;
486
487 if(hasEdge(iA, iB))
488 {
489 fixEdge(edge, originalEdge);
490 return;
491 }
492
493 const V2d<T>& a = vertices[iA];
494 const V2d<T>& b = vertices[iB];
495 const T distanceTolerance =
496 m_minDistToConstraintEdge == T(0)
497 ? T(0)
498 : m_minDistToConstraintEdge * distance(a, b);
499
500 TriInd iT;
501 // Note: 'L' is left and 'R' is right of the inserted constraint edge
502 VertInd iVL, iVR;
503 tie(iT, iVL, iVR) = intersectedTriangle(iA, a, b, distanceTolerance);
504 // if one of the triangle vertices is on the edge, move edge start
505 if(iT == noNeighbor)
506 {
507 const Edge edgePart(iA, iVL);
508 fixEdge(edgePart, originalEdge);
509 remaining.push_back(Edge(iVL, iB));
510 return;
511 }
512 Triangle t = triangles[iT];
513 std::vector<TriInd> intersected(1, iT);
514 std::vector<VertInd> polyL, polyR;
515 polyL.reserve(2);
516 polyL.push_back(iA);
517 polyL.push_back(iVL);
518 polyR.reserve(2);
519 polyR.push_back(iA);
520 polyR.push_back(iVR);
521 unordered_map<Edge, TriInd> outerTris;
522 outerTris[Edge(iA, iVL)] = edgeNeighbor(t, iA, iVL);
523 outerTris[Edge(iA, iVR)] = edgeNeighbor(t, iA, iVR);
524 VertInd iV = iA;
525
526 while(!t.containsVertex(iB))
527 {
528 const TriInd iTopo = opposedTriangle(t, iV);
529 const Triangle& tOpo = triangles[iTopo];
530 const VertInd iVopo = opposedVertex(tOpo, iT);
531
532 switch(m_intersectingEdgesStrategy)
533 {
535 if(fixedEdges.count(Edge(iVL, iVR)))
536 {
537 // make sure to report original input edges in the exception
538 Edge e1 = originalEdge;
539 Edge e2 = Edge(iVL, iVR);
540 e2 = pieceToOriginals.count(e2)
541 ? pieceToOriginals.at(e2).front()
542 : e2;
543 // don't count super-triangle vertices
544 e1 = Edge(
545 VertInd(e1.v1() - m_nTargetVerts),
546 VertInd(e1.v2() - m_nTargetVerts));
547 e2 = Edge(
548 VertInd(e2.v1() - m_nTargetVerts),
549 VertInd(e2.v2() - m_nTargetVerts));
550 handleException(IntersectingConstraintsError(
551 e1,
552 pieceToOriginals.count(e2) ? pieceToOriginals.at(e2).front()
553 : e2,
554 CDT_SOURCE_LOCATION));
555 }
556 break;
558 if(!fixedEdges.count(Edge(iVL, iVR)))
559 break;
560 // split edge at the intersection of two constraint edges
561 const V2d<T> newV = detail::intersectionPosition(
562 vertices[iA], vertices[iB], vertices[iVL], vertices[iVR]);
563 const VertInd iNewVert =
564 splitFixedEdgeAt(Edge(iVL, iVR), newV, iT, iTopo);
565 // TODO: is it's possible to re-use pseudo-polygons
566 // for inserting [iA, iNewVert] edge half?
567 remaining.push_back(Edge(iA, iNewVert));
568 remaining.push_back(Edge(iNewVert, iB));
569 return;
570 }
572 assert(!fixedEdges.count(Edge(iVL, iVR)));
573 break;
574 }
575
576 const PtLineLocation::Enum loc =
577 locatePointLine(vertices[iVopo], a, b, distanceTolerance);
578 if(loc == PtLineLocation::Left)
579 {
580 const Edge e(polyL.back(), iVopo);
581 const TriInd outer = edgeNeighbor(tOpo, e.v1(), e.v2());
582 if(!outerTris.insert(std::make_pair(e, outer)).second)
583 outerTris.at(e) = noNeighbor; // hanging edge detected
584 polyL.push_back(iVopo);
585 iV = iVL;
586 iVL = iVopo;
587 }
588 else if(loc == PtLineLocation::Right)
589 {
590 const Edge e(polyR.back(), iVopo);
591 const TriInd outer = edgeNeighbor(tOpo, e.v1(), e.v2());
592 if(!outerTris.insert(std::make_pair(e, outer)).second)
593 outerTris.at(e) = noNeighbor; // hanging edge detected
594 polyR.push_back(iVopo);
595 iV = iVR;
596 iVR = iVopo;
597 }
598 else // encountered point on the edge
599 iB = iVopo;
600
601 intersected.push_back(iTopo);
602 iT = iTopo;
603 t = triangles[iT];
604 }
605 outerTris[Edge(polyL.back(), iB)] = edgeNeighbor(t, polyL.back(), iB);
606 outerTris[Edge(polyR.back(), iB)] = edgeNeighbor(t, polyR.back(), iB);
607 polyL.push_back(iB);
608 polyR.push_back(iB);
609
610 assert(!intersected.empty());
611 // make sure start/end vertices have a valid adjacent triangle
612 // that is not intersected by an edge
613 if(m_vertTris[iA] == intersected.front())
614 pivotVertexTriangleCW(iA);
615 if(m_vertTris[iB] == intersected.back())
616 pivotVertexTriangleCW(iB);
617
618 {
619#ifdef CDT_ENABLE_CALLBACK_HANDLER
620 if(m_callbackHandler)
621 {
622 m_callbackHandler->onReTriangulatePolygon(intersected);
623 }
624#endif
625
626 // Triangulate pseudo-polygons on both sides
627 std::reverse(polyR.begin(), polyR.end());
628
629 // note: intersected triangles are re-used for new triangles
630 // every triangulation of an n-gon has n − 2 triangles
631 // even if outer polygon has hanging edges it holds
632 assert(intersected.size() >= 2);
633 const TriInd iTL = intersected.back();
634 intersected.pop_back();
635 const TriInd iTR = intersected.back();
636 intersected.pop_back();
637
638 triangulatePseudoPolygon(
639 polyL, outerTris, iTL, iTR, intersected, tppIterations);
640 triangulatePseudoPolygon(
641 polyR, outerTris, iTR, iTL, intersected, tppIterations);
642 assert(intersected.empty());
643 }
644
645 if(iB != edge.v2()) // encountered point on the edge
646 {
647 // fix edge part
648 const Edge edgePart(iA, iB);
649 fixEdge(edgePart, originalEdge);
650 remaining.push_back(Edge(iB, edge.v2()));
651 return;
652 }
653 else
654 {
655 fixEdge(edge, originalEdge);
656 }
657}
658
659template <typename T, typename TNearPointLocator>
660void Triangulation<T, TNearPointLocator>::insertEdge(
661 Edge edge,
662 const Edge originalEdge,
663 EdgeVec& remaining,
664 std::vector<TriangulatePseudoPolygonTask>& tppIterations)
665{
666#ifdef CDT_ENABLE_CALLBACK_HANDLER
667 if(m_callbackHandler)
668 {
669 m_callbackHandler->onAddEdgeStart(edge);
670 }
671#endif
672
673 // use iteration over recursion to avoid stack overflows
674 remaining.clear();
675 remaining.push_back(edge);
676 while(!remaining.empty())
677 {
678 edge = remaining.back();
679 remaining.pop_back();
680 insertEdgeIteration(edge, originalEdge, remaining, tppIterations);
681 }
682}
683
684template <typename T, typename TNearPointLocator>
685void Triangulation<T, TNearPointLocator>::conformToEdgeIteration(
686 Edge edge,
687 const EdgeVec& originals,
688 BoundaryOverlapCount overlaps,
689 std::vector<ConformToEdgeTask>& remaining)
690{
691 const VertInd iA = edge.v1();
692 VertInd iB = edge.v2();
693 if(iA == iB) // edge connects a vertex to itself
694 return;
695
696 if(hasEdge(iA, iB))
697 {
698 fixEdge(edge);
699 if(overlaps > 0)
700 overlapCount[edge] = overlaps;
701 // avoid marking edge as a part of itself
702 if(!originals.empty() && edge != originals.front())
703 {
704 detail::insert_unique(pieceToOriginals[edge], originals);
705 }
706 return;
707 }
708
709 const V2d<T>& a = vertices[iA];
710 const V2d<T>& b = vertices[iB];
711 const T distanceTolerance =
712 m_minDistToConstraintEdge == T(0)
713 ? T(0)
714 : m_minDistToConstraintEdge * distance(a, b);
715 TriInd iT;
716 VertInd iVleft, iVright;
717 tie(iT, iVleft, iVright) = intersectedTriangle(iA, a, b, distanceTolerance);
718 // if one of the triangle vertices is on the edge, move edge start
719 if(iT == noNeighbor)
720 {
721 const Edge edgePart(iA, iVleft);
722 fixEdge(edgePart);
723 if(overlaps > 0)
724 overlapCount[edgePart] = overlaps;
725 detail::insert_unique(pieceToOriginals[edgePart], originals);
726#ifdef CDT_CXX11_IS_SUPPORTED
727 remaining.emplace_back(Edge(iVleft, iB), originals, overlaps);
728#else
729 remaining.push_back(make_tuple(Edge(iVleft, iB), originals, overlaps));
730#endif
731 return;
732 }
733
734 VertInd iV = iA;
735 Triangle t = triangles[iT];
736 while(std::find(t.vertices.begin(), t.vertices.end(), iB) ==
737 t.vertices.end())
738 {
739 const TriInd iTopo = opposedTriangle(t, iV);
740 const Triangle& tOpo = triangles[iTopo];
741 const VertInd iVopo = opposedVertex(tOpo, iT);
742 const V2d<T> vOpo = vertices[iVopo];
743
744 switch(m_intersectingEdgesStrategy)
745 {
747 if(fixedEdges.count(Edge(iVleft, iVright)))
748 {
749 // make sure to report original input edges in the exception
750 Edge e1 = pieceToOriginals.count(edge)
751 ? pieceToOriginals.at(edge).front()
752 : edge;
753 Edge e2(iVleft, iVright);
754 e2 = pieceToOriginals.count(e2)
755 ? pieceToOriginals.at(e2).front()
756 : e2;
757 // don't count super-triangle vertices
758 e1 = Edge(
759 VertInd(e1.v1() - m_nTargetVerts),
760 VertInd(e1.v2() - m_nTargetVerts));
761 e2 = Edge(
762 VertInd(e2.v1() - m_nTargetVerts),
763 VertInd(e2.v2() - m_nTargetVerts));
764 handleException(
765 IntersectingConstraintsError(e1, e2, CDT_SOURCE_LOCATION));
766 }
767 break;
769 if(!fixedEdges.count(Edge(iVleft, iVright)))
770 break;
771 // split edge at the intersection of two constraint edges
772 const V2d<T> newV = detail::intersectionPosition(
773 vertices[iA],
774 vertices[iB],
775 vertices[iVleft],
776 vertices[iVright]);
777 const VertInd iNewVert =
778 splitFixedEdgeAt(Edge(iVleft, iVright), newV, iT, iTopo);
779#ifdef CDT_CXX11_IS_SUPPORTED
780 remaining.emplace_back(Edge(iNewVert, iB), originals, overlaps);
781 remaining.emplace_back(Edge(iA, iNewVert), originals, overlaps);
782#else
783 remaining.push_back(
784 make_tuple(Edge(iNewVert, iB), originals, overlaps));
785 remaining.push_back(
786 make_tuple(Edge(iA, iNewVert), originals, overlaps));
787#endif
788 return;
789 }
791 assert(!fixedEdges.count(Edge(iVleft, iVright)));
792 break;
793 }
794
795 iT = iTopo;
796 t = triangles[iT];
797
798 const PtLineLocation::Enum loc =
799 locatePointLine(vOpo, a, b, distanceTolerance);
800 if(loc == PtLineLocation::Left)
801 {
802 iV = iVleft;
803 iVleft = iVopo;
804 }
805 else if(loc == PtLineLocation::Right)
806 {
807 iV = iVright;
808 iVright = iVopo;
809 }
810 else // encountered point on the edge
811 iB = iVopo;
812 }
813
814 // encountered one or more points on the edge: add remaining edge part
815 if(iB != edge.v2())
816 {
817#ifdef CDT_CXX11_IS_SUPPORTED
818 remaining.emplace_back(Edge(iB, edge.v2()), originals, overlaps);
819#else
820 remaining.push_back(
821 make_tuple(Edge(iB, edge.v2()), originals, overlaps));
822#endif
823 }
824
825 // add mid-point to triangulation
826 const VertInd iMid = static_cast<VertInd>(vertices.size());
827 const V2d<T>& start = vertices[iA];
828 const V2d<T>& end = vertices[iB];
829 addNewVertex(
830 V2d<T>((start.x + end.x) / T(2), (start.y + end.y) / T(2)), noNeighbor);
831
832#ifdef CDT_ENABLE_CALLBACK_HANDLER
833 if(m_callbackHandler)
834 {
835 m_callbackHandler->onAddVertexStart(
837 }
838#endif
839
840 const std::vector<Edge> flippedFixedEdges =
841 insertVertex_FlipFixedEdges(iMid);
842
843#ifdef CDT_CXX11_IS_SUPPORTED
844 remaining.emplace_back(Edge(iMid, iB), originals, overlaps);
845 remaining.emplace_back(Edge(iA, iMid), originals, overlaps);
846#else
847 remaining.push_back(make_tuple(Edge(iMid, iB), originals, overlaps));
848 remaining.push_back(make_tuple(Edge(iA, iMid), originals, overlaps));
849#endif
850
851 // re-introduce fixed edges that were flipped
852 // and make sure overlap count is preserved
853 for(std::vector<Edge>::const_iterator it = flippedFixedEdges.begin();
854 it != flippedFixedEdges.end();
855 ++it)
856 {
857 const Edge& flippedFixedEdge = *it;
858 fixedEdges.erase(flippedFixedEdge);
859
860 BoundaryOverlapCount prevOverlaps = 0;
861 const unordered_map<Edge, BoundaryOverlapCount>::const_iterator
862 overlapsIt = overlapCount.find(flippedFixedEdge);
863 if(overlapsIt != overlapCount.end())
864 {
865 prevOverlaps = overlapsIt->second;
866 overlapCount.erase(overlapsIt);
867 }
868 // override overlapping boundaries count when re-inserting an edge
869 EdgeVec prevOriginals(1, flippedFixedEdge);
870 const unordered_map<Edge, EdgeVec>::const_iterator originalsIt =
871 pieceToOriginals.find(flippedFixedEdge);
872 if(originalsIt != pieceToOriginals.end())
873 {
874 prevOriginals = originalsIt->second;
875 }
876#ifdef CDT_CXX11_IS_SUPPORTED
877 remaining.emplace_back(flippedFixedEdge, prevOriginals, prevOverlaps);
878#else
879 remaining.push_back(
880 make_tuple(flippedFixedEdge, prevOriginals, prevOverlaps));
881#endif
882 }
883}
884
885template <typename T, typename TNearPointLocator>
886void Triangulation<T, TNearPointLocator>::conformToEdge(
887 Edge edge,
888 EdgeVec originals,
889 BoundaryOverlapCount overlaps,
890 std::vector<ConformToEdgeTask>& remaining)
891{
892#ifdef CDT_ENABLE_CALLBACK_HANDLER
893 if(m_callbackHandler)
894 {
895 m_callbackHandler->onAddEdgeStart(edge);
896 }
897#endif
898
899 // use iteration over recursion to avoid stack overflows
900 remaining.clear();
901#ifdef CDT_CXX11_IS_SUPPORTED
902 remaining.emplace_back(edge, originals, overlaps);
903#else
904 remaining.push_back(make_tuple(edge, originals, overlaps));
905#endif
906 while(!remaining.empty())
907 {
908 tie(edge, originals, overlaps) = remaining.back();
909 remaining.pop_back();
910 conformToEdgeIteration(edge, originals, overlaps, remaining);
911 }
912}
913
924template <typename T, typename TNearPointLocator>
925tuple<TriInd, VertInd, VertInd>
926Triangulation<T, TNearPointLocator>::intersectedTriangle(
927 const VertInd iA,
928 const V2d<T>& a,
929 const V2d<T>& b,
930 const T orientationTolerance) const
931{
932 const TriInd startTri = m_vertTris[iA];
933 TriInd iT = startTri;
934 do
935 {
936 const Triangle t = triangles[iT];
937 const Index i = vertexInd(t.vertices, iA);
938 const VertInd iP2 = t.vertices[ccw(i)];
939 const T orientP2 = orient2D(vertices[iP2], a, b);
940 const PtLineLocation::Enum locP2 = classifyOrientation(orientP2);
941 if(locP2 == PtLineLocation::Right)
942 {
943 const VertInd iP1 = t.vertices[cw(i)];
944 const T orientP1 = orient2D(vertices[iP1], a, b);
945 const PtLineLocation::Enum locP1 = classifyOrientation(orientP1);
946 if(locP1 == PtLineLocation::OnLine)
947 {
948 return make_tuple(noNeighbor, iP1, iP1);
949 }
950 if(locP1 == PtLineLocation::Left)
951 {
952 if(orientationTolerance)
953 {
954 T closestOrient;
955 VertInd iClosestP;
956 if(std::abs(orientP1) <= std::abs(orientP2))
957 {
958 closestOrient = orientP1;
959 iClosestP = iP1;
960 }
961 else
962 {
963 closestOrient = orientP2;
964 iClosestP = iP2;
965 }
967 closestOrient, orientationTolerance) ==
968 PtLineLocation::OnLine)
969 {
970 return make_tuple(noNeighbor, iClosestP, iClosestP);
971 }
972 }
973 return make_tuple(iT, iP1, iP2);
974 }
975 }
976 iT = t.next(iA).first;
977 } while(iT != startTri);
978
979 handleException(Error(
980 "Could not find vertex triangle intersected by an edge.",
981 CDT_SOURCE_LOCATION));
982 return make_tuple(noNeighbor, noVertex, noVertex);
983}
984
985template <typename T, typename TNearPointLocator>
986void Triangulation<T, TNearPointLocator>::addSuperTriangle(const Box2d<T>& box)
987{
988 m_nTargetVerts = nSuperTriangleVertices;
989 m_superGeomType = SuperGeometryType::SuperTriangle;
990
991 const V2d<T> center = {
992 (box.min.x + box.max.x) / T(2), (box.min.y + box.max.y) / T(2)};
993 const T w = box.max.x - box.min.x;
994 const T h = box.max.y - box.min.y;
995 T r = std::max(w, h); // incircle radius upper bound
996
997 // Note: make sure radius is big enough. Constants chosen experimentally.
998 // - for tiny bounding boxes: use 1.0 as the smallest radius
999 // - multiply radius by 2.0 for extra safety margin
1000 r = std::max(T(2) * r, T(1));
1001
1002 // Note: for very large floating point numbers rounding can lead to wrong
1003 // super-triangle coordinates. This is a very rare corner-case so the
1004 // handling is very primitive.
1005 { // note: '<=' means '==' but avoids the warning
1006 while(center.y <= center.y - r)
1007 r *= T(2);
1008 }
1009
1010 const T R = T(2) * r; // excircle radius
1011 const T cos_30_deg = T(0.8660254037844386); // note: (std::sqrt(3.0) / 2.0)
1012 const T shiftX = R * cos_30_deg;
1013 const V2d<T> posV1 = {center.x - shiftX, center.y - r};
1014 const V2d<T> posV2 = {center.x + shiftX, center.y - r};
1015 const V2d<T> posV3 = {center.x, center.y + R};
1016 addNewVertex(posV1, TriInd(0));
1017 addNewVertex(posV2, TriInd(0));
1018 addNewVertex(posV3, TriInd(0));
1019
1020#ifdef CDT_ENABLE_CALLBACK_HANDLER
1021 if(m_callbackHandler)
1022 {
1023 m_callbackHandler->onAddSuperTriangle();
1024 }
1025#endif
1026
1027 addTriangle(
1028 Triangle(arr3(VertInd(0), VertInd(1), VertInd(2)), arr3(noNeighbor)));
1029
1030 if(m_vertexInsertionOrder != VertexInsertionOrder::Auto)
1031 {
1032 m_nearPtLocator.initialize(vertices);
1033 }
1034}
1035
1036template <typename T, typename TNearPointLocator>
1037void Triangulation<T, TNearPointLocator>::addNewVertex(
1038 const V2d<T>& pos,
1039 const TriInd iT)
1040{
1041 vertices.push_back(pos);
1042 m_vertTris.push_back(iT);
1043}
1044
1045template <typename T, typename TNearPointLocator>
1046std::vector<Edge>
1047Triangulation<T, TNearPointLocator>::insertVertex_FlipFixedEdges(
1048 const VertInd iV1)
1049{
1050 std::vector<Edge> flippedFixedEdges;
1051
1052 const V2d<T>& v1 = vertices[iV1];
1053 const VertInd startVertex = m_nearPtLocator.nearPoint(v1, vertices);
1054 array<TriInd, 2> trisAt = walkingSearchTrianglesAt(iV1, startVertex);
1055 std::stack<TriInd> triStack =
1056 trisAt[1] == noNeighbor ? insertVertexInsideTriangle(iV1, trisAt[0])
1057 : insertVertexOnEdge(iV1, trisAt[0], trisAt[1]);
1058
1059 TriInd iTopo, n1, n2, n3, n4;
1060 VertInd iV2, iV3, iV4;
1061 while(!triStack.empty())
1062 {
1063 const TriInd iT = triStack.top();
1064 triStack.pop();
1065
1066 edgeFlipInfo(iT, iV1, iTopo, iV2, iV3, iV4, n1, n2, n3, n4);
1067 if(iTopo != noNeighbor && isFlipNeeded(iV1, iV2, iV3, iV4))
1068 {
1069 // if flipped edge is fixed, remember it
1070 const Edge flippedEdge(iV2, iV4);
1071 if(!fixedEdges.empty() &&
1072 fixedEdges.find(flippedEdge) != fixedEdges.end())
1073 {
1074 flippedFixedEdges.push_back(flippedEdge);
1075 }
1076
1077 flipEdge(iT, iTopo, iV1, iV2, iV3, iV4, n1, n2, n3, n4);
1078 triStack.push(iT);
1079 triStack.push(iTopo);
1080 }
1081 }
1082
1083 tryAddVertexToLocator(iV1);
1084 return flippedFixedEdges;
1085}
1086
1087template <typename T, typename TNearPointLocator>
1088void Triangulation<T, TNearPointLocator>::insertVertex(
1089 const VertInd iVert,
1090 const VertInd walkStart)
1091{
1092#ifdef CDT_ENABLE_CALLBACK_HANDLER
1093 if(m_callbackHandler)
1094 {
1095 m_callbackHandler->onAddVertexStart(iVert, AddVertexType::UserInput);
1096 }
1097#endif
1098
1099 const array<TriInd, 2> trisAt = walkingSearchTrianglesAt(iVert, walkStart);
1100 std::stack<TriInd> triStack =
1101 trisAt[1] == noNeighbor
1102 ? insertVertexInsideTriangle(iVert, trisAt[0])
1103 : insertVertexOnEdge(iVert, trisAt[0], trisAt[1]);
1104 ensureDelaunayByEdgeFlips(iVert, triStack);
1105}
1106
1107template <typename T, typename TNearPointLocator>
1108void Triangulation<T, TNearPointLocator>::insertVertex(const VertInd iVert)
1109{
1110 const V2d<T>& v = vertices[iVert];
1111 const VertInd walkStart = m_nearPtLocator.nearPoint(v, vertices);
1112 insertVertex(iVert, walkStart);
1113 tryAddVertexToLocator(iVert);
1114}
1115
1116template <typename T, typename TNearPointLocator>
1117void Triangulation<T, TNearPointLocator>::ensureDelaunayByEdgeFlips(
1118 const VertInd iV1,
1119 std::stack<TriInd>& triStack)
1120{
1121 TriInd iTopo, n1, n2, n3, n4;
1122 VertInd iV2, iV3, iV4;
1123 while(!triStack.empty())
1124 {
1125 const TriInd iT = triStack.top();
1126 triStack.pop();
1127
1128 edgeFlipInfo(iT, iV1, iTopo, iV2, iV3, iV4, n1, n2, n3, n4);
1129 if(iTopo != noNeighbor && isFlipNeeded(iV1, iV2, iV3, iV4))
1130 {
1131 flipEdge(iT, iTopo, iV1, iV2, iV3, iV4, n1, n2, n3, n4);
1132 triStack.push(iT);
1133 triStack.push(iTopo);
1134 }
1135 }
1136}
1137
1138/*
1139 * v4 original edge: (v1, v3)
1140 * /|\ flip-candidate edge: (v, v2)
1141 * / | \
1142 * n3 / | \ n4
1143 * / | \
1144 * new vertex--> v1 T | Topo v3
1145 * \ | /
1146 * n1 \ | / n2
1147 * \ | /
1148 * \|/
1149 * v2
1150 */
1151template <typename T, typename TNearPointLocator>
1152void Triangulation<T, TNearPointLocator>::edgeFlipInfo(
1153 const TriInd iT,
1154 const VertInd iV1,
1155 TriInd& iTopo,
1156 VertInd& iV2,
1157 VertInd& iV3,
1158 VertInd& iV4,
1159 TriInd& n1,
1160 TriInd& n2,
1161 TriInd& n3,
1162 TriInd& n4)
1163{
1164 /* v[2]
1165 / \
1166 n[2]/ \n[1]
1167 /_____\
1168 v[0] n[0] v[1] */
1169 const Triangle& t = triangles[iT];
1170 if(t.vertices[0] == iV1)
1171 {
1172 iV2 = t.vertices[1];
1173 iV4 = t.vertices[2];
1174 n1 = t.neighbors[0];
1175 n3 = t.neighbors[2];
1176 iTopo = t.neighbors[1];
1177 }
1178 else if(t.vertices[1] == iV1)
1179 {
1180 iV2 = t.vertices[2];
1181 iV4 = t.vertices[0];
1182 n1 = t.neighbors[1];
1183 n3 = t.neighbors[0];
1184 iTopo = t.neighbors[2];
1185 }
1186 else
1187 {
1188 iV2 = t.vertices[0];
1189 iV4 = t.vertices[1];
1190 n1 = t.neighbors[2];
1191 n3 = t.neighbors[1];
1192 iTopo = t.neighbors[0];
1193 }
1194 if(iTopo == noNeighbor)
1195 return;
1196 const Triangle& tOpo = triangles[iTopo];
1197 if(tOpo.neighbors[0] == iT)
1198 {
1199 iV3 = tOpo.vertices[2];
1200 n2 = tOpo.neighbors[1];
1201 n4 = tOpo.neighbors[2];
1202 }
1203 else if(tOpo.neighbors[1] == iT)
1204 {
1205 iV3 = tOpo.vertices[0];
1206 n2 = tOpo.neighbors[2];
1207 n4 = tOpo.neighbors[0];
1208 }
1209 else
1210 {
1211 iV3 = tOpo.vertices[1];
1212 n2 = tOpo.neighbors[0];
1213 n4 = tOpo.neighbors[1];
1214 }
1215}
1216
1227/*
1228 * v4 original edge: (v2, v4)
1229 * /|\ flip-candidate edge: (v1, v3)
1230 * / | \
1231 * / | \
1232 * / | \
1233 * new vertex--> v1 | v3
1234 * \ | /
1235 * \ | /
1236 * \ | /
1237 * \|/
1238 * v2
1239 */
1240template <typename T, typename TNearPointLocator>
1241bool Triangulation<T, TNearPointLocator>::isFlipNeeded(
1242 const VertInd iV1,
1243 const VertInd iV2,
1244 const VertInd iV3,
1245 const VertInd iV4) const
1246{
1247 if(fixedEdges.count(Edge(iV2, iV4)))
1248 return false; // flip not needed if the original edge is fixed
1249 const V2d<T>& v1 = vertices[iV1];
1250 const V2d<T>& v2 = vertices[iV2];
1251 const V2d<T>& v3 = vertices[iV3];
1252 const V2d<T>& v4 = vertices[iV4];
1253 if(m_superGeomType == SuperGeometryType::SuperTriangle)
1254 {
1255 // If flip-candidate edge touches super-triangle in-circumference
1256 // test has to be replaced with orient2d test against the line
1257 // formed by two non-artificial vertices (that don't belong to
1258 // super-triangle)
1259 if(iV1 < 3) // flip-candidate edge touches super-triangle
1260 {
1261 // does original edge also touch super-triangle?
1262 if(iV2 < 3)
1263 return locatePointLine(v2, v3, v4) ==
1264 locatePointLine(v1, v3, v4);
1265 if(iV4 < 3)
1266 return locatePointLine(v4, v2, v3) ==
1267 locatePointLine(v1, v2, v3);
1268 return false; // original edge does not touch super-triangle
1269 }
1270 if(iV3 < 3) // flip-candidate edge touches super-triangle
1271 {
1272 // does original edge also touch super-triangle?
1273 if(iV2 < 3)
1274 {
1275 return locatePointLine(v2, v1, v4) ==
1276 locatePointLine(v3, v1, v4);
1277 }
1278 if(iV4 < 3)
1279 {
1280 return locatePointLine(v4, v2, v1) ==
1281 locatePointLine(v3, v2, v1);
1282 }
1283 return false; // original edge does not touch super-triangle
1284 }
1285 // flip-candidate edge does not touch super-triangle
1286 if(iV2 < 3)
1287 return locatePointLine(v2, v3, v4) == locatePointLine(v1, v3, v4);
1288 if(iV4 < 3)
1289 return locatePointLine(v4, v2, v3) == locatePointLine(v1, v2, v3);
1290 }
1291 return isInCircumcircle(v1, v2, v3, v4);
1292}
1293
1294/* Flip edge between T and Topo:
1295 *
1296 * v4 | - old edge
1297 * /|\ ~ - new edge
1298 * / | \
1299 * n3 / T' \ n4
1300 * / | \
1301 * / | \
1302 * T -> v1~~~~~~~~~v3 <- Topo
1303 * \ | /
1304 * \ | /
1305 * n1 \Topo'/ n2
1306 * \ | /
1307 * \|/
1308 * v2
1309 */
1310template <typename T, typename TNearPointLocator>
1312 const TriInd iT,
1313 const TriInd iTopo)
1314{
1315#ifdef CDT_ENABLE_CALLBACK_HANDLER
1316 if(m_callbackHandler)
1317 {
1318 m_callbackHandler->onFlipEdge(iT, iTopo);
1319 }
1320#endif
1321
1322 Triangle& t = triangles[iT];
1323 Triangle& tOpo = triangles[iTopo];
1324 const array<TriInd, 3>& triNs = t.neighbors;
1325 const array<TriInd, 3>& triOpoNs = tOpo.neighbors;
1326 const array<VertInd, 3>& triVs = t.vertices;
1327 const array<VertInd, 3>& triOpoVs = tOpo.vertices;
1328 // find vertices and neighbors
1329 Index i = opposedVertexInd(t.neighbors, iTopo);
1330 const VertInd v1 = triVs[i];
1331 const VertInd v2 = triVs[ccw(i)];
1332 const TriInd n1 = triNs[i];
1333 const TriInd n3 = triNs[cw(i)];
1334 i = opposedVertexInd(tOpo.neighbors, iT);
1335 const VertInd v3 = triOpoVs[i];
1336 const VertInd v4 = triOpoVs[ccw(i)];
1337 const TriInd n4 = triOpoNs[i];
1338 const TriInd n2 = triOpoNs[cw(i)];
1339 // change vertices and neighbors
1340 t = Triangle(arr3(v4, v1, v3), arr3(n3, iTopo, n4));
1341 tOpo = Triangle(arr3(v2, v3, v1), arr3(n2, iT, n1));
1342 // adjust neighboring triangles and vertices
1343 changeNeighbor(n1, iT, iTopo);
1344 changeNeighbor(n4, iTopo, iT);
1345 // only adjust adjacent triangles if triangulation is not finalized:
1346 // can happen when called from outside on an already finalized
1347 // triangulation
1348 if(!isFinalized())
1349 {
1350 setAdjacentTriangle(v4, iT);
1351 setAdjacentTriangle(v2, iTopo);
1352 }
1353}
1354
1355/* Flip edge between T and Topo:
1356 *
1357 * v4 | - old edge
1358 * /|\ ~ - new edge
1359 * / | \
1360 * n3 / T' \ n4
1361 * / | \
1362 * / | \
1363 * T -> v1 ~~~~~~~~ v3 <- Topo
1364 * \ | /
1365 * \ | /
1366 * n1 \Topo'/ n2
1367 * \ | /
1368 * \|/
1369 * v2
1370 */
1371template <typename T, typename TNearPointLocator>
1373 const TriInd iT,
1374 const TriInd iTopo,
1375 const VertInd v1,
1376 const VertInd v2,
1377 const VertInd v3,
1378 const VertInd v4,
1379 const TriInd n1,
1380 const TriInd n2,
1381 const TriInd n3,
1382 const TriInd n4)
1383{
1384#ifdef CDT_ENABLE_CALLBACK_HANDLER
1385 if(m_callbackHandler)
1386 {
1387 m_callbackHandler->onFlipEdge(iT, iTopo);
1388 }
1389#endif
1390
1391 // change vertices and neighbors
1392 triangles[iT] = Triangle(arr3(v4, v1, v3), arr3(n3, iTopo, n4));
1393 triangles[iTopo] = Triangle(arr3(v2, v3, v1), arr3(n2, iT, n1));
1394 // adjust neighboring triangles and vertices
1395 changeNeighbor(n1, iT, iTopo);
1396 changeNeighbor(n4, iTopo, iT);
1397 // only adjust adjacent triangles if triangulation is not finalized:
1398 // can happen when called from outside on an already finalized
1399 // triangulation
1400 if(!isFinalized())
1401 {
1402 setAdjacentTriangle(v4, iT);
1403 setAdjacentTriangle(v2, iTopo);
1404 }
1405}
1406
1407/* Insert point into triangle: split into 3 triangles:
1408 * - create 2 new triangles
1409 * - re-use old triangle for the 3rd
1410 * v3
1411 * / | \
1412 * / | \ <-- original triangle (t)
1413 * / | \
1414 * n3 / | \ n2
1415 * /newT2|newT1\
1416 * / v \
1417 * / __/ \__ \
1418 * / __/ \__ \
1419 * / _/ t' \_ \
1420 * v1 ___________________ v2
1421 * n1
1422 */
1423template <typename T, typename TNearPointLocator>
1424std::stack<TriInd>
1425Triangulation<T, TNearPointLocator>::insertVertexInsideTriangle(
1426 VertInd v,
1427 TriInd iT)
1428{
1429 const TriInd iNewT1 = addTriangle();
1430 const TriInd iNewT2 = addTriangle();
1431
1432#ifdef CDT_ENABLE_CALLBACK_HANDLER
1433 if(m_callbackHandler)
1434 {
1435 m_callbackHandler->onInsertVertexInsideTriangle(iT, iNewT1, iNewT2);
1436 }
1437#endif
1438
1439 Triangle& t = triangles[iT];
1440 const array<VertInd, 3> vv = t.vertices;
1441 const array<TriInd, 3> nn = t.neighbors;
1442 const VertInd v1 = vv[0], v2 = vv[1], v3 = vv[2];
1443 const TriInd n1 = nn[0], n2 = nn[1], n3 = nn[2];
1444 // make two new triangles and convert current triangle to 3rd new
1445 // triangle
1446 triangles[iNewT1] = Triangle(arr3(v2, v3, v), arr3(n2, iNewT2, iT));
1447 triangles[iNewT2] = Triangle(arr3(v3, v1, v), arr3(n3, iT, iNewT1));
1448 t = Triangle(arr3(v1, v2, v), arr3(n1, iNewT1, iNewT2));
1449 // adjust adjacent triangles
1450 setAdjacentTriangle(v, iT);
1451 setAdjacentTriangle(v3, iNewT1);
1452 // change triangle neighbor's neighbors to new triangles
1453 changeNeighbor(n2, iT, iNewT1);
1454 changeNeighbor(n3, iT, iNewT2);
1455 // return newly added triangles
1456 std::stack<TriInd> newTriangles;
1457 newTriangles.push(iT);
1458 newTriangles.push(iNewT1);
1459 newTriangles.push(iNewT2);
1460 return newTriangles;
1461}
1462
1463/* Inserting a point on the edge between two triangles
1464 * T1 (top) v1
1465 * /|\
1466 * n1 / | \ n4
1467 * / | \
1468 * / T1' | Tnew1\
1469 * v2-------v-------v4
1470 * \ T2' | Tnew2/
1471 * \ | /
1472 * n2 \ | / n3
1473 * \|/
1474 * T2 (bottom) v3
1475 */
1476template <typename T, typename TNearPointLocator>
1477std::stack<TriInd> Triangulation<T, TNearPointLocator>::insertVertexOnEdge(
1478 VertInd v,
1479 TriInd iT1,
1480 TriInd iT2)
1481{
1482 const TriInd iTnew1 = addTriangle();
1483 const TriInd iTnew2 = addTriangle();
1484
1485#ifdef CDT_ENABLE_CALLBACK_HANDLER
1486 if(m_callbackHandler)
1487 {
1488 m_callbackHandler->onInsertVertexOnEdge(iT1, iT2, iTnew1, iTnew2);
1489 }
1490#endif
1491
1492 Triangle& t1 = triangles[iT1];
1493 Triangle& t2 = triangles[iT2];
1494 Index i = opposedVertexInd(t1.neighbors, iT2);
1495 const VertInd v1 = t1.vertices[i];
1496 const VertInd v2 = t1.vertices[ccw(i)];
1497 const TriInd n1 = t1.neighbors[i];
1498 const TriInd n4 = t1.neighbors[cw(i)];
1499 i = opposedVertexInd(t2.neighbors, iT1);
1500 const VertInd v3 = t2.vertices[i];
1501 const VertInd v4 = t2.vertices[ccw(i)];
1502 const TriInd n3 = t2.neighbors[i];
1503 const TriInd n2 = t2.neighbors[cw(i)];
1504 // add new triangles and change existing ones
1505 t1 = Triangle(arr3(v, v1, v2), arr3(iTnew1, n1, iT2));
1506 t2 = Triangle(arr3(v, v2, v3), arr3(iT1, n2, iTnew2));
1507 triangles[iTnew1] = Triangle(arr3(v, v4, v1), arr3(iTnew2, n4, iT1));
1508 triangles[iTnew2] = Triangle(arr3(v, v3, v4), arr3(iT2, n3, iTnew1));
1509 // adjust adjacent triangles
1510 setAdjacentTriangle(v, iT1);
1511 setAdjacentTriangle(v4, iTnew1);
1512 // adjust neighboring triangles and vertices
1513 changeNeighbor(n4, iT1, iTnew1);
1514 changeNeighbor(n3, iT2, iTnew2);
1515 // return newly added triangles
1516 std::stack<TriInd> newTriangles;
1517 newTriangles.push(iT1);
1518 newTriangles.push(iTnew2);
1519 newTriangles.push(iT2);
1520 newTriangles.push(iTnew1);
1521 return newTriangles;
1522}
1523
1524template <typename T, typename TNearPointLocator>
1525array<TriInd, 2>
1526Triangulation<T, TNearPointLocator>::trianglesAt(const V2d<T>& pos) const
1527{
1528 array<TriInd, 2> out = {noNeighbor, noNeighbor};
1529 for(TriInd i = TriInd(0); i < TriInd(triangles.size()); ++i)
1530 {
1531 const Triangle& t = triangles[i];
1532 const V2d<T>& v1 = vertices[t.vertices[0]];
1533 const V2d<T>& v2 = vertices[t.vertices[1]];
1534 const V2d<T>& v3 = vertices[t.vertices[2]];
1535 const PtTriLocation::Enum loc = locatePointTriangle(pos, v1, v2, v3);
1536 if(loc == PtTriLocation::Outside)
1537 continue;
1538 out[0] = i;
1539 if(isOnEdge(loc))
1540 out[1] = t.neighbors[edgeNeighbor(loc)];
1541 return out;
1542 }
1543 handleException(
1544 Error("No triangle was found at position", CDT_SOURCE_LOCATION));
1545 return out;
1546}
1547
1548template <typename T, typename TNearPointLocator>
1549TriInd Triangulation<T, TNearPointLocator>::walkTriangles(
1550 const VertInd startVertex,
1551 const V2d<T>& pos) const
1552{
1553 // begin walk in search of triangle at pos
1554 TriInd currTri = m_vertTris[startVertex];
1555 bool found = false;
1556 detail::SplitMix64RandGen prng;
1557 while(!found)
1558 {
1559 const Triangle& t = triangles[currTri];
1560 found = true;
1561 // stochastic offset to randomize which edge we check first
1562 const Index offset(prng() % 3);
1563 for(Index i_(0); i_ < Index(3); ++i_)
1564 {
1565 const Index i((i_ + offset) % 3);
1566 const V2d<T>& vStart = vertices[t.vertices[i]];
1567 const V2d<T>& vEnd = vertices[t.vertices[ccw(i)]];
1568 const PtLineLocation::Enum edgeCheck =
1569 locatePointLine(pos, vStart, vEnd);
1570 const TriInd iN = t.neighbors[i];
1571 if(edgeCheck == PtLineLocation::Right && iN != noNeighbor)
1572 {
1573 found = false;
1574 currTri = iN;
1575 break;
1576 }
1577 }
1578 }
1579 return currTri;
1580}
1581
1582template <typename T, typename TNearPointLocator>
1583array<TriInd, 2> Triangulation<T, TNearPointLocator>::walkingSearchTrianglesAt(
1584 const VertInd iV,
1585 const VertInd startVertex) const
1586{
1587 const V2d<T> v = vertices[iV];
1588 array<TriInd, 2> out = {noNeighbor, noNeighbor};
1589 const TriInd iT = walkTriangles(startVertex, v);
1590 // Finished walk, locate point in current triangle
1591 const Triangle& t = triangles[iT];
1592 const V2d<T>& v1 = vertices[t.vertices[0]];
1593 const V2d<T>& v2 = vertices[t.vertices[1]];
1594 const V2d<T>& v3 = vertices[t.vertices[2]];
1595 const PtTriLocation::Enum loc = locatePointTriangle(v, v1, v2, v3);
1596
1597 if(loc == PtTriLocation::Outside)
1598 {
1599 handleException(
1600 Error("No triangle was found at position", CDT_SOURCE_LOCATION));
1601 }
1602 if(loc == PtTriLocation::OnVertex)
1603 {
1604 const VertInd iDupe = v1 == v ? t.vertices[0]
1605 : v2 == v ? t.vertices[1]
1606 : t.vertices[2];
1607 handleException(DuplicateVertexError(
1608 VertInd(iV - m_nTargetVerts),
1609 VertInd(iDupe - m_nTargetVerts),
1610 CDT_SOURCE_LOCATION));
1611 }
1612
1613 out[0] = iT;
1614 if(isOnEdge(loc))
1615 out[1] = t.neighbors[edgeNeighbor(loc)];
1616 return out;
1617}
1618
1619template <typename T, typename TNearPointLocator>
1620void Triangulation<T, TNearPointLocator>::changeNeighbor(
1621 const TriInd iT,
1622 const TriInd oldNeighbor,
1623 const TriInd newNeighbor)
1624{
1625 if(iT == noNeighbor)
1626 return;
1627 NeighborsArr3& nn = triangles[iT].neighbors;
1628 assert(
1629 nn[0] == oldNeighbor || nn[1] == oldNeighbor || nn[2] == oldNeighbor);
1630 if(nn[0] == oldNeighbor)
1631 nn[0] = newNeighbor;
1632 else if(nn[1] == oldNeighbor)
1633 nn[1] = newNeighbor;
1634 else
1635 nn[2] = newNeighbor;
1636}
1637
1638template <typename T, typename TNearPointLocator>
1639void Triangulation<T, TNearPointLocator>::changeNeighbor(
1640 const TriInd iT,
1641 const VertInd iVedge1,
1642 const VertInd iVedge2,
1643 const TriInd newNeighbor)
1644{
1645 assert(iT != noNeighbor);
1646 Triangle& t = triangles[iT];
1647 t.neighbors[edgeNeighborInd(t.vertices, iVedge1, iVedge2)] = newNeighbor;
1648}
1649
1650template <typename T, typename TNearPointLocator>
1651void Triangulation<T, TNearPointLocator>::triangulatePseudoPolygon(
1652 const std::vector<VertInd>& poly,
1653 unordered_map<Edge, TriInd>& outerTris,
1654 TriInd iT,
1655 TriInd iN,
1656 std::vector<TriInd>& trianglesToReuse,
1657 std::vector<TriangulatePseudoPolygonTask>& iterations)
1658{
1659 assert(poly.size() > 2);
1660 // note: uses iteration instead of recursion to avoid stack overflows
1661 iterations.clear();
1662 iterations.push_back(make_tuple(
1663 IndexSizeType(0),
1664 static_cast<IndexSizeType>(poly.size() - 1),
1665 iT,
1666 iN,
1667 Index(0)));
1668 while(!iterations.empty())
1669 {
1670 triangulatePseudoPolygonIteration(
1671 poly, outerTris, trianglesToReuse, iterations);
1672 }
1673}
1674
1675template <typename T, typename TNearPointLocator>
1676void Triangulation<T, TNearPointLocator>::triangulatePseudoPolygonIteration(
1677 const std::vector<VertInd>& poly,
1678 unordered_map<Edge, TriInd>& outerTris,
1679 std::vector<TriInd>& trianglesToReuse,
1680 std::vector<TriangulatePseudoPolygonTask>& iterations)
1681{
1682 IndexSizeType iA, iB;
1683 TriInd iT, iParent;
1684 Index iInParent;
1685 assert(!iterations.empty());
1686 tie(iA, iB, iT, iParent, iInParent) = iterations.back();
1687 iterations.pop_back();
1688 assert(iB - iA > 1 && iT != noNeighbor && iParent != noNeighbor);
1689 Triangle& t = triangles[iT];
1690 // find Delaunay point
1691 const IndexSizeType iC = findDelaunayPoint(poly, iA, iB);
1692
1693 const VertInd a = poly[iA];
1694 const VertInd b = poly[iB];
1695 const VertInd c = poly[iC];
1696
1697 // split pseudo-polygon in two parts and triangulate them
1698 // note: second part needs to be pushed on stack first to be processed first
1699
1700 // second part: points after the Delaunay point
1701 if(iB - iC > 1)
1702 {
1703 assert(!trianglesToReuse.empty());
1704 const TriInd iNext = trianglesToReuse.back();
1705 trianglesToReuse.pop_back();
1706 iterations.push_back(make_tuple(iC, iB, iNext, iT, Index(1)));
1707 }
1708 else // pseudo-poly is reduced to a single outer edge
1709 {
1710 const Edge outerEdge(b, c);
1711 const TriInd outerTri = outerTris.at(outerEdge);
1712 if(outerTri != noNeighbor)
1713 {
1714 assert(outerTri != iT);
1715 t.neighbors[1] = outerTri;
1716 changeNeighbor(outerTri, c, b, iT);
1717 }
1718 else
1719 outerTris.at(outerEdge) = iT;
1720 }
1721 // first part: points before the Delaunay point
1722 if(iC - iA > 1)
1723 { // add next triangle and add another iteration
1724 assert(!trianglesToReuse.empty());
1725 const TriInd iNext = trianglesToReuse.back();
1726 trianglesToReuse.pop_back();
1727 iterations.push_back(make_tuple(iA, iC, iNext, iT, Index(2)));
1728 }
1729 else
1730 { // pseudo-poly is reduced to a single outer edge
1731 const Edge outerEdge(c, a);
1732 const TriInd outerTri = outerTris.at(outerEdge);
1733 if(outerTri != noNeighbor)
1734 {
1735 assert(outerTri != iT);
1736 t.neighbors[2] = outerTri;
1737 changeNeighbor(outerTri, c, a, iT);
1738 }
1739 else
1740 outerTris.at(outerEdge) = iT;
1741 }
1742 // Finalize triangle
1743 // note: only when triangle is finalized to we add it as a neighbor to
1744 // parent to maintain triangulation topology consistency
1745 triangles[iParent].neighbors[iInParent] = iT;
1746 t.neighbors[0] = iParent;
1747 t.vertices = arr3(a, b, c);
1748 setAdjacentTriangle(c, iT);
1749}
1750
1751template <typename T, typename TNearPointLocator>
1752IndexSizeType Triangulation<T, TNearPointLocator>::findDelaunayPoint(
1753 const std::vector<VertInd>& poly,
1754 const IndexSizeType iA,
1755 const IndexSizeType iB) const
1756{
1757 assert(iB - iA > 1);
1758 const V2d<T>& a = vertices[poly[iA]];
1759 const V2d<T>& b = vertices[poly[iB]];
1760 IndexSizeType out = iA + 1;
1761 const V2d<T>* c = &vertices[poly[out]]; // caching for better performance
1762 for(IndexSizeType i = iA + 1; i < iB; ++i)
1763 {
1764 const V2d<T>& v = vertices[poly[i]];
1765 if(isInCircumcircle(v, a, b, *c))
1766 {
1767 out = i;
1768 c = &v;
1769 }
1770 }
1771 assert(out > iA && out < iB); // point is between ends
1772 return out;
1773}
1774
1775template <typename T, typename TNearPointLocator>
1777 const std::vector<V2d<T> >& newVertices)
1778{
1779 return insertVertices(
1780 newVertices.begin(), newVertices.end(), getX_V2d<T>, getY_V2d<T>);
1781}
1782
1783template <typename T, typename TNearPointLocator>
1785{
1786 return m_vertTris.empty() && !vertices.empty();
1787}
1788
1789template <typename T, typename TNearPointLocator>
1790unordered_map<TriInd, LayerDepth>
1791Triangulation<T, TNearPointLocator>::peelLayer(
1792 std::stack<TriInd> seeds,
1793 const LayerDepth layerDepth,
1794 std::vector<LayerDepth>& triDepths) const
1795{
1796 unordered_map<TriInd, LayerDepth> behindBoundary;
1797 while(!seeds.empty())
1798 {
1799 const TriInd iT = seeds.top();
1800 seeds.pop();
1801 triDepths[iT] = std::min(triDepths[iT], layerDepth);
1802 behindBoundary.erase(iT);
1803 const Triangle& t = triangles[iT];
1804 for(Index i(0); i < Index(3); ++i)
1805 {
1806 const Edge opEdge(t.vertices[ccw(i)], t.vertices[cw(i)]);
1807 const TriInd iN = t.neighbors[opoNbr(i)];
1808 if(iN == noNeighbor || triDepths[iN] <= layerDepth)
1809 continue;
1810 if(fixedEdges.count(opEdge))
1811 {
1812 const unordered_map<Edge, LayerDepth>::const_iterator cit =
1813 overlapCount.find(opEdge);
1814 const LayerDepth triDepth = cit == overlapCount.end()
1815 ? layerDepth + 1
1816 : layerDepth + cit->second + 1;
1817 behindBoundary[iN] = triDepth;
1818 continue;
1819 }
1820 seeds.push(iN);
1821 }
1822 }
1823 return behindBoundary;
1824}
1825
1826template <typename T, typename TNearPointLocator>
1827std::vector<LayerDepth>
1829{
1830 std::vector<LayerDepth> triDepths(
1831 triangles.size(), std::numeric_limits<LayerDepth>::max());
1832 std::stack<TriInd> seeds(TriDeque(1, m_vertTris[0]));
1833 LayerDepth layerDepth = 0;
1834 LayerDepth deepestSeedDepth = 0;
1835
1836 unordered_map<LayerDepth, TriIndUSet> seedsByDepth;
1837 do
1838 {
1839 const unordered_map<TriInd, LayerDepth>& newSeeds =
1840 peelLayer(seeds, layerDepth, triDepths);
1841
1842 seedsByDepth.erase(layerDepth);
1843 typedef unordered_map<TriInd, LayerDepth>::const_iterator Iter;
1844 for(Iter it = newSeeds.begin(); it != newSeeds.end(); ++it)
1845 {
1846 deepestSeedDepth = std::max(deepestSeedDepth, it->second);
1847 seedsByDepth[it->second].insert(it->first);
1848 }
1849 const TriIndUSet& nextLayerSeeds = seedsByDepth[layerDepth + 1];
1850 seeds = std::stack<TriInd>(
1851 TriDeque(nextLayerSeeds.begin(), nextLayerSeeds.end()));
1852 ++layerDepth;
1853 } while(!seeds.empty() || deepestSeedDepth > layerDepth);
1854
1855 return triDepths;
1856}
1857
1858#ifdef CDT_ENABLE_CALLBACK_HANDLER
1859template <typename T, typename TNearPointLocator>
1861 ICallbackHandler* callbackHandler)
1862{
1863 m_callbackHandler = callbackHandler;
1864}
1865#endif
1866
1867template <typename T, typename TNearPointLocator>
1868void Triangulation<T, TNearPointLocator>::insertVertices_AsProvided(
1869 VertInd superGeomVertCount)
1870{
1871 for(VertInd iV = superGeomVertCount; iV < vertices.size(); ++iV)
1872 {
1873#ifdef CDT_ENABLE_CALLBACK_HANDLER
1874 if(m_callbackHandler && m_callbackHandler->isAbortCalculation())
1875 {
1876 return;
1877 }
1878#endif
1879 insertVertex(iV);
1880 }
1881}
1882
1883template <typename T, typename TNearPointLocator>
1884void Triangulation<T, TNearPointLocator>::insertVertices_Randomized(
1885 VertInd superGeomVertCount)
1886{
1887 std::size_t vertexCount = vertices.size() - superGeomVertCount;
1888 std::vector<VertInd> ii(vertexCount);
1889 detail::iota(ii.begin(), ii.end(), superGeomVertCount);
1890 detail::random_shuffle(ii.begin(), ii.end());
1891 for(std::vector<VertInd>::iterator it = ii.begin(); it != ii.end(); ++it)
1892 {
1893#ifdef CDT_ENABLE_CALLBACK_HANDLER
1894 if(m_callbackHandler && m_callbackHandler->isAbortCalculation())
1895 {
1896 return;
1897 }
1898#endif
1899 insertVertex(*it);
1900 }
1901}
1902
1903namespace detail
1904{
1905
1906// log2 implementation backwards compatible with pre c++11
1907template <typename T>
1908inline double log2_bc(T x)
1909{
1910#ifdef CDT_CXX11_IS_SUPPORTED
1911 return std::log2(x);
1912#else
1913 static double log2_constant = std::log(2.0);
1914 return std::log(static_cast<double>(x)) / log2_constant;
1915#endif
1916}
1917
1922inline std::size_t maxQueueLengthBFSKDTree(const std::size_t vertexCount)
1923{
1924 const int filledLayerPow2 =
1925 static_cast<int>(std::floor(log2_bc(vertexCount)) - 1);
1926 const std::size_t nodesInFilledTree =
1927 static_cast<std::size_t>(std::pow(2., filledLayerPow2 + 1) - 1);
1928 const std::size_t nodesInLastFilledLayer =
1929 static_cast<std::size_t>(std::pow(2., filledLayerPow2));
1930 const std::size_t nodesInLastLayer = vertexCount - nodesInFilledTree;
1931 return nodesInLastLayer >= nodesInLastFilledLayer
1932 ? nodesInLastFilledLayer + nodesInLastLayer -
1933 nodesInLastFilledLayer
1934 : nodesInLastFilledLayer;
1935}
1936
1937template <typename T>
1938class FixedCapacityQueue
1939{
1940public:
1941 FixedCapacityQueue(const std::size_t capacity)
1942 : m_vec(capacity)
1943 , m_front(m_vec.begin())
1944 , m_back(m_vec.begin())
1945 , m_size(0)
1946 {}
1947 bool empty() const
1948 {
1949 return m_size == 0;
1950 }
1951 const T& front() const
1952 {
1953 return *m_front;
1954 }
1955 void pop()
1956 {
1957 assert(m_size > 0);
1958 ++m_front;
1959 if(m_front == m_vec.end())
1960 m_front = m_vec.begin();
1961 --m_size;
1962 }
1963 void push(const T& t)
1964 {
1965 assert(m_size < m_vec.size());
1966 *m_back = t;
1967 ++m_back;
1968 if(m_back == m_vec.end())
1969 m_back = m_vec.begin();
1970 ++m_size;
1971 }
1972#ifdef CDT_CXX11_IS_SUPPORTED
1973 void push(const T&& t)
1974 {
1975 assert(m_size < m_vec.size());
1976 *m_back = t;
1977 ++m_back;
1978 if(m_back == m_vec.end())
1979 m_back = m_vec.begin();
1980 ++m_size;
1981 }
1982#endif
1983private:
1984 std::vector<T> m_vec;
1985 typename std::vector<T>::iterator m_front;
1986 typename std::vector<T>::iterator m_back;
1987 std::size_t m_size;
1988};
1989
1990template <typename T>
1991class less_than_x
1992{
1993 const std::vector<V2d<T> >& m_vertices;
1994
1995public:
1996 less_than_x(const std::vector<V2d<T> >& vertices)
1997 : m_vertices(vertices)
1998 {}
1999 bool operator()(const VertInd a, const VertInd b) const
2000 {
2001 return m_vertices[a].x < m_vertices[b].x;
2002 }
2003};
2004
2005template <typename T>
2006class less_than_y
2007{
2008 const std::vector<V2d<T> >& m_vertices;
2009
2010public:
2011 less_than_y(const std::vector<V2d<T> >& vertices)
2012 : m_vertices(vertices)
2013 {}
2014 bool operator()(const VertInd a, const VertInd b) const
2015 {
2016 return m_vertices[a].y < m_vertices[b].y;
2017 }
2018};
2019
2020} // namespace detail
2021
2022template <typename T, typename TNearPointLocator>
2023void Triangulation<T, TNearPointLocator>::insertVertices_KDTreeBFS(
2024 VertInd superGeomVertCount,
2025 Box2d<T> box)
2026{
2027 // calculate original indices
2028 const VertInd vertexCount(vertices.size() - superGeomVertCount);
2029 if(vertexCount <= VertInd(0))
2030 return;
2031 std::vector<VertInd> ii(vertexCount);
2032 detail::iota(ii.begin(), ii.end(), superGeomVertCount);
2033
2034 typedef std::vector<VertInd>::iterator It;
2036 detail::maxQueueLengthBFSKDTree(vertexCount));
2037 queue.push(make_tuple(ii.begin(), ii.end(), box.min, box.max, VertInd(0)));
2038
2039 It first, last;
2040 V2d<T> newBoxMin, newBoxMax;
2041 VertInd parent, mid;
2042
2043 const detail::less_than_x<T> cmpX(vertices);
2044 const detail::less_than_y<T> cmpY(vertices);
2045
2046 while(!queue.empty())
2047 {
2048#ifdef CDT_ENABLE_CALLBACK_HANDLER
2049 if(m_callbackHandler && m_callbackHandler->isAbortCalculation())
2050 {
2051 return;
2052 }
2053#endif
2054 tie(first, last, box.min, box.max, parent) = queue.front();
2055 queue.pop();
2056 assert(first != last);
2057
2058 const std::ptrdiff_t len = std::distance(first, last);
2059 if(len == 1)
2060 {
2061 insertVertex(*first, parent);
2062 continue;
2063 }
2064 const It midIt = first + len / 2;
2065 if(box.max.x - box.min.x >= box.max.y - box.min.y)
2066 {
2067 detail::portable_nth_element(first, midIt, last, cmpX);
2068 mid = *midIt;
2069 const T split = vertices[mid].x;
2070 newBoxMin.x = split;
2071 newBoxMin.y = box.min.y;
2072 newBoxMax.x = split;
2073 newBoxMax.y = box.max.y;
2074 }
2075 else
2076 {
2077 detail::portable_nth_element(first, midIt, last, cmpY);
2078 mid = *midIt;
2079 const T split = vertices[mid].y;
2080 newBoxMin.x = box.min.x;
2081 newBoxMin.y = split;
2082 newBoxMax.x = box.max.x;
2083 newBoxMax.y = split;
2084 }
2085 insertVertex(mid, parent);
2086 if(first != midIt)
2087 {
2088 queue.push(make_tuple(first, midIt, box.min, newBoxMax, mid));
2089 }
2090 if(midIt + 1 != last)
2091 {
2092 queue.push(make_tuple(midIt + 1, last, newBoxMin, box.max, mid));
2093 }
2094 }
2095}
2096
2097template <typename T, typename TNearPointLocator>
2098std::pair<TriInd, TriInd> Triangulation<T, TNearPointLocator>::edgeTriangles(
2099 const VertInd a,
2100 const VertInd b) const
2101{
2102 const TriInd triStart = m_vertTris[a];
2103 assert(triStart != noNeighbor);
2104 TriInd iT = triStart, iTNext = triStart;
2105 VertInd iV = noVertex;
2106 do
2107 {
2108 const Triangle& t = triangles[iT];
2109 tie(iTNext, iV) = t.next(a);
2110 assert(iTNext != noNeighbor);
2111 if(iV == b)
2112 {
2113 return std::make_pair(iT, iTNext);
2114 }
2115 iT = iTNext;
2116 } while(iT != triStart);
2117 return std::make_pair(noNeighbor, noNeighbor);
2118}
2119
2120template <typename T, typename TNearPointLocator>
2121bool Triangulation<T, TNearPointLocator>::hasEdge(
2122 const VertInd a,
2123 const VertInd b) const
2124{
2125 return edgeTriangles(a, b).first != invalidIndexSizeType;
2126}
2127
2128template <typename T, typename TNearPointLocator>
2129void Triangulation<T, TNearPointLocator>::setAdjacentTriangle(
2130 const VertInd v,
2131 const TriInd t)
2132{
2133 assert(t != noNeighbor);
2134 m_vertTris[v] = t;
2135 assert(
2136 triangles[t].vertices[0] == v || triangles[t].vertices[1] == v ||
2137 triangles[t].vertices[2] == v);
2138}
2139
2140template <typename T, typename TNearPointLocator>
2141void Triangulation<T, TNearPointLocator>::pivotVertexTriangleCW(const VertInd v)
2142{
2143 assert(m_vertTris[v] != noNeighbor);
2144 m_vertTris[v] = triangles[m_vertTris[v]].next(v).first;
2145 assert(m_vertTris[v] != noNeighbor);
2146 assert(
2147 triangles[m_vertTris[v]].vertices[0] == v ||
2148 triangles[m_vertTris[v]].vertices[1] == v ||
2149 triangles[m_vertTris[v]].vertices[2] == v);
2150}
2151
2152template <typename T, typename TNearPointLocator>
2153void Triangulation<T, TNearPointLocator>::tryAddVertexToLocator(const VertInd v)
2154{
2155 if(!m_nearPtLocator.empty()) // only if locator is initialized already
2156 m_nearPtLocator.addPoint(v, vertices);
2157}
2158
2159template <typename T, typename TNearPointLocator>
2160void Triangulation<T, TNearPointLocator>::tryInitNearestPointLocator()
2161{
2162 if(!vertices.empty() && m_nearPtLocator.empty())
2163 {
2164 m_nearPtLocator.initialize(vertices);
2165 }
2166}
2167
2168} // namespace CDT
Triangulation class.
void iota(ForwardIt first, ForwardIt last, T value)
backport from c++11
std::size_t maxQueueLengthBFSKDTree(const std::size_t vertexCount)
Since KD-tree bulk load builds a balanced tree the maximum length of a queue can be pre-calculated: i...
Base class for errors.
Interface for the callback handler that user can derive from and inject into the triangulation to mon...
Error thrown when intersecting constraint edges are detected, but triangulation is not configured to ...
void eraseOuterTriangles()
Erase triangles outside of constrained boundary using growing.
void conformToEdges(TEdgeIter first, TEdgeIter last, TGetEdgeVertexStart getStart, TGetEdgeVertexEnd getEnd)
Insert constraint edges into triangulation for Conforming Delaunay Triangulation (for example see fig...
std::vector< LayerDepth > calculateTriangleDepths() const
Calculate depth of each triangle in constraint triangulation.
bool isFinalized() const
Check if the triangulation was finalized with erase... method and super-triangle was removed.
void eraseOuterTrianglesAndHoles()
Erase triangles outside of constrained boundary and auto-detected holes.
V2dVec vertices
triangulation's vertices
void insertEdges(TEdgeIter first, TEdgeIter last, TGetEdgeVertexStart getStart, TGetEdgeVertexEnd getEnd)
Insert constraint edges into triangulation for Constrained Delaunay Triangulation (for example see fi...
void insertVertices(TVertexIter first, TVertexIter last, TGetVertexCoordX getX, TGetVertexCoordY getY)
Insert custom point-types specified by iterator range and X/Y-getters.
void eraseSuperTriangle()
Erase triangles adjacent to super triangle.
TriangleVec triangles
triangulation's triangles
Triangulation()
Default constructor.
void setCallbackHandler(ICallbackHandler *callbackHandler)
Set user-provided callback handler.
void initializedWithCustomSuperGeometry()
Call this method after directly setting custom super-geometry via vertices and triangles members.
unsigned short LayerDepth
Type used for storing layer depths for triangles.
Definition CDT.h:40
void flipEdge(TriInd iT, TriInd iTopo)
Flip an edge between two triangle.
TriIndVec & VertTrisInternal()
Access internal vertex adjacent triangles.
void removeTriangles(const TriIndUSet &removedTriangles)
Remove triangles with specified indices.
Namespace containing triangulation functionality.
std::vector< Edge > EdgeVec
Vector of edges.
Definition CDTUtils.h:330
CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY VertInd opposedVertex(const Triangle &tri, TriInd iTopo)
Given two triangles, return vertex of first triangle opposed to the second.
Definition CDTUtils.hpp:201
unordered_set< Edge > EdgeUSet
Hash table of edges.
Definition CDTUtils.h:331
VertInd edge_get_v2(const Edge &e)
Get edge second vertex.
Definition CDTUtils.h:319
std::vector< TriInd > TriIndVec
Vector of triangle indices.
Definition CDTUtils.h:209
CDT_INLINE_IF_HEADER_ONLY Index edgeNeighborInd(const VerticesArr3 &vv, VertInd iVedge1, VertInd iVedge2)
Index of triangle's neighbor opposed to an edge.
Definition CDTUtils.hpp:139
VertInd edge_get_v1(const Edge &e)
Get edge first vertex.
Definition CDTUtils.h:313
CDT_EXPORT PtLineLocation::Enum classifyOrientation(T orientation, T orientationTolerance=T(0))
Classify value of orient2d predicate.
Definition CDTUtils.hpp:60
array< TriInd, 3 > NeighborsArr3
array of three neighbors
Definition CDTUtils.h:211
CDT_EXPORT T distance(const V2d< T > &a, const V2d< T > &b)
Distance between two 2D points.
Definition CDTUtils.hpp:248
IndexSizeType VertInd
Vertex index.
Definition CDTUtils.h:191
CDT_EXPORT Index cw(Index i)
Advance vertex or neighbor index clockwise.
Definition CDTUtils.hpp:24
array< VertInd, 3 > VerticesArr3
array of three vertex indices
Definition CDTUtils.h:210
CDT_INLINE_IF_HEADER_ONLY bool touchesSuperTriangle(const Triangle &t)
Check if any of triangle's vertices belongs to a super-triangle.
Definition CDTUtils.hpp:259
array< T, 3 > arr3(const T &v0, const T &v1, const T &v2)
Needed for c++03 compatibility (no uniform initialization available)
Definition CDTUtils.h:119
CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index opoNbr(Index vertIndex)
Opposed neighbor index from vertex index.
Definition CDTUtils.hpp:102
CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index vertexInd(const VerticesArr3 &vv, VertInd iV)
If triangle has a given vertex return vertex-index.
Definition CDTUtils.hpp:184
CDT_EXPORT T orient2D(const V2d< T > &p, const V2d< T > &v1, const V2d< T > &v2)
Orient p against line v1-v2 2D: robust geometric predicate.
Definition CDTUtils.hpp:43
CDT_EXPORT PtLineLocation::Enum locatePointLine(const V2d< T > &p, const V2d< T > &v1, const V2d< T > &v2, T orientationTolerance=T(0))
Check if point lies to the left of, to the right of, or on a line.
Definition CDTUtils.hpp:49
unordered_set< TriInd > TriIndUSet
Hash table of triangles.
Definition CDTUtils.h:332
CDT_EXPORT Index ccw(Index i)
Advance vertex or neighbor index counter-clockwise.
Definition CDTUtils.hpp:19
CDT_EXPORT Index edgeNeighbor(PtTriLocation::Enum location)
Neighbor index from a on-edge location.
Definition CDTUtils.hpp:36
unordered_map< TriInd, TriInd > TriIndUMap
Triangle hash map.
Definition CDTUtils.h:333
const T & getX_V2d(const V2d< T > &v)
X- coordinate getter for V2d.
Definition CDTUtils.h:155
unsigned char Index
Index in triangle.
Definition CDTUtils.h:189
CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY TriInd opposedTriangle(const Triangle &tri, VertInd iVert)
Given triangle and a vertex find opposed triangle.
Definition CDTUtils.hpp:195
CDT_EXPORT bool isInCircumcircle(const V2d< T > &p, const V2d< T > &v1, const V2d< T > &v2, const V2d< T > &v3)
Test if point lies in a circumscribed circle of a triangle.
Definition CDTUtils.hpp:214
Edge RemapNoSuperTriangle(const Edge &e)
Remap removing super-triangle: subtract 3 from vertices.
CDT_EXPORT bool isOnEdge(PtTriLocation::Enum location)
Check if location is classified as on any of three edges.
Definition CDTUtils.hpp:29
IndexSizeType TriInd
Triangle index.
Definition CDTUtils.h:193
CDT_EXPORT PtTriLocation::Enum locatePointTriangle(const V2d< T > &p, const V2d< T > &v1, const V2d< T > &v2, const V2d< T > &v3)
Check if point a lies inside of, outside of, or on an edge of a triangle.
Definition CDTUtils.hpp:70
CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index opposedVertexInd(const NeighborsArr3 &nn, TriInd iTopo)
Index of triangle's vertex opposed to a triangle.
Definition CDTUtils.hpp:173
const T & getY_V2d(const V2d< T > &v)
Y-coordinate getter for V2d.
Definition CDTUtils.h:162
@ FixedEdgeMidpoint
During conforming triangulation edge mid-point is added.
@ FixedEdgesIntersection
Resolving fixed/constraint edges' intersection.
2D bounding box
Definition CDTUtils.h:216
V2d< T > max
max box corner
Definition CDTUtils.h:218
V2d< T > min
min box corner
Definition CDTUtils.h:217
Edge connecting two vertices: vertex with smaller index is always first.
Definition CDTUtils.h:271
VertInd v1() const
V1 getter.
Definition CDTUtils.h:291
VertInd v2() const
V2 getter.
Definition CDTUtils.h:297
@ TryResolve
attempt to resolve constraint edge intersections
@ NotAllowed
constraint edge intersections are not allowed
@ DontCheck
No checks: slightly faster but less safe.
Enum
The Enum itself.
@ SuperTriangle
conventional super-triangle
@ Custom
user-specified custom geometry (e.g., grid)
Triangulation triangle (counter-clockwise winding)
Definition CDTUtils.h:344
VerticesArr3 vertices
triangle's three vertices
Definition CDTUtils.h:345
NeighborsArr3 neighbors
triangle's three neighbors
Definition CDTUtils.h:346
std::pair< TriInd, VertInd > next(const VertInd i) const
Next triangle adjacent to a vertex (clockwise)
Definition CDTUtils.h:362
2D vector
Definition CDTUtils.h:136
T y
Y-coordinate.
Definition CDTUtils.h:138
T x
X-coordinate.
Definition CDTUtils.h:137
@ Auto
Automatic insertion order optimized for better performance.