copilot/infinipaint
1
0
Fork 0
Code Issues Pull requests Projects Releases 1 Packages Wiki Activity Actions 1

Using a proper polygon triangulation library instead of delaunator for lasso tool

Browse source
This commit is contained in:
Yousef Khadadeh 2025-08-02 14:42:53 +03:00
parent a4ba969e11
commit d07eaf4207
5 changed files with 879 additions and 580 deletions
Download patch file Download diff file Expand all files Collapse all files

561
include/delaunator.hpp
View file

@ -1,561 +0,0 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <exception>
#include <iostream>
#include <limits>
#include <memory>
#include <utility>
#include <vector>
namespace delaunator {
// Kahan and Babuska summation, Neumaier variant; accumulates less FP error
inline double sum(const std::vector<double>& x) {
double sum = x[0];
double err = 0.0;
for (size_t i = 1; i < x.size(); i++) {
const double k = x[i];
const double m = sum + k;
err += std::fabs(sum) >= std::fabs(k) ? sum - m + k : k - m + sum;
sum = m;
}
return sum + err;
}
inline double dist(
const double ax,
const double ay,
const double bx,
const double by) {
const double dx = ax - bx;
const double dy = ay - by;
return dx * dx + dy * dy;
}
inline double circumradius(
const double ax,
const double ay,
const double bx,
const double by,
const double cx,
const double cy) {
const double dx = bx - ax;
const double dy = by - ay;
const double ex = cx - ax;
const double ey = cy - ay;
const double bl = dx * dx + dy * dy;
const double cl = ex * ex + ey * ey;
const double d = dx * ey - dy * ex;
const double x = (ey * bl - dy * cl) * 0.5 / d;
const double y = (dx * cl - ex * bl) * 0.5 / d;
if ((bl > 0.0 || bl < 0.0) && (cl > 0.0 || cl < 0.0) && (d > 0.0 || d < 0.0)) {
return x * x + y * y;
} else {
return std::numeric_limits<double>::max();
}
}
inline bool orient(
const double px,
const double py,
const double qx,
const double qy,
const double rx,
const double ry) {
return (qy - py) * (rx - qx) - (qx - px) * (ry - qy) < 0.0;
}
inline std::pair<double, double> circumcenter(
const double ax,
const double ay,
const double bx,
const double by,
const double cx,
const double cy) {
const double dx = bx - ax;
const double dy = by - ay;
const double ex = cx - ax;
const double ey = cy - ay;
const double bl = dx * dx + dy * dy;
const double cl = ex * ex + ey * ey;
const double d = dx * ey - dy * ex;
const double x = ax + (ey * bl - dy * cl) * 0.5 / d;
const double y = ay + (dx * cl - ex * bl) * 0.5 / d;
return std::make_pair(x, y);
}
inline double compare(
std::vector<double> const& coords,
std::size_t i,
std::size_t j,
double cx,
double cy) {
const double d1 = dist(coords[2 * i], coords[2 * i + 1], cx, cy);
const double d2 = dist(coords[2 * j], coords[2 * j + 1], cx, cy);
const double diff1 = d1 - d2;
const double diff2 = coords[2 * i] - coords[2 * j];
const double diff3 = coords[2 * i + 1] - coords[2 * j + 1];
if (diff1 > 0.0 || diff1 < 0.0) {
return diff1;
} else if (diff2 > 0.0 || diff2 < 0.0) {
return diff2;
} else {
return diff3;
}
}
struct sort_to_center {
std::vector<double> const& coords;
double cx;
double cy;
bool operator()(std::size_t i, std::size_t j) {
return compare(coords, i, j, cx, cy) < 0;
}
};
inline bool in_circle(
double ax,
double ay,
double bx,
double by,
double cx,
double cy,
double px,
double py) {
const double dx = ax - px;
const double dy = ay - py;
const double ex = bx - px;
const double ey = by - py;
const double fx = cx - px;
const double fy = cy - py;
const double ap = dx * dx + dy * dy;
const double bp = ex * ex + ey * ey;
const double cp = fx * fx + fy * fy;
return (dx * (ey * cp - bp * fy) -
dy * (ex * cp - bp * fx) +
ap * (ex * fy - ey * fx)) < 0.0;
}
constexpr double EPSILON = std::numeric_limits<double>::epsilon();
constexpr std::size_t INVALID_INDEX = std::numeric_limits<std::size_t>::max();
inline bool check_pts_equal(double x1, double y1, double x2, double y2) {
return std::fabs(x1 - x2) <= EPSILON &&
std::fabs(y1 - y2) <= EPSILON;
}
// monotonically increases with real angle, but doesn't need expensive trigonometry
inline double pseudo_angle(double dx, double dy) {
const double p = dx / (std::abs(dx) + std::abs(dy));
return (dy > 0.0 ? 3.0 - p : 1.0 + p) / 4.0; // [0..1)
}
struct DelaunatorPoint {
std::size_t i;
double x;
double y;
std::size_t t;
std::size_t prev;
std::size_t next;
bool removed;
};
class Delaunator {
public:
std::vector<double> const& coords;
std::vector<std::size_t> triangles;
std::vector<std::size_t> halfedges;
std::vector<std::size_t> hull_prev;
std::vector<std::size_t> hull_next;
std::vector<std::size_t> hull_tri;
std::size_t hull_start;
Delaunator(std::vector<double> const& in_coords);
double get_hull_area();
private:
std::vector<std::size_t> m_hash;
double m_center_x;
double m_center_y;
std::size_t m_hash_size;
std::size_t legalize(std::size_t a);
std::size_t hash_key(double x, double y);
std::size_t add_triangle(
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t a,
std::size_t b,
std::size_t c);
void link(std::size_t a, std::size_t b);
};
Delaunator::Delaunator(std::vector<double> const& in_coords)
: coords(in_coords),
triangles(),
halfedges(),
hull_prev(),
hull_next(),
hull_tri(),
hull_start(),
m_hash(),
m_center_x(),
m_center_y(),
m_hash_size() {
std::size_t n = coords.size() >> 1;
double max_x = std::numeric_limits<double>::min();
double max_y = std::numeric_limits<double>::min();
double min_x = std::numeric_limits<double>::max();
double min_y = std::numeric_limits<double>::max();
std::vector<std::size_t> ids;
ids.reserve(n);
for (std::size_t i = 0; i < n; i++) {
const double x = coords[2 * i];
const double y = coords[2 * i + 1];
if (x < min_x) min_x = x;
if (y < min_y) min_y = y;
if (x > max_x) max_x = x;
if (y > max_y) max_y = y;
ids.push_back(i);
}
const double cx = (min_x + max_x) / 2;
const double cy = (min_y + max_y) / 2;
double min_dist = std::numeric_limits<double>::max();
std::size_t i0 = INVALID_INDEX;
std::size_t i1 = INVALID_INDEX;
std::size_t i2 = INVALID_INDEX;
// pick a seed point close to the centroid
for (std::size_t i = 0; i < n; i++) {
const double d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]);
if (d < min_dist) {
i0 = i;
min_dist = d;
}
}
const double i0x = coords[2 * i0];
const double i0y = coords[2 * i0 + 1];
min_dist = std::numeric_limits<double>::max();
// find the point closest to the seed
for (std::size_t i = 0; i < n; i++) {
if (i == i0) continue;
const double d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1]);
if (d < min_dist && d > 0.0) {
i1 = i;
min_dist = d;
}
}
double i1x = coords[2 * i1];
double i1y = coords[2 * i1 + 1];
double min_radius = std::numeric_limits<double>::max();
// find the third point which forms the smallest circumcircle with the first two
for (std::size_t i = 0; i < n; i++) {
if (i == i0 || i == i1) continue;
const double r = circumradius(
i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1]);
if (r < min_radius) {
i2 = i;
min_radius = r;
}
}
if (!(min_radius < std::numeric_limits<double>::max())) {
throw std::runtime_error("not triangulation");
}
double i2x = coords[2 * i2];
double i2y = coords[2 * i2 + 1];
if (orient(i0x, i0y, i1x, i1y, i2x, i2y)) {
std::swap(i1, i2);
std::swap(i1x, i2x);
std::swap(i1y, i2y);
}
std::tie(m_center_x, m_center_y) = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);
// sort the points by distance from the seed triangle circumcenter
std::sort(ids.begin(), ids.end(), sort_to_center{ coords, m_center_x, m_center_y });
// initialize a hash table for storing edges of the advancing convex hull
m_hash_size = static_cast<std::size_t>(std::llround(std::ceil(std::sqrt(n))));
m_hash.resize(m_hash_size);
std::fill(m_hash.begin(), m_hash.end(), INVALID_INDEX);
// initialize arrays for tracking the edges of the advancing convex hull
hull_prev.resize(n);
hull_next.resize(n);
hull_tri.resize(n);
hull_start = i0;
size_t hull_size = 3;
hull_next[i0] = hull_prev[i2] = i1;
hull_next[i1] = hull_prev[i0] = i2;
hull_next[i2] = hull_prev[i1] = i0;
hull_tri[i0] = 0;
hull_tri[i1] = 1;
hull_tri[i2] = 2;
m_hash[hash_key(i0x, i0y)] = i0;
m_hash[hash_key(i1x, i1y)] = i1;
m_hash[hash_key(i2x, i2y)] = i2;
std::size_t max_triangles = n < 3 ? 1 : 2 * n - 5;
triangles.reserve(max_triangles * 3);
halfedges.reserve(max_triangles * 3);
add_triangle(i0, i1, i2, INVALID_INDEX, INVALID_INDEX, INVALID_INDEX);
double xp = std::numeric_limits<double>::quiet_NaN();
double yp = std::numeric_limits<double>::quiet_NaN();
for (std::size_t k = 0; k < n; k++) {
const std::size_t i = ids[k];
const double x = coords[2 * i];
const double y = coords[2 * i + 1];
// skip near-duplicate points
if (k > 0 && check_pts_equal(x, y, xp, yp)) continue;
xp = x;
yp = y;
// skip seed triangle points
if (
check_pts_equal(x, y, i0x, i0y) ||
check_pts_equal(x, y, i1x, i1y) ||
check_pts_equal(x, y, i2x, i2y)) continue;
// find a visible edge on the convex hull using edge hash
std::size_t start = 0;
size_t key = hash_key(x, y);
for (size_t j = 0; j < m_hash_size; j++) {
start = m_hash[(key + j) % m_hash_size];
if (start != INVALID_INDEX && start != hull_next[start]) break;
}
start = hull_prev[start];
size_t e = start;
size_t q;
while (q = hull_next[e], !orient(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1])) { //TODO: does it works in a same way as in JS
e = q;
if (e == start) {
e = INVALID_INDEX;
break;
}
}
if (e == INVALID_INDEX) continue; // likely a near-duplicate point; skip it
// add the first triangle from the point
std::size_t t = add_triangle(
e,
i,
hull_next[e],
INVALID_INDEX,
INVALID_INDEX,
hull_tri[e]);
hull_tri[i] = legalize(t + 2);
hull_tri[e] = t;
hull_size++;
// walk forward through the hull, adding more triangles and flipping recursively
std::size_t next = hull_next[e];
while (
q = hull_next[next],
orient(x, y, coords[2 * next], coords[2 * next + 1], coords[2 * q], coords[2 * q + 1])) {
t = add_triangle(next, i, q, hull_tri[i], INVALID_INDEX, hull_tri[next]);
hull_tri[i] = legalize(t + 2);
hull_next[next] = next; // mark as removed
hull_size--;
next = q;
}
// walk backward from the other side, adding more triangles and flipping
if (e == start) {
while (
q = hull_prev[e],
orient(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1])) {
t = add_triangle(q, i, e, INVALID_INDEX, hull_tri[e], hull_tri[q]);
legalize(t + 2);
hull_tri[q] = t;
hull_next[e] = e; // mark as removed
hull_size--;
e = q;
}
}
// update the hull indices
hull_prev[i] = e;
hull_start = e;
hull_prev[next] = i;
hull_next[e] = i;
hull_next[i] = next;
m_hash[hash_key(x, y)] = i;
m_hash[hash_key(coords[2 * e], coords[2 * e + 1])] = e;
}
}
double Delaunator::get_hull_area() {
std::vector<double> hull_area;
size_t e = hull_start;
do {
hull_area.push_back((coords[2 * e] - coords[2 * hull_prev[e]]) * (coords[2 * e + 1] + coords[2 * hull_prev[e] + 1]));
e = hull_next[e];
} while (e != hull_start);
return sum(hull_area);
}
std::size_t Delaunator::legalize(std::size_t a) {
const std::size_t b = halfedges[a];
/* if the pair of triangles doesn't satisfy the Delaunay condition
* (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
* then do the same check/flip recursively for the new pair of triangles
*
* pl pl
* /||\ / \
* al/ || \bl al/ \a
* / || \ / \
* / a||b \ flip /___ar___\
* p0\ || /p1 => p0\---bl---/p1
* \ || / \ /
* ar\ || /br b\ /br
* \||/ \ /
* pr pr
*/
const std::size_t a0 = a - a % 3;
const std::size_t b0 = b - b % 3;
const std::size_t al = a0 + (a + 1) % 3;
const std::size_t ar = a0 + (a + 2) % 3;
const std::size_t bl = b0 + (b + 2) % 3;
const std::size_t p0 = triangles[ar];
const std::size_t pr = triangles[a];
const std::size_t pl = triangles[al];
const std::size_t p1 = triangles[bl];
if (b == INVALID_INDEX) {
return ar;
}
const bool illegal = in_circle(
coords[2 * p0],
coords[2 * p0 + 1],
coords[2 * pr],
coords[2 * pr + 1],
coords[2 * pl],
coords[2 * pl + 1],
coords[2 * p1],
coords[2 * p1 + 1]);
if (illegal) {
triangles[a] = p1;
triangles[b] = p0;
auto hbl = halfedges[bl];
// edge swapped on the other side of the hull (rare); fix the halfedge reference
if (hbl == INVALID_INDEX) {
std::size_t e = hull_start;
do {
if (hull_tri[e] == bl) {
hull_tri[e] = a;
break;
}
e = hull_next[e];
} while (e != hull_start);
}
link(a, hbl);
link(b, halfedges[ar]);
link(ar, bl);
std::size_t br = b0 + (b + 1) % 3;
legalize(a);
return legalize(br);
}
return ar;
}
std::size_t Delaunator::hash_key(double x, double y) {
const double dx = x - m_center_x;
const double dy = y - m_center_y;
return static_cast<std::size_t>(std::llround(
std::floor(pseudo_angle(dx, dy) * static_cast<double>(m_hash_size)))) %
m_hash_size;
}
std::size_t Delaunator::add_triangle(
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t a,
std::size_t b,
std::size_t c) {
std::size_t t = triangles.size();
triangles.push_back(i0);
triangles.push_back(i1);
triangles.push_back(i2);
link(t, a);
link(t + 1, b);
link(t + 2, c);
return t;
}
void Delaunator::link(std::size_t a, std::size_t b) {
std::size_t s = halfedges.size();
if (a == s) {
halfedges.push_back(b);
} else if (a < s) {
halfedges[a] = b;
} else {
throw std::runtime_error("Cannot link edge");
}
if (b != INVALID_INDEX) {
std::size_t s2 = halfedges.size();
if (b == s2) {
halfedges.push_back(a);
} else if (b < s2) {
halfedges[b] = a;
} else {
throw std::runtime_error("Cannot link edge");
}
}
}
} //namespace delaunator

