juli/OrcaSlicer
1
1
Fork 0
mirror of https://github.com/OrcaSlicer/OrcaSlicer.git synced 2026-10-04 14:20:58 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
5704fa348f9fd598f28b660ce608780e29e45de1
OrcaSlicer/deps_src/clipper2/Clipper2Lib/src/clipper.triangulation.cpp
T

1223 lines
37 KiB
C++
Raw Blame History

/*******************************************************************************
* Author : Angus Johnson *
* Date : 13 December 2025 *
* Release : BETA RELEASE *
* Website : https://www.angusj.com *
* Copyright : Angus Johnson 2010-2025 *
* Purpose : Constrained Delaunay Triangulation *
* License : https://www.boost.org/LICENSE_1_0.txt *
*******************************************************************************/
#include "clipper2/clipper.h"
#include "clipper2/clipper.triangulation.h"
#ifdef USINGZ
namespace Clipper2Lib_Z {
#else
namespace Clipper2Lib {
#endif
enum class EdgeKind { loose, ascend, descend }; // ascend & descend are 'fixed' edges
enum class IntersectKind { none, collinear, intersect };
enum class EdgeContainsResult { neither, left, right };
//forward definitions
class Vertex2;
class Edge;
class Triangle;
typedef std::vector<Vertex2*> VertexList;
typedef std::vector<Edge*> EdgeList;
class Vertex2 {
public:
Point64 pt;
EdgeList edges;
bool innerLM = false;
Vertex2(const Point64& p64) : pt(p64) { edges.reserve(2); };
};
class Edge {
public:
Vertex2* vL = nullptr;
Vertex2* vR = nullptr;
Vertex2* vB = nullptr;
Vertex2* vT = nullptr;
EdgeKind kind = EdgeKind::loose;
Triangle* triA = nullptr;
Triangle* triB = nullptr;
bool isActive = false;
Edge* nextE = nullptr;
Edge* prevE = nullptr;
};
class Triangle {
public:
Edge* edges[3];
Triangle(Edge* e1, Edge* e2, Edge* e3)
{
edges[0] = e1;
edges[1] = e2;
edges[2] = e3;
}
};
/////////////////////////////////////////////////////////////////////////////
// Delaunay class declaration
/////////////////////////////////////////////////////////////////////////////
class Delaunay {
private:
VertexList allVertices;
EdgeList allEdges;
std::vector<Triangle*> allTriangles;
std::stack<Edge*> pendingDelaunayStack;
std::stack<Edge*> horzEdgeStack;
std::stack<Vertex2*> locMinStack;
bool useDelaunay = true;
Vertex2* lowermostVertex = nullptr;
Edge* firstActive = nullptr;
void AddPath(const Path64& path);
bool AddPaths(const Paths64& paths);
void CleanUp();
bool FixupEdgeIntersects();
void MergeDupOrCollinearVertices();
void SplitEdge(Edge* longE, Edge* shortE);
bool RemoveIntersection(Edge* e1, Edge* e2);
Edge* CreateEdge(Vertex2* v1, Vertex2* v2, EdgeKind k);
Triangle* CreateTriangle(Edge* e1, Edge* e2, Edge* e3);
Edge* CreateInnerLocMinLooseEdge(Vertex2* vAbove);
Edge* HorizontalBetween(Vertex2* v1, Vertex2* v2);
void DoTriangulateLeft(Edge* edge, Vertex2* pivot, int64_t minY);
void DoTriangulateRight(Edge* edge, Vertex2* pivot, int64_t minY);
void AddEdgeToActives(Edge* edge);
void RemoveEdgeFromActives(Edge* edge);
void ForceLegal(Edge* edge);
public:
explicit Delaunay(bool delaunay = true) : useDelaunay(delaunay) {};
~Delaunay() { CleanUp(); };
Paths64 Execute(const Paths64& paths, TriangulateResult& triResult);
};
/////////////////////////////////////////////////////////////////////////////
// Miscellaneous functions
/////////////////////////////////////////////////////////////////////////////
static bool VertexListSort(const Vertex2* a, const Vertex2* b)
{
return (a->pt.y == b->pt.y) ? (a->pt.x < b->pt.x) : (a->pt.y > b->pt.y);
}
static bool EdgeListSort(const Edge* a, const Edge* b)
{
return (a->vL->pt.x < b->vL->pt.x);
}
static bool IsLooseEdge(const Edge& e)
{
return e.kind == EdgeKind::loose;
}
static bool IsLeftEdge(const Edge& e)
{
// left edge (bound) of a fill region
// ie fills on the **right** side of the edge
// precondition - e is never a 'loose' edge
return e.kind == EdgeKind::ascend;
}
static bool IsRightEdge(const Edge& e)
{
// right edge (bound) of a fill region
// but still fills on the **right** side of the edge
// precondition - e is never a 'loose' edge
return e.kind == EdgeKind::descend;
}
static bool IsHorizontal(const Edge& e)
{
return e.vB->pt.y == e.vT->pt.y;
}
static bool LeftTurning(const Point64& p1, const Point64& p2, const Point64& p3)
{
return CrossProductSign(p1, p2, p3) < 0;
}
static bool RightTurning(const Point64& p1, const Point64& p2, const Point64& p3)
{
return CrossProductSign(p1, p2, p3) > 0;
}
static bool EdgeCompleted(Edge* edge)
{
if (!edge->triA) return false;
if (edge->triB) return true;
return edge->kind != EdgeKind::loose;
}
static EdgeContainsResult EdgeContains(const Edge* edge, const Vertex2* v)
{
if (edge->vL == v) return EdgeContainsResult::left;
else if (edge->vR == v) return EdgeContainsResult::right;
else return EdgeContainsResult::neither;
}
static double GetAngle(const Point64& a, const Point64& b, const Point64& c)
{
//https://stackoverflow.com/a/3487062/359538
double abx = static_cast<double>(b.x - a.x);
double aby = static_cast<double>(b.y - a.y);
double bcx = static_cast<double>(b.x - c.x);
double bcy = static_cast<double>(b.y - c.y);
double dp = (abx * bcx + aby * bcy);
double cp = (abx * bcy - aby * bcx);
return std::atan2(cp, dp); //range between -Pi and Pi
}
static double GetLocMinAngle(Vertex2* v)
{
// todo - recheck the result's sign compared to left vs right turning
// (currently assumes left turning => positive values)
// precondition - this function is called before processing locMin.
//Assert(v->edges.size() == 2);
int asc, des;
if (v->edges[0]->kind == EdgeKind::ascend)
{
asc = 0;
des = 1;
}
else
{
des = 0;
asc = 1;
}
// winding direction - descending to ascending
return GetAngle(v->edges[des]->vT->pt, v->pt, v->edges[asc]->vT->pt);
}
static void RemoveEdgeFromVertex(Vertex2* vert, Edge* edge)
{
auto it = std::find(vert->edges.begin(), vert->edges.end(), edge);
if (it == vert->edges.end()) throw "oops!";
vert->edges.erase(it);
}
static bool FindLocMinIdx(const Path64& path, size_t len, size_t& idx)
{
if (len < 3) return false;
size_t i0 = idx, n = (idx + 1) % len;
while (path[n].y <= path[idx].y)
{
idx = n; n = (n + 1) % len;
if (idx == i0) return false; // fails if the path is completely horizontal
}
while (path[n].y >= path[idx].y)
{
idx = n; n = (n + 1) % len;
}
return true;
}
static size_t Prev(size_t& idx, size_t len)
{
if (idx == 0) return len - 1; else return idx - 1;
}
static size_t Next(size_t& idx, size_t len)
{
return (idx + 1) % len;
}
static Edge* FindLinkingEdge(const Vertex2* vert1, const Vertex2* vert2, bool preferAscending)
{
Edge* res = nullptr;
for (auto e : vert1->edges)
{
if (e->vL == vert2 || e->vR == vert2)
{
if (e->kind == EdgeKind::loose ||
((e->kind == EdgeKind::ascend) == preferAscending)) return e;
res = e;
}
}
return res;
}
static Path64 PathFromTriangle(Triangle tri)
{
Path64 res;
res.reserve(3);
res.push_back(tri.edges[0]->vL->pt);
res.push_back(tri.edges[0]->vR->pt);
const Edge& e = *tri.edges[1];
if (e.vL->pt == res[0] || e.vL->pt == res[1])
res.push_back(e.vR->pt);
else
res.push_back(e.vL->pt);
return res;
}
static double InCircleTest(const Point64& ptA, const Point64& ptB,
const Point64& ptC, const Point64& ptD)
{
// Return the determinant value of 3 x 3 matrix ...
// | ax-dx ay-dy Sqr(ax-dx)+Sqr(ay-dy) |
// | bx-dx by-dy Sqr(bx-dx)+Sqr(by-dy) |
// | cx-dx cy-dy Sqr(cx-dx)+Sqr(cy-dy) |
// The *sign* of the return value is determined by
// the orientation (CW vs CCW) of ptA, ptB & ptC.
// precondition - ptA, ptB & ptC make a *non-empty* triangle
double m00 = static_cast<double>(ptA.x - ptD.x);
double m01 = static_cast<double>(ptA.y - ptD.y);
double m02 = (Sqr(m00) + Sqr(m01));
double m10 = static_cast<double>(ptB.x - ptD.x);
double m11 = static_cast<double>(ptB.y - ptD.y);
double m12 = (Sqr(m10) + Sqr(m11));
double m20 = static_cast<double>(ptC.x - ptD.x);
double m21 = static_cast<double>(ptC.y - ptD.y);
double m22 = (Sqr(m20) + Sqr(m21));
return m00 * (m11 * m22 - m21 * m12) -
m10 * (m01 * m22 - m21 * m02) +
m20 * (m01 * m12 - m11 * m02);
}
static double ShortestDistFromSegment(const Point64& pt, const Point64& segPt1, const Point64& segPt2)
{
// precondition: segPt1 <> segPt2
double dx = static_cast<double>(segPt2.x - segPt1.x);
double dy = static_cast<double>(segPt2.y - segPt1.y);
//Assert((dx < > 0) or (dy < > 0)); // ie segPt1 <> segPt2
double ax = static_cast<double>(pt.x - segPt1.x);
double ay = static_cast<double>(pt.y - segPt1.y);
//q = (ax * dx + ay * dy) / (dx * dx + dy * dy);
double qNum = ax * dx + ay * dy;
if (qNum < 0) // pt is closest to seq1
return DistanceSqr(pt, segPt1);
else if (qNum > (Sqr(dx) + Sqr(dy))) // pt is closest to seq2
return DistanceSqr(pt, segPt2);
else
return Sqr(ax * dy - dx * ay) / (dx * dx + dy * dy);
}
static IntersectKind SegsIntersect(const Point64 s1a, const Point64 s1b,
const Point64 s2a, const Point64 s2b)
{
double dy1 = static_cast<double>(s1b.y - s1a.y);
double dx1 = static_cast<double>(s1b.x - s1a.x);
double dy2 = static_cast<double>(s2b.y - s2a.y);
double dx2 = static_cast<double>(s2b.x - s2a.x);
double cp = dy1 * dx2 - dy2 * dx1;
if (cp == 0) return IntersectKind::collinear;
double t = (static_cast<double>(s1a.x - s2a.x) * dy2 -
static_cast<double>(s1a.y - s2a.y) * dx2);
//ignore segments that 'intersect' at an end-point
if (t == 0) return IntersectKind::none;
if (t > 0)
{
if (cp < 0 || t >= cp) return IntersectKind::none;
}
else
{
if (cp > 0 || t <= cp) return IntersectKind::none;
}
// so far, the *segment* 's1' intersects the *line* through 's2',
// but now make sure it also intersects the *segment* 's2'
t = ((s1a.x - s2a.x) * dy1 - (s1a.y - s2a.y) * dx1);
if (t == 0) return IntersectKind::none;
if (t > 0)
{
if (cp > 0 && t < cp) return IntersectKind::intersect;
}
else
{
if (cp < 0 && t > cp) return IntersectKind::intersect;
}
return IntersectKind::none;
}
/////////////////////////////////////////////////////////////////////////////
// Delaunay class definitions
/////////////////////////////////////////////////////////////////////////////
void Delaunay::CleanUp()
{
for (auto v : allVertices) delete v;
allVertices.resize(0);
for (auto e : allEdges) delete e;
allEdges.resize(0);
for (auto t : allTriangles) delete t;
allTriangles.resize(0);
firstActive = nullptr;
lowermostVertex = nullptr;
}
void Delaunay::ForceLegal(Edge* edge)
{
// don't try to make empty triangles legal
if (!edge->triA || !edge->triB) return;
// vertA will be assigned the vertex in edge's triangleA
// that is NOT an end vertex of edge
// Likewise, vertB will be assigned the vertex in edge's
// triangleB that is NOT an end vertex of edge
// If edge is rotated, vertA & vertB will become its end vertices.
Vertex2* vertA = nullptr;
Vertex2* vertB = nullptr;
// Excluding 'edge', edgesA will contain two edges (one from
// triangleA and one from triangleB) that touch edge.vL.
// And edgesB will contain the two edges that touch edge.vR.
Edge* edgesA[3] = { nullptr, nullptr, nullptr };
Edge* edgesB[3] = { nullptr, nullptr, nullptr };
for (int i = 0; i < 3; ++i)
{
if (edge->triA->edges[i] == edge) continue;
switch (EdgeContains(edge->triA->edges[i], edge->vL))
{
case EdgeContainsResult::left:
edgesA[1] = edge->triA->edges[i];
vertA = edge->triA->edges[i]->vR;
break;
case EdgeContainsResult::right:
edgesA[1] = edge->triA->edges[i];
vertA = edge->triA->edges[i]->vL;
break;
default:
edgesB[1] = edge->triA->edges[i];
}
}
for (int i = 0; i < 3; ++i)
{
if (edge->triB->edges[i] == edge) continue;
switch (EdgeContains(edge->triB->edges[i], edge->vL))
{
case EdgeContainsResult::left:
edgesA[2] = edge->triB->edges[i];
vertB = edge->triB->edges[i]->vR;
break;
case EdgeContainsResult::right:
edgesA[2] = edge->triB->edges[i];
vertB = edge->triB->edges[i]->vL;
break;
default:
edgesB[2] = edge->triB->edges[i];
}
}
// InCircleTest reqires edge.triangleA to be a valid triangle
// if IsEmptyTriangle(edge.triangleA) then Exit; // slower
if (CrossProductSign(vertA->pt, edge->vL->pt, edge->vR->pt) == 0) return;
// ictResult - result sign is dependant on triangleA's orientation
double ictResult = InCircleTest(vertA->pt, edge->vL->pt, edge->vR->pt, vertB->pt);
if (ictResult == 0 || // if on or out of circle then exit
(RightTurning(vertA->pt, edge->vL->pt, edge->vR->pt) == (ictResult < 0))) return;
// TRIANGLES HERE ARE **NOT** DELAUNAY COMPLIANT, SO MAKE THEM SO.
// NOTE: ONCE WE BEGIN DELAUNAY COMPLIANCE, vL & vR WILL
// NO LONGER REPRESENT LEFT AND RIGHT VERTEX ORIENTATION.
// THIS IS MINOR PERFORMANCE EFFICIENCY IS SAFE AS LONG AS
// THE TRIANGULATE() METHOD IS CALLED ONCE ONLY ON A GIVEN
// SET OF PATHS
edge->vL = vertA;
edge->vR = vertB;
edge->triA->edges[0] = edge;
for (int i = 1; i < 3; ++i)
{
edge->triA->edges[i] = edgesA[i];
if (!edgesA[i]) throw "oops"; // stops compiler warnings
if (IsLooseEdge(*edgesA[i]))
pendingDelaunayStack.push(edgesA[i]);
// since each edge has its own triangleA and triangleB, we have to be careful
// to update the correct one ...
if (edgesA[i]->triA == edge->triA || edgesA[i]->triB == edge->triA) continue;
if (edgesA[i]->triA == edge->triB)
edgesA[i]->triA = edge->triA;
else if (edgesA[i]->triB == edge->triB)
edgesA[i]->triB = edge->triA;
else throw "oops";
}
edge->triB->edges[0] = edge;
for (int i = 1; i < 3; ++i)
{
edge->triB->edges[i] = edgesB[i];
if (!edgesB[i]) throw "oops"; // stops compiler warnings
if (IsLooseEdge(*edgesB[i]))
pendingDelaunayStack.push(edgesB[i]);
// since each edge has its own triangleA and triangleB, we have to be careful
// to update the correct one ...
if (edgesB[i]->triA == edge->triB || edgesB[i]->triB == edge->triB) continue;
if (edgesB[i]->triA == edge->triA)
edgesB[i]->triA = edge->triB;
else if (edgesB[i]->triB == edge->triA)
edgesB[i]->triB = edge->triB;
else throw "oops";
}
}
Edge* Delaunay::CreateEdge(Vertex2* v1, Vertex2* v2, EdgeKind k)
{
Edge* res = allEdges.emplace_back(new Edge());
if (v1->pt.y == v2->pt.y)
{
res->vB = v1; res->vT = v2;
}
else if (v1->pt.y < v2->pt.y)
{
res->vB = v2; res->vT = v1;
}
else
{
res->vB = v1; res->vT = v2;
}
if (v1->pt.x <= v2->pt.x)
{
res->vL = v1; res->vR = v2;
}
else
{
res->vL = v2; res->vR = v1;
}
res->kind = k;
v1->edges.push_back(res);
v2->edges.push_back(res);
if (k == EdgeKind::loose)
{
pendingDelaunayStack.push(res);
AddEdgeToActives(res);
}
return res;
}
Triangle* Delaunay::CreateTriangle(Edge* e1, Edge* e2, Edge* e3)
{
Triangle* res = allTriangles.emplace_back(new Triangle(e1, e2, e3));
// nb: only expire loose edges when both sides of these edges have triangles.
for (int i = 0; i < 3; ++i)
if (res->edges[i]->triA)
{
res->edges[i]->triB = res;
// this is the edge's second triangle hence no longer active
RemoveEdgeFromActives(res->edges[i]);
}
else
{
res->edges[i]->triA = res;
// this is the edge's first triangle, so only remove
// this edge from actives if it's a fixed edge.
if (!IsLooseEdge(*res->edges[i]))
RemoveEdgeFromActives(res->edges[i]);
}
return res;
}
bool Delaunay::RemoveIntersection(Edge* e1, Edge* e2)
{
// find which vertex is closest to the other segment
// (ie not the vertex closest to the intersection point) ...
Vertex2* v = e1->vL;
Edge* tmpE = e2;
double d = ShortestDistFromSegment(e1->vL->pt, e2->vL->pt, e2->vR->pt);
double d2 = ShortestDistFromSegment(e1->vR->pt, e2->vL->pt, e2->vR->pt);
if (d2 < d) { d = d2; v = e1->vR; }
d2 = ShortestDistFromSegment(e2->vL->pt, e1->vL->pt, e1->vR->pt);
if (d2 < d) { d = d2; tmpE = e1; v = e2->vL; }
d2 = ShortestDistFromSegment(e2->vR->pt, e1->vL->pt, e1->vR->pt);
if (d2 < d) { d = d2; tmpE = e1; v = e2->vR; }
if (d > 1.000)
return false; // Oops - this is not just a simple 'rounding' intersection
// split 'tmpE' into 2 edges at 'v'
Vertex2* v2 = tmpE->vT;
RemoveEdgeFromVertex(v2, tmpE);
// replace v2 in tmpE with v
if (tmpE->vL == v2) tmpE->vL = v;
else tmpE->vR = v;
tmpE->vT = v;
v->edges.push_back(tmpE);
v->innerLM = false; // #47
// left turning is angle positive
if (tmpE->vB->innerLM && GetLocMinAngle(tmpE->vB) <= 0)
tmpE->vB->innerLM = false; // #44, 52
// finally create a new edge between v and v2 ...
CreateEdge(v, v2, tmpE->kind);
return true;
}
bool Delaunay::FixupEdgeIntersects()
{
// precondition - edgeList must be sorted - ascending on edge.vL.pt.x
for (size_t i1 = 0; i1 < allEdges.size(); ++i1)
{
Edge* e1 = allEdges[i1];
// nb: we can safely ignore edges newly created inside this for loop
for (size_t i2 = i1 + 1; i2 < allEdges.size(); ++i2)
{
Edge* e2 = allEdges[i2];
if (e2->vL->pt.x >= e1->vR->pt.x)
break; // all 'e' from now on are too far right
// 'e2' is inside e1's horizontal region. If 'e2' is also inside
// e1's vertical region, only then check for an intersection ...
if (e2->vT->pt.y < e1->vB->pt.y && e2->vB->pt.y > e1->vT->pt.y &&
(SegsIntersect(e2->vL->pt, e2->vR->pt,
e1->vL->pt, e1->vR->pt) == IntersectKind::intersect))
{
if (!RemoveIntersection(e2, e1))
return false; // oops!!
}
// nb: collinear edges are managed in MergeDupOrCollinearVertices below
}
}
return true;
}
void Delaunay::SplitEdge(Edge* longE, Edge* shortE)
{
auto oldT = longE->vT, newT = shortE->vT;
// remove longEdge from longEdge.vT.edges
RemoveEdgeFromVertex(oldT, longE);
// shorten longEdge
longE->vT = newT;
if (longE->vL == oldT) longE->vL = newT;
else longE->vR = newT;
// add shortened longEdge to newT.edges
newT->edges.push_back(longE);
// and create a new edge betweem newV, oldT
CreateEdge(newT, oldT, longE->kind);
}
void Delaunay::MergeDupOrCollinearVertices()
{
// note: this procedure may add new edges and change the
// number of edges connected with a given vertex, but it
// won't add or delete vertices (so it's safe to use iterators)
auto vIter1 = allVertices.begin();
for (auto vIter2 = allVertices.begin() + 1; vIter2 != allVertices.end(); ++vIter2)
{
if ((*vIter1)->pt != (*vIter2)->pt)
{
vIter1 = vIter2;
continue;
}
// merge v1 & v2
Vertex2* v1 = *vIter1, * v2 = *vIter2;
if (!v1->innerLM || !v2->innerLM)
v1->innerLM = false;
// in all of v2's edges, replace links to v2 with links to v1
for (auto e : v2->edges)
{
if (e->vB == v2) e->vB = v1;
else e->vT = v1;
if (e->vL == v2) e->vL = v1;
else e->vR = v1;
}
// move all of v2's edges to v1
std::copy(v2->edges.begin(), v2->edges.end(), back_inserter(v1->edges));
v2->edges.resize(0);
// excluding horizontals, if pv.edges contains two edges
// that are *collinear* and share the same bottom coords
// but have different lengths, split the longer edge at
// the top of the shorter edge ...
for (auto itE = v1->edges.begin(); itE != v1->edges.end(); ++itE)
{
if (IsHorizontal(*(*itE)) || (*itE)->vB != v1) continue;
for (auto itE2 = itE + 1; itE2 != v1->edges.end(); ++itE2)
{
auto e1 = *itE, e2 = *itE2;
if (e2->vB != v1 || e1->vT->pt.y == e2->vT->pt.y ||
(CrossProductSign(e1->vT->pt, v1->pt, e2->vT->pt) != 0)) continue;
// we have parallel edges, both heading up from v1.pt.
// split the longer edge at the top of the shorter edge.
if (e1->vT->pt.y < e2->vT->pt.y) SplitEdge(e1, e2);
else SplitEdge(e2, e1);
break; // because only two edges can be collinear
}
}
}
}
Edge* Delaunay::CreateInnerLocMinLooseEdge(Vertex2* vAbove)
{
if (!firstActive) return nullptr; // oops!!
int64_t xAbove = vAbove->pt.x;
int64_t yAbove = vAbove->pt.y;
// find the closest 'active' edge with a vertex that's not above vAbove
Edge* e = firstActive, *eBelow = nullptr;
double bestD = -1.0;
while (e)
{
if (e->vL->pt.x <= xAbove && e->vR->pt.x >= xAbove &&
e->vB->pt.y >= yAbove && e->vB != vAbove && e->vT != vAbove &&
!LeftTurning(e->vL->pt, vAbove->pt, e->vR->pt))
{
double d = ShortestDistFromSegment(vAbove->pt, e->vL->pt, e->vR->pt);
if (!eBelow || d < bestD) // compare e with eBelow
{
eBelow = e;
bestD = d;
}
}
e = e->nextE;
}
if (!eBelow) return nullptr; // oops!!
// get the best vertex from 'eBelow'
Vertex2* vBest = (eBelow->vT->pt.y <= yAbove) ? eBelow->vB : eBelow->vT;
int64_t xBest = vBest->pt.x;
int64_t yBest = vBest->pt.y;
// make sure no edges intersect 'vAbove' and 'vBest' ...
// todo: fActives is currently *unsorted* but consider making it
// a tree structure based on each edge's left and right bounds
e = firstActive;
if (xBest < xAbove)
{
while (e)
{
if (e->vR->pt.x > xBest && e->vL->pt.x < xAbove &&
e->vB->pt.y > yAbove && e->vT->pt.y < yBest &&
(SegsIntersect(e->vB->pt, e->vT->pt,
vBest->pt, vAbove->pt) == IntersectKind::intersect))
{
vBest = (e->vT->pt.y > yAbove) ? e->vT : e->vB;
xBest = vBest->pt.x;
yBest = vBest->pt.y;
}
e = e->nextE;
}
}
else
{
while (e)
{
if (e->vR->pt.x < xBest && e->vL->pt.x > xAbove &&
e->vB->pt.y > yAbove && e->vT->pt.y < yBest &&
(SegsIntersect(e->vB->pt, e->vT->pt,
vBest->pt, vAbove->pt) == IntersectKind::intersect))
{
vBest = e->vT->pt.y > yAbove ? e->vT : e->vB;
xBest = vBest->pt.x;
yBest = vBest->pt.y;
}
e = e->nextE;
}
}
return CreateEdge(vBest, vAbove, EdgeKind::loose);
}
Edge* Delaunay::HorizontalBetween(Vertex2* v1, Vertex2* v2)
{
int64_t y = v1->pt.y, l, r;
if (v1->pt.x > v2->pt.x)
{
l = v2->pt.x;
r = v1->pt.x;
}
else
{
l = v1->pt.x;
r = v2->pt.x;
}
Edge* res = firstActive;
while (res)
{
if (res->vL->pt.y == y && res->vR->pt.y == y &&
res->vL->pt.x >= l && res->vR->pt.x <= r &&
(res->vL->pt.x != l || res->vL->pt.x != r)) break;
res = res->nextE;
}
return res;
}
void Delaunay::DoTriangulateLeft(Edge* edge, Vertex2* pivot, int64_t minY)
{
// precondition - pivot must be one end of edge (Usually .vB)
//Assert(!EdgeCompleted(edge));
Vertex2* vAlt = nullptr;
Edge* eAlt = nullptr;
Vertex2* v = (edge->vB == pivot) ? edge->vT : edge->vB;
for (auto e : pivot->edges)
{
if (e == edge || !e->isActive) continue;
Vertex2* vX = e->vT == pivot ? e->vB : e->vT;
if (vX == v) continue;
int cps = CrossProductSign(v->pt, pivot->pt, vX->pt);
if (cps == 0) //collinear paths
{
// if pivot is between v and vX then continue;
// nb: this is important for both horiz and non-horiz collinear
if ((v->pt.x > pivot->pt.x) == (pivot->pt.x > vX->pt.x)) continue;
}
// else if right-turning or not the best edge, then continue
else if (cps > 0 || (vAlt && !LeftTurning(vX->pt, pivot->pt, vAlt->pt)))
continue;
vAlt = vX;
eAlt = e;
}
if (!vAlt || vAlt->pt.y < minY) return;
// Don't triangulate **across** fixed edges
if (vAlt->pt.y < pivot->pt.y)
{
if (IsLeftEdge(*eAlt)) return;
}
else if (vAlt->pt.y > pivot->pt.y)
{
if (IsRightEdge(*eAlt)) return;
}
Edge* eX = FindLinkingEdge(vAlt, v, (vAlt->pt.y < v->pt.y));
if (!eX)
{
// be very careful creating loose horizontals at minY
if (vAlt->pt.y == v->pt.y && v->pt.y == minY &&
HorizontalBetween(vAlt, v)) return;
eX = CreateEdge(vAlt, v, EdgeKind::loose);
}
CreateTriangle(edge, eAlt, eX);
if (!EdgeCompleted(eX))
DoTriangulateLeft(eX, vAlt, minY);
}
void Delaunay::DoTriangulateRight(Edge* edge, Vertex2* pivot, int64_t minY)
{
// precondition - pivot must be one end of edge (Usually .vB)
//Assert(!EdgeCompleted(edge));
Vertex2* vAlt = nullptr;
Edge* eAlt = nullptr;
Vertex2* v = (edge->vB == pivot) ? edge->vT : edge->vB;
for (auto e : pivot->edges)
{
if (e == edge || !e->isActive) continue;
Vertex2* vX = e->vT == pivot ? e->vB : e->vT;
if (vX == v) continue;
int cps = CrossProductSign(v->pt, pivot->pt, vX->pt);
if (cps == 0) //collinear paths
{
// if pivot is between v and vX then continue;
// nb: this is important for both horiz and non-horiz collinear
if ((v->pt.x > pivot->pt.x) == (pivot->pt.x > vX->pt.x)) continue;
}
// else if right-turning or not the best edge, then continue
else if (cps < 0 || (vAlt && !RightTurning(vX->pt, pivot->pt, vAlt->pt)))
continue;
vAlt = vX;
eAlt = e;
}
if (!vAlt || vAlt->pt.y < minY) return;
// Don't triangulate **across** fixed edges
if (vAlt->pt.y < pivot->pt.y)
{
if (IsRightEdge(*eAlt)) return;
}
else if (vAlt->pt.y > pivot->pt.y)
{
if (IsLeftEdge(*eAlt)) return;
}
Edge* eX = FindLinkingEdge(vAlt, v, (vAlt->pt.y > v->pt.y));
if (!eX)
{
// be very careful creating loose horizontals at minY
if (vAlt->pt.y == v->pt.y && v->pt.y == minY &&
HorizontalBetween(vAlt, v)) return;
eX = CreateEdge(vAlt, v, EdgeKind::loose);
}
CreateTriangle(edge, eX, eAlt);
if (!EdgeCompleted(eX))
DoTriangulateRight(eX, vAlt, minY);
}
void Delaunay::AddEdgeToActives(Edge* edge)
{
// nb: on occassions this method can get called twice for a given edge
// This is because, in the Triangulate() method where vertex 'edges'
// arrays are being parsed, edges can can be removed from the array
// which changes the index of following edges.
if (edge->isActive) return;
edge->prevE = nullptr;
edge->nextE = firstActive;
edge->isActive = true;
if (firstActive)
firstActive->prevE = edge;
firstActive = edge;
}
void Delaunay::RemoveEdgeFromActives(Edge* edge)
{
// first, remove the edge from its vertices
RemoveEdgeFromVertex(edge->vB, edge);
RemoveEdgeFromVertex(edge->vT, edge);
// now remove the edge from double linked list (AEL)
Edge* prev = edge->prevE;
Edge* next = edge->nextE;
if (next) next->prevE = prev;
if (prev) prev->nextE = next;
edge->isActive = false;
if (firstActive == edge) firstActive = next;
}
Paths64 Delaunay::Execute(const Paths64& paths, TriangulateResult& triResult)
{
if (!AddPaths(paths))
{
triResult = TriangulateResult::no_polygons;
return Paths64(); // oops!
}
// if necessary fix path orientation because the algorithm
// expects clockwise outer paths and counter-clockwise inner paths
if (lowermostVertex->innerLM)
{
// the orientation of added paths must be wrong, so
// 1. reverse innerLM flags ...
Vertex2* lm;
while (!locMinStack.empty())
{
lm = locMinStack.top();
lm->innerLM = !lm->innerLM;
locMinStack.pop();
}
// 2. swap edge kinds
for (Edge* e : allEdges)
if (e->kind == EdgeKind::ascend)
e->kind = EdgeKind::descend;
else
e->kind = EdgeKind::ascend;
}
else
{
// path orientation is fine so ...
while (!locMinStack.empty())
locMinStack.pop();
}
std::sort(allEdges.begin(), allEdges.end(), EdgeListSort);
if (!FixupEdgeIntersects())
{
CleanUp();
triResult = TriangulateResult::paths_intersect;
return Paths64(); // oops!
}
std::sort(allVertices.begin(), allVertices.end(), VertexListSort);
MergeDupOrCollinearVertices();
int64_t currY = allVertices[0]->pt.y;
for (auto vIt = allVertices.begin(); vIt != allVertices.end(); ++vIt)
{
Vertex2* v = *vIt;
if (v->edges.empty()) continue;
if (v->pt.y != currY)
{
// JOIN AN INNER LOCMIN WITH A SUITABLE EDGE BELOW
while (!locMinStack.empty())
{
Vertex2* lm = locMinStack.top();
locMinStack.pop();
Edge* e = CreateInnerLocMinLooseEdge(lm);
if (!e)
{
CleanUp();
triResult = TriangulateResult::fail;
return Paths64(); // oops!
}
if (IsHorizontal(*e))
{
if (e->vL == e->vB)
DoTriangulateLeft(e, e->vB, currY); else
DoTriangulateRight(e, e->vB, currY);
}
else
{
DoTriangulateLeft(e, e->vB, currY);
if (!EdgeCompleted(e))
DoTriangulateRight(e, e->vB, currY);
}
// and because adding locMin edges to Actives was delayed ..
AddEdgeToActives(lm->edges[0]);
AddEdgeToActives(lm->edges[1]);
}
while (!horzEdgeStack.empty())
{
Edge* e = horzEdgeStack.top();
horzEdgeStack.pop();
if (EdgeCompleted(e)) continue;
if (e->vB == e->vL) // #45
{
if (IsLeftEdge(*e))
DoTriangulateLeft(e, e->vB, currY);
}
else
if (IsRightEdge(*e))
DoTriangulateRight(e, e->vB, currY);
}
currY = v->pt.y;
}
for (int i = static_cast<int>(v->edges.size()) - 1; i >= 0; --i)
{
// the following line may look superfluous, but within this loop
// v->edges may be altered with additions and or deletions.
// So this line *is* necessary (and why we can't use an iterator).
// Also, we need to use a *descending* index which is safe because
// any additions will be loose edges which are ignored here.
if (i >= static_cast<int>(v->edges.size())) continue;
Edge* e = v->edges[i];
if (EdgeCompleted(e) || IsLooseEdge(*e)) continue;
if (v == e->vB)
{
if (IsHorizontal(*e))
horzEdgeStack.push(e);
// delay adding locMin edges to actives
if (!v->innerLM)
AddEdgeToActives(e);
}
else
{
if (IsHorizontal(*e))
horzEdgeStack.push(e);
else if (IsLeftEdge(*e))
DoTriangulateLeft(e, e->vB, v->pt.y);
else
DoTriangulateRight(e, e->vB, v->pt.y);
}
} // for v->edges loop
if (v->innerLM) locMinStack.push(v);
} // for allVertices loop
while (!horzEdgeStack.empty())
{
Edge* e = horzEdgeStack.top();
horzEdgeStack.pop();
if (!EdgeCompleted(e) && e->vB == e->vL)
DoTriangulateLeft(e, e->vB, currY);
}
if (useDelaunay)
{
// Convert triangles to Delaunay conforming
while (!pendingDelaunayStack.empty())
{
Edge* e = pendingDelaunayStack.top();
pendingDelaunayStack.pop();
ForceLegal(e);
}
}
Paths64 res;
res.reserve(allTriangles.size());
for (auto tri : allTriangles)
{
Path64 p = PathFromTriangle(*tri);
int cps = CrossProductSign(p[0], p[1], p[2]);
if (cps == 0) continue; // skip any empty triangles
if (cps < 0) // ccw
std::reverse(p.begin(), p.end());
res.push_back(p);
}
CleanUp();
triResult = TriangulateResult::success;
return res;
}
static double DistSqr(const Point64& pt1, const Point64& pt2)
{
return Sqr(pt1.x - pt2.x) + Sqr(pt1.y - pt2.y);
}
void Delaunay::AddPath(const Path64& path)
{
size_t len = path.size(), i0 = 0, iPrev, iNext;
// find the first locMin for the current path
if (!FindLocMinIdx(path, len, i0)) return;
iPrev = Prev(i0, len);
while (path[iPrev] == path[i0]) iPrev = Prev(iPrev, len);
iNext = Next(i0, len);
// it is possible for a locMin here to simply be a
// collinear spike that should be ignored, so ...
size_t i = i0;
while (CrossProductSign(path[iPrev], path[i], path[iNext]) == 0)
{
FindLocMinIdx(path, len, i);
if (i == i0) return; // this is an entirely collinear path
iPrev = Prev(i, len);
while (path[iPrev] == path[i]) iPrev = Prev(iPrev, len);
iNext = Next(i, len);
}
size_t vert_cnt = allVertices.size();
// we are now at the first legitimate locMin
Vertex2* v0 = allVertices.emplace_back(new Vertex2(path[i]));
if (LeftTurning(path[iPrev], path[i], path[iNext]))
v0->innerLM = true;
Vertex2* vPrev = v0;
i = iNext;
for (;;)
{
// vPrev is a locMin here
locMinStack.push(vPrev);
// ? update lowermostVertex ...
if (!lowermostVertex ||
vPrev->pt.y > lowermostVertex->pt.y ||
(vPrev->pt.y == lowermostVertex->pt.y &&
vPrev->pt.x < lowermostVertex->pt.x))
lowermostVertex = vPrev;
iNext = Next(i, len);
if (CrossProductSign(vPrev->pt, path[i], path[iNext]) == 0)
{
i = iNext;
continue;
}
// ascend up next bound to LocMax
while (path[i].y <= vPrev->pt.y)
{
Vertex2* v = allVertices.emplace_back(new Vertex2(path[i]));
CreateEdge(vPrev, v, EdgeKind::ascend);
vPrev = v;
i = iNext;
iNext = Next(i, len);
while (CrossProductSign(vPrev->pt, path[i], path[iNext]) == 0)
{
i = iNext;
iNext = Next(i, len);
}
}
// Now at a locMax, so descend to next locMin
Vertex2* vPrevPrev = vPrev;
while (i != i0 && path[i].y >= vPrev->pt.y)
{
Vertex2* v = allVertices.emplace_back(new Vertex2(path[i]));
CreateEdge(v, vPrev, EdgeKind::descend);
vPrevPrev = vPrev;
vPrev = v;
i = iNext;
iNext = Next(i, len);
while (CrossProductSign(vPrev->pt, path[i], path[iNext]) == 0)
{
i = iNext;
iNext = Next(i, len);
}
}
// now at the next locMin
if (i == i0) break; // break for(;;) loop
if (LeftTurning(vPrevPrev->pt, vPrev->pt, path[i]))
vPrev->innerLM = true;
}
CreateEdge(v0, vPrev, EdgeKind::descend);
// finally, ignore this path if is not a polygon or too small
len = allVertices.size() - vert_cnt;
i = vert_cnt;
if (len < 3 || (len == 3 && // or just a very tiny triangle
((DistSqr(allVertices[i]->pt, allVertices[i + 1]->pt) <= 1) ||
(DistSqr(allVertices[i + 1]->pt, allVertices[i + 2]->pt) <= 1) ||
(DistSqr(allVertices[i + 2]->pt, allVertices[i]->pt) <= 1))))
{
for (size_t j = vert_cnt; j < allVertices.size(); ++j)
allVertices[j]->edges.clear(); // flag to ignore
}
}
bool Delaunay::AddPaths(const Paths64& paths)
{
const auto total_vertex_count =
std::accumulate(paths.begin(), paths.end(), size_t(0),
[](const auto& a, const Path64& path)
{return a + path.size(); });
if (total_vertex_count == 0) return false;
allVertices.reserve(allVertices.capacity() + total_vertex_count);
allEdges.reserve(allEdges.capacity() + total_vertex_count);
for (const Path64& path : paths)
AddPath(path);
return (allVertices.size() > 2);
}
TriangulateResult Triangulate(const Paths64& pp, Paths64& solution, bool useDelaunay)
{
TriangulateResult result;
Delaunay d(useDelaunay);
solution = d.Execute(pp, result);
return result;
}
TriangulateResult Triangulate(const PathsD& pp, int decPlaces, PathsD& solution, bool useDelaunay)
{
int ec;
double scale;
TriangulateResult result;
if (decPlaces <= 0) scale = 1;
else if (decPlaces > 8) scale = std::pow(10, 8);
else scale = std::pow(10, decPlaces);
Paths64 pp64 = ScalePaths<int64_t, double>(pp, scale, ec);
Delaunay d(useDelaunay);
Paths64 sol64 = d.Execute(pp64, result);
solution = ScalePaths<double, int64_t>(sol64, 1 / scale, ec);
return result;
}
} // Clipper2Lib namespace
Reference in New Issue View Git Blame Copy Permalink