814
include/earcut.hpp Normal file
View file

@ -0,0 +1,814 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <memory>
#include <utility>
#include <vector>
namespace mapbox {
namespace util {
template <std::size_t I, typename T> struct nth {
inline static typename std::tuple_element<I, T>::type
get(const T& t) { return std::get<I>(t); };
};
}
namespace detail {
template <typename N = uint32_t>
class Earcut {
public:
std::vector<N> indices;
std::size_t vertices = 0;
template <typename Polygon>
void operator()(const Polygon& points);
private:
struct Node {
Node(N index, double x_, double y_) : i(index), x(x_), y(y_) {}
Node(const Node&) = delete;
Node& operator=(const Node&) = delete;
Node(Node&&) = delete;
Node& operator=(Node&&) = delete;
const N i;
const double x;
const double y;
// previous and next vertice nodes in a polygon ring
Node* prev = nullptr;
Node* next = nullptr;
// z-order curve value
int32_t z = 0;
// previous and next nodes in z-order
Node* prevZ = nullptr;
Node* nextZ = nullptr;
// indicates whether this is a steiner point
bool steiner = false;
};
template <typename Ring> Node* linkedList(const Ring& points, const bool clockwise);
Node* filterPoints(Node* start, Node* end = nullptr);
void earcutLinked(Node* ear, int pass = 0);
bool isEar(Node* ear);
bool isEarHashed(Node* ear);
Node* cureLocalIntersections(Node* start);
void splitEarcut(Node* start);
template <typename Polygon> Node* eliminateHoles(const Polygon& points, Node* outerNode);
Node* eliminateHole(Node* hole, Node* outerNode);
Node* findHoleBridge(Node* hole, Node* outerNode);
bool sectorContainsSector(const Node* m, const Node* p);
void indexCurve(Node* start);
Node* sortLinked(Node* list);
int32_t zOrder(const double x_, const double y_);
Node* getLeftmost(Node* start);
bool pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const;
bool isValidDiagonal(Node* a, Node* b);
double area(const Node* p, const Node* q, const Node* r) const;
bool equals(const Node* p1, const Node* p2);
bool intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2);
bool onSegment(const Node* p, const Node* q, const Node* r);
int sign(double val);
bool intersectsPolygon(const Node* a, const Node* b);
bool locallyInside(const Node* a, const Node* b);
bool middleInside(const Node* a, const Node* b);
Node* splitPolygon(Node* a, Node* b);
template <typename Point> Node* insertNode(std::size_t i, const Point& p, Node* last);
void removeNode(Node* p);
bool hashing;
double minX, maxX;
double minY, maxY;
double inv_size = 0;
template <typename T, typename Alloc = std::allocator<T>>
class ObjectPool {
public:
ObjectPool() { }
ObjectPool(std::size_t blockSize_) {
reset(blockSize_);
}
~ObjectPool() {
clear();
}
template <typename... Args>
T* construct(Args&&... args) {
if (currentIndex >= blockSize) {
currentBlock = alloc_traits::allocate(alloc, blockSize);
allocations.emplace_back(currentBlock);
currentIndex = 0;
}
T* object = &currentBlock[currentIndex++];
alloc_traits::construct(alloc, object, std::forward<Args>(args)...);
return object;
}
void reset(std::size_t newBlockSize) {
for (auto allocation : allocations) {
alloc_traits::deallocate(alloc, allocation, blockSize);
}
allocations.clear();
blockSize = std::max<std::size_t>(1, newBlockSize);
currentBlock = nullptr;
currentIndex = blockSize;
}
void clear() { reset(blockSize); }
private:
T* currentBlock = nullptr;
std::size_t currentIndex = 1;
std::size_t blockSize = 1;
std::vector<T*> allocations;
Alloc alloc;
typedef typename std::allocator_traits<Alloc> alloc_traits;
};
ObjectPool<Node> nodes;
};
template <typename N> template <typename Polygon>
void Earcut<N>::operator()(const Polygon& points) {
// reset
indices.clear();
vertices = 0;
if (points.empty()) return;
double x;
double y;
int threshold = 80;
std::size_t len = 0;
for (size_t i = 0; threshold >= 0 && i < points.size(); i++) {
threshold -= static_cast<int>(points[i].size());
len += points[i].size();
}
//estimate size of nodes and indices
nodes.reset(len * 3 / 2);
indices.reserve(len + points[0].size());
Node* outerNode = linkedList(points[0], true);
if (!outerNode || outerNode->prev == outerNode->next) return;
if (points.size() > 1) outerNode = eliminateHoles(points, outerNode);
// if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox
hashing = threshold < 0;
if (hashing) {
Node* p = outerNode->next;
minX = maxX = outerNode->x;
minY = maxY = outerNode->y;
do {
x = p->x;
y = p->y;
minX = std::min<double>(minX, x);
minY = std::min<double>(minY, y);
maxX = std::max<double>(maxX, x);
maxY = std::max<double>(maxY, y);
p = p->next;
} while (p != outerNode);
// minX, minY and inv_size are later used to transform coords into integers for z-order calculation
inv_size = std::max<double>(maxX - minX, maxY - minY);
inv_size = inv_size != .0 ? (32767. / inv_size) : .0;
}
earcutLinked(outerNode);
nodes.clear();
}
// create a circular doubly linked list from polygon points in the specified winding order
template <typename N> template <typename Ring>
typename Earcut<N>::Node*
Earcut<N>::linkedList(const Ring& points, const bool clockwise) {
using Point = typename Ring::value_type;
double sum = 0;
const std::size_t len = points.size();
std::size_t i, j;
Node* last = nullptr;
// calculate original winding order of a polygon ring
for (i = 0, j = len > 0 ? len - 1 : 0; i < len; j = i++) {
const auto& p1 = points[i];
const auto& p2 = points[j];
const double p20 = util::nth<0, Point>::get(p2);
const double p10 = util::nth<0, Point>::get(p1);
const double p11 = util::nth<1, Point>::get(p1);
const double p21 = util::nth<1, Point>::get(p2);
sum += (p20 - p10) * (p11 + p21);
}
// link points into circular doubly-linked list in the specified winding order
if (clockwise == (sum > 0)) {
for (i = 0; i < len; i++) last = insertNode(vertices + i, points[i], last);
} else {
for (i = len; i-- > 0;) last = insertNode(vertices + i, points[i], last);
}
if (last && equals(last, last->next)) {
removeNode(last);
last = last->next;
}
vertices += len;
return last;
}
// eliminate colinear or duplicate points
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::filterPoints(Node* start, Node* end) {
if (!end) end = start;
Node* p = start;
bool again;
do {
again = false;
if (!p->steiner && (equals(p, p->next) || area(p->prev, p, p->next) == 0)) {
removeNode(p);
p = end = p->prev;
if (p == p->next) break;
again = true;
} else {
p = p->next;
}
} while (again || p != end);
return end;
}
// main ear slicing loop which triangulates a polygon (given as a linked list)
template <typename N>
void Earcut<N>::earcutLinked(Node* ear, int pass) {
if (!ear) return;
// interlink polygon nodes in z-order
if (!pass && hashing) indexCurve(ear);
Node* stop = ear;
Node* prev;
Node* next;
// iterate through ears, slicing them one by one
while (ear->prev != ear->next) {
prev = ear->prev;
next = ear->next;
if (hashing ? isEarHashed(ear) : isEar(ear)) {
// cut off the triangle
indices.emplace_back(prev->i);
indices.emplace_back(ear->i);
indices.emplace_back(next->i);
removeNode(ear);
// skipping the next vertice leads to less sliver triangles
ear = next->next;
stop = next->next;
continue;
}
ear = next;
// if we looped through the whole remaining polygon and can't find any more ears
if (ear == stop) {
// try filtering points and slicing again
if (!pass) earcutLinked(filterPoints(ear), 1);
// if this didn't work, try curing all small self-intersections locally
else if (pass == 1) {
ear = cureLocalIntersections(filterPoints(ear));
earcutLinked(ear, 2);
// as a last resort, try splitting the remaining polygon into two
} else if (pass == 2) splitEarcut(ear);
break;
}
}
}
// check whether a polygon node forms a valid ear with adjacent nodes
template <typename N>
bool Earcut<N>::isEar(Node* ear) {
const Node* a = ear->prev;
const Node* b = ear;
const Node* c = ear->next;
if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
// now make sure we don't have other points inside the potential ear
Node* p = ear->next->next;
while (p != ear->prev) {
if (pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->next;
}
return true;
}
template <typename N>
bool Earcut<N>::isEarHashed(Node* ear) {
const Node* a = ear->prev;
const Node* b = ear;
const Node* c = ear->next;
if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
// triangle bbox; min & max are calculated like this for speed
const double minTX = std::min<double>(a->x, std::min<double>(b->x, c->x));
const double minTY = std::min<double>(a->y, std::min<double>(b->y, c->y));
const double maxTX = std::max<double>(a->x, std::max<double>(b->x, c->x));
const double maxTY = std::max<double>(a->y, std::max<double>(b->y, c->y));
// z-order range for the current triangle bbox;
const int32_t minZ = zOrder(minTX, minTY);
const int32_t maxZ = zOrder(maxTX, maxTY);
// first look for points inside the triangle in increasing z-order
Node* p = ear->nextZ;
while (p && p->z <= maxZ) {
if (p != ear->prev && p != ear->next &&
pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->nextZ;
}
// then look for points in decreasing z-order
p = ear->prevZ;
while (p && p->z >= minZ) {
if (p != ear->prev && p != ear->next &&
pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->prevZ;
}
return true;
}
// go through all polygon nodes and cure small local self-intersections
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::cureLocalIntersections(Node* start) {
Node* p = start;
do {
Node* a = p->prev;
Node* b = p->next->next;
// a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2])
if (!equals(a, b) && intersects(a, p, p->next, b) && locallyInside(a, b) && locallyInside(b, a)) {
indices.emplace_back(a->i);
indices.emplace_back(p->i);
indices.emplace_back(b->i);
// remove two nodes involved
removeNode(p);
removeNode(p->next);
p = start = b;
}
p = p->next;
} while (p != start);
return filterPoints(p);
}
// try splitting polygon into two and triangulate them independently
template <typename N>
void Earcut<N>::splitEarcut(Node* start) {
// look for a valid diagonal that divides the polygon into two
Node* a = start;
do {
Node* b = a->next->next;
while (b != a->prev) {
if (a->i != b->i && isValidDiagonal(a, b)) {
// split the polygon in two by the diagonal
Node* c = splitPolygon(a, b);
// filter colinear points around the cuts
a = filterPoints(a, a->next);
c = filterPoints(c, c->next);
// run earcut on each half
earcutLinked(a);
earcutLinked(c);
return;
}
b = b->next;
}
a = a->next;
} while (a != start);
}
// link every hole into the outer loop, producing a single-ring polygon without holes
template <typename N> template <typename Polygon>
typename Earcut<N>::Node*
Earcut<N>::eliminateHoles(const Polygon& points, Node* outerNode) {
const size_t len = points.size();
std::vector<Node*> queue;
for (size_t i = 1; i < len; i++) {
Node* list = linkedList(points[i], false);
if (list) {
if (list == list->next) list->steiner = true;
queue.push_back(getLeftmost(list));
}
}
std::sort(queue.begin(), queue.end(), [](const Node* a, const Node* b) {
return a->x < b->x;
});
// process holes from left to right
for (size_t i = 0; i < queue.size(); i++) {
outerNode = eliminateHole(queue[i], outerNode);
}
return outerNode;
}
// find a bridge between vertices that connects hole with an outer ring and and link it
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::eliminateHole(Node* hole, Node* outerNode) {
Node* bridge = findHoleBridge(hole, outerNode);
if (!bridge) {
return outerNode;
}
Node* bridgeReverse = splitPolygon(bridge, hole);
// filter collinear points around the cuts
filterPoints(bridgeReverse, bridgeReverse->next);
// Check if input node was removed by the filtering
return filterPoints(bridge, bridge->next);
}
// David Eberly's algorithm for finding a bridge between hole and outer polygon
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::findHoleBridge(Node* hole, Node* outerNode) {
Node* p = outerNode;
double hx = hole->x;
double hy = hole->y;
double qx = -std::numeric_limits<double>::infinity();
Node* m = nullptr;
// find a segment intersected by a ray from the hole's leftmost Vertex to the left;
// segment's endpoint with lesser x will be potential connection Vertex
do {
if (hy <= p->y && hy >= p->next->y && p->next->y != p->y) {
double x = p->x + (hy - p->y) * (p->next->x - p->x) / (p->next->y - p->y);
if (x <= hx && x > qx) {
qx = x;
m = p->x < p->next->x ? p : p->next;
if (x == hx) return m; // hole touches outer segment; pick leftmost endpoint
}
}
p = p->next;
} while (p != outerNode);
if (!m) return 0;
// look for points inside the triangle of hole Vertex, segment intersection and endpoint;
// if there are no points found, we have a valid connection;
// otherwise choose the Vertex of the minimum angle with the ray as connection Vertex
const Node* stop = m;
double tanMin = std::numeric_limits<double>::infinity();
double tanCur = 0;
p = m;
double mx = m->x;
double my = m->y;
do {
if (hx >= p->x && p->x >= mx && hx != p->x &&
pointInTriangle(hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy, p->x, p->y)) {
tanCur = std::abs(hy - p->y) / (hx - p->x); // tangential
if (locallyInside(p, hole) &&
(tanCur < tanMin || (tanCur == tanMin && (p->x > m->x || sectorContainsSector(m, p))))) {
m = p;
tanMin = tanCur;
}
}
p = p->next;
} while (p != stop);
return m;
}
// whether sector in vertex m contains sector in vertex p in the same coordinates
template <typename N>
bool Earcut<N>::sectorContainsSector(const Node* m, const Node* p) {
return area(m->prev, m, p->prev) < 0 && area(p->next, m, m->next) < 0;
}
// interlink polygon nodes in z-order
template <typename N>
void Earcut<N>::indexCurve(Node* start) {
assert(start);
Node* p = start;
do {
p->z = p->z ? p->z : zOrder(p->x, p->y);
p->prevZ = p->prev;
p->nextZ = p->next;
p = p->next;
} while (p != start);
p->prevZ->nextZ = nullptr;
p->prevZ = nullptr;
sortLinked(p);
}
// Simon Tatham's linked list merge sort algorithm
// http://www.chiark.greenend.org.uk/~sgtatham/algorithms/listsort.html
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::sortLinked(Node* list) {
assert(list);
Node* p;
Node* q;
Node* e;
Node* tail;
int i, numMerges, pSize, qSize;
int inSize = 1;
for (;;) {
p = list;
list = nullptr;
tail = nullptr;
numMerges = 0;
while (p) {
numMerges++;
q = p;
pSize = 0;
for (i = 0; i < inSize; i++) {
pSize++;
q = q->nextZ;
if (!q) break;
}
qSize = inSize;
while (pSize > 0 || (qSize > 0 && q)) {
if (pSize == 0) {
e = q;
q = q->nextZ;
qSize--;
} else if (qSize == 0 || !q) {
e = p;
p = p->nextZ;
pSize--;
} else if (p->z <= q->z) {
e = p;
p = p->nextZ;
pSize--;
} else {
e = q;
q = q->nextZ;
qSize--;
}
if (tail) tail->nextZ = e;
else list = e;
e->prevZ = tail;
tail = e;
}
p = q;
}
tail->nextZ = nullptr;
if (numMerges <= 1) return list;
inSize *= 2;
}
}
// z-order of a Vertex given coords and size of the data bounding box
template <typename N>
int32_t Earcut<N>::zOrder(const double x_, const double y_) {
// coords are transformed into non-negative 15-bit integer range
int32_t x = static_cast<int32_t>((x_ - minX) * inv_size);
int32_t y = static_cast<int32_t>((y_ - minY) * inv_size);
x = (x | (x << 8)) & 0x00FF00FF;
x = (x | (x << 4)) & 0x0F0F0F0F;
x = (x | (x << 2)) & 0x33333333;
x = (x | (x << 1)) & 0x55555555;
y = (y | (y << 8)) & 0x00FF00FF;
y = (y | (y << 4)) & 0x0F0F0F0F;
y = (y | (y << 2)) & 0x33333333;
y = (y | (y << 1)) & 0x55555555;
return x | (y << 1);
}
// find the leftmost node of a polygon ring
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::getLeftmost(Node* start) {
Node* p = start;
Node* leftmost = start;
do {
if (p->x < leftmost->x || (p->x == leftmost->x && p->y < leftmost->y))
leftmost = p;
p = p->next;
} while (p != start);
return leftmost;
}
// check if a point lies within a convex triangle
template <typename N>
bool Earcut<N>::pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const {
return (cx - px) * (ay - py) >= (ax - px) * (cy - py) &&
(ax - px) * (by - py) >= (bx - px) * (ay - py) &&
(bx - px) * (cy - py) >= (cx - px) * (by - py);
}
// check if a diagonal between two polygon nodes is valid (lies in polygon interior)
template <typename N>
bool Earcut<N>::isValidDiagonal(Node* a, Node* b) {
return a->next->i != b->i && a->prev->i != b->i && !intersectsPolygon(a, b) && // dones't intersect other edges
((locallyInside(a, b) && locallyInside(b, a) && middleInside(a, b) && // locally visible
(area(a->prev, a, b->prev) != 0.0 || area(a, b->prev, b) != 0.0)) || // does not create opposite-facing sectors
(equals(a, b) && area(a->prev, a, a->next) > 0 && area(b->prev, b, b->next) > 0)); // special zero-length case
}
// signed area of a triangle
template <typename N>
double Earcut<N>::area(const Node* p, const Node* q, const Node* r) const {
return (q->y - p->y) * (r->x - q->x) - (q->x - p->x) * (r->y - q->y);
}
// check if two points are equal
template <typename N>
bool Earcut<N>::equals(const Node* p1, const Node* p2) {
return p1->x == p2->x && p1->y == p2->y;
}
// check if two segments intersect
template <typename N>
bool Earcut<N>::intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2) {
int o1 = sign(area(p1, q1, p2));
int o2 = sign(area(p1, q1, q2));
int o3 = sign(area(p2, q2, p1));
int o4 = sign(area(p2, q2, q1));
if (o1 != o2 && o3 != o4) return true; // general case
if (o1 == 0 && onSegment(p1, p2, q1)) return true; // p1, q1 and p2 are collinear and p2 lies on p1q1
if (o2 == 0 && onSegment(p1, q2, q1)) return true; // p1, q1 and q2 are collinear and q2 lies on p1q1
if (o3 == 0 && onSegment(p2, p1, q2)) return true; // p2, q2 and p1 are collinear and p1 lies on p2q2
if (o4 == 0 && onSegment(p2, q1, q2)) return true; // p2, q2 and q1 are collinear and q1 lies on p2q2
return false;
}
// for collinear points p, q, r, check if point q lies on segment pr
template <typename N>
bool Earcut<N>::onSegment(const Node* p, const Node* q, const Node* r) {
return q->x <= std::max<double>(p->x, r->x) &&
q->x >= std::min<double>(p->x, r->x) &&
q->y <= std::max<double>(p->y, r->y) &&
q->y >= std::min<double>(p->y, r->y);
}
template <typename N>
int Earcut<N>::sign(double val) {
return (0.0 < val) - (val < 0.0);
}
// check if a polygon diagonal intersects any polygon segments
template <typename N>
bool Earcut<N>::intersectsPolygon(const Node* a, const Node* b) {
const Node* p = a;
do {
if (p->i != a->i && p->next->i != a->i && p->i != b->i && p->next->i != b->i &&
intersects(p, p->next, a, b)) return true;
p = p->next;
} while (p != a);
return false;
}
// check if a polygon diagonal is locally inside the polygon
template <typename N>
bool Earcut<N>::locallyInside(const Node* a, const Node* b) {
return area(a->prev, a, a->next) < 0 ?
area(a, b, a->next) >= 0 && area(a, a->prev, b) >= 0 :
area(a, b, a->prev) < 0 || area(a, a->next, b) < 0;
}
// check if the middle Vertex of a polygon diagonal is inside the polygon
template <typename N>
bool Earcut<N>::middleInside(const Node* a, const Node* b) {
const Node* p = a;
bool inside = false;
double px = (a->x + b->x) / 2;
double py = (a->y + b->y) / 2;
do {
if (((p->y > py) != (p->next->y > py)) && p->next->y != p->y &&
(px < (p->next->x - p->x) * (py - p->y) / (p->next->y - p->y) + p->x))
inside = !inside;
p = p->next;
} while (p != a);
return inside;
}
// link two polygon vertices with a bridge; if the vertices belong to the same ring, it splits
// polygon into two; if one belongs to the outer ring and another to a hole, it merges it into a
// single ring
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::splitPolygon(Node* a, Node* b) {
Node* a2 = nodes.construct(a->i, a->x, a->y);
Node* b2 = nodes.construct(b->i, b->x, b->y);
Node* an = a->next;
Node* bp = b->prev;
a->next = b;
b->prev = a;
a2->next = an;
an->prev = a2;
b2->next = a2;
a2->prev = b2;
bp->next = b2;
b2->prev = bp;
return b2;
}
// create a node and util::optionally link it with previous one (in a circular doubly linked list)
template <typename N> template <typename Point>
typename Earcut<N>::Node*
Earcut<N>::insertNode(std::size_t i, const Point& pt, Node* last) {
Node* p = nodes.construct(static_cast<N>(i), util::nth<0, Point>::get(pt), util::nth<1, Point>::get(pt));
if (!last) {
p->prev = p;
p->next = p;
} else {
assert(last);
p->next = last->next;
p->prev = last;
last->next->prev = p;
last->next = p;
}
return p;
}
template <typename N>
void Earcut<N>::removeNode(Node* p) {
p->next->prev = p->prev;
p->prev->next = p->next;
if (p->prevZ) p->prevZ->nextZ = p->nextZ;
if (p->nextZ) p->nextZ->prevZ = p->prevZ;
}
}
template <typename N = uint32_t, typename Polygon>
std::vector<N> earcut(const Polygon& poly) {
mapbox::detail::Earcut<N> earcut;
earcut(poly);
return std::move(earcut.indices);
}
}