Skip to content

Commit 8cc95fb

Browse files
AlePaluAlePalu
committed
stable support for polygon triangulation
1 parent 1b1a86e commit 8cc95fb

File tree

1 file changed

+115
-120
lines changed

1 file changed

+115
-120
lines changed

‎fdaPDE/geometry/polygon.h

Lines changed: 115 additions & 120 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,23 @@ template <int LocalDim, int EmbedDim> class Polygon {
3232

3333
// constructors
3434
Polygon() noexcept = default;
35-
Polygon(const DMatrix<double>& nodes) noexcept : triangulation_() { triangulate_(nodes); }
35+
Polygon(const DMatrix<double>& nodes) noexcept : triangulation_() {
36+
fdapde_assert(nodes.rows() > 0 && nodes.cols() == embed_dim);
37+
// check if nodes are sorted counter-clocwise
38+
double area = 0;
39+
int n_nodes = nodes.rows();
40+
for (int i = 0; i < n_nodes - 1; ++i) {
41+
area += (nodes(i, 0) + nodes(i + 1, 0)) * (nodes(i + 1, 1) - nodes(i, 1));
42+
}
43+
area += (nodes(n_nodes - 1, 0) + nodes(0, 0)) * (nodes(0, 1) - nodes(n_nodes - 1, 1)); // last edge
44+
if (area > 0) {
45+
triangulate_(nodes);
46+
} else { // nodes are in clocwise order, reverse node ordering
47+
DMatrix<double> reversed_nodes(n_nodes, embed_dim);
48+
for (int i = 0; i < n_nodes; ++i) { reversed_nodes.row(i) = nodes.row(n_nodes - 1 - i); }
49+
triangulate_(reversed_nodes);
50+
}
51+
}
3652
Polygon(const Polygon&) noexcept = default;
3753
Polygon(Polygon&&) noexcept = default;
3854
// observers
@@ -43,12 +59,6 @@ template <int LocalDim, int EmbedDim> class Polygon {
4359
int n_nodes() const { return triangulation_.n_nodes(); }
4460
int n_cells() const { return triangulation_.n_cells(); }
4561
int n_edges() const { return triangulation_.n_edges(); }
46-
// point location
47-
template <int Rows, int Cols>
48-
std::conditional_t<Rows == Dynamic || Cols == Dynamic, DVector<int>, int>
49-
locate(const Eigen::Matrix<double, Rows, Cols>& p) const {
50-
return triangulation_.locate(p);
51-
}
5262
// test whether point p is contained in polygon
5363
template <int Rows, int Cols>
5464
requires((Rows == embed_dim && Cols == 1) || (Cols == embed_dim && Rows == 1))
@@ -76,161 +86,147 @@ template <int LocalDim, int EmbedDim> class Polygon {
7686

7787
// partition an arbitrary polygon P into a set of monotone polygons (plane sweep approach, section 3.2 of De Berg,
7888
// M. (2000). Computational geometry: algorithms and applications. Springer Science & Business Media.)
79-
std::vector<std::vector<int>> monotone_partition_(const DMatrix<double>& nodes) {
89+
std::vector<std::vector<int>> monotone_partition_(const DMatrix<double>& coords) {
8090
using poly_t = DCEL<local_dim, embed_dim>;
81-
using node_t = typename poly_t::node_t;
82-
using edge_t = typename poly_t::edge_t;
83-
using node_ptr_t = std::add_pointer_t<node_t>;
84-
using edge_ptr_t = std::add_pointer_t<edge_t>;
91+
using halfedge_t = typename poly_t::halfedge_t;
92+
using halfedge_ptr_t = std::add_pointer_t<halfedge_t>;
8593
// data structure inducing an ad-hoc edge ordering for monotone partitioning
86-
struct edge_t_ {
94+
struct edge_t {
8795
double p1x, p1y; // p1 coordinates
8896
double p2x, p2y; // p2 coordinates
8997
int id;
9098
// constructor
91-
edge_t_(double p1x_, double p1y_, double p2x_, double p2y_) noexcept :
99+
edge_t(double p1x_, double p1y_, double p2x_, double p2y_) noexcept :
92100
p1x(p1x_), p1y(p1y_), p2x(p2x_), p2y(p2y_) { }
93-
edge_t_(int id_, double p1x_, double p1y_, double p2x_, double p2y_) noexcept :
101+
edge_t(int id_, double p1x_, double p1y_, double p2x_, double p2y_) noexcept :
94102
p1x(p1x_), p1y(p1y_), p2x(p2x_), p2y(p2y_), id(id_) { }
95103
// right-to-left edge ordering relation along x-coordinate
96-
bool operator<(const edge_t_& rhs) const { // for monotone partitioning, rhs is always below the p1 point
97-
if (rhs.p1y == rhs.p2y) { // rhs is horizontal
98-
if (p1y == p2y) { return (p1x < rhs.p1x); } // both edges are horizontal lines
104+
bool operator<(const edge_t& rhs) const { // for monotone partitioning, rhs is always below the p1 point
105+
if (rhs.p1y == rhs.p2y) { // rhs is horizontal
106+
if (p1y == p2y) { return (p1y < rhs.p1y); } // both edges are horizontal lines
99107
return internals::orientation(
100108
std::array {p2x, p2y}, std::array {p1x, p1y}, std::array {rhs.p1x, rhs.p1y}) ==
101109
internals::Orientation::LEFT;
102-
} else if (p1y <= p2y) {
110+
} else if (p1y == p2y || p1y < rhs.p1y) {
103111
return internals::orientation(
104112
std::array {rhs.p2x, rhs.p2y}, std::array {rhs.p1x, rhs.p1y}, std::array {p1x, p1y}) !=
105113
internals::Orientation::LEFT;
106114
} else {
107115
return internals::orientation(
108-
std::array {p2x, p2y}, std::array {rhs.p1x, rhs.p1y}, std::array {p1x, p1y}) ==
116+
std::array {p2x, p2y}, std::array {p1x, p1y}, std::array {rhs.p1x, rhs.p1y}) ==
109117
internals::Orientation::LEFT;
110118
}
111119
}
112120
};
121+
// given nodes p1, p2, asserts true if p1 is below p2 (induces a y-decresing ordering on polygon nodes)
122+
auto below = []<typename node_t>(const node_t& p1, const node_t& p2) {
123+
if (p1[1] < p2[1]) {
124+
return true;
125+
} else if (p1[1] == p2[1]) {
126+
if (p1[0] < p2[0]) { return true; }
127+
}
128+
return false;
129+
};
130+
113131
// O(n) polygon construction as Doubly Connected Edge List
114-
poly_t poly = DCEL<local_dim, embed_dim>::make_polygon(nodes);
115-
int n_nodes = poly.n_nodes();
116-
int n_edges = poly.n_edges();
117-
std::set<edge_t_> sweep_line; // edges pierced by sweep line, sorted by x-coord
118-
std::vector<node_ptr_t> helper(n_edges, nullptr);
119-
std::vector<node_ptr_t> node_ptr(n_nodes);
120-
for (auto it = std::make_pair(poly.nodes_begin(), 0); it.first != poly.nodes_end(); ++it.first, ++it.second) {
121-
node_ptr[it.second] = std::addressof(*it.first);
122-
}
123-
// given node id i, return the previous and next node ids on the counterclockwise sorted polygon border
124-
auto node_prev = [&](int node_id) { return node_id == 0 ? n_nodes - 1 : node_id - 1; };
125-
auto node_next = [&](int node_id) { return node_id == n_nodes - 1 ? 0 : node_id + 1; };
132+
poly_t dcel = DCEL<local_dim, embed_dim>::make_polygon(coords);
133+
int n_nodes = dcel.n_nodes();
134+
int n_edges = dcel.n_edges();
135+
std::set<edge_t> sweep_line; // edges pierced by sweep line, sorted by x-coord
136+
std::vector<halfedge_ptr_t> helper(n_nodes, nullptr);
137+
std::vector<halfedge_ptr_t> nodes(n_nodes); // maps node id to one of its halfedges
126138

127139
// O(n) node type detection
128140
enum node_category_t { start = 0, split = 1, end = 2, merge = 3, regular = 4 };
129-
std::vector<node_category_t> node_category(n_nodes);
130-
for (int i = 0; i < n_nodes; ++i) { // counterclocwise order loop, the interior of P is always on the left
131-
auto prev = nodes.row(node_prev(i));
132-
auto curr = nodes.row(i);
133-
auto next = nodes.row(node_next(i));
141+
std::unordered_map<halfedge_ptr_t, node_category_t> node_category;
142+
for (auto it = dcel.nodes_begin(); it != dcel.nodes_end(); ++it) { // counterclocwise order loop
143+
nodes[it->id()] = it->halfedge();
144+
auto prev = it->prev()->coords();
145+
auto curr = it->coords();
146+
auto next = it->next()->coords();
134147
double m_signed = internals::signed_measure_2d_tri(prev, curr, next);
135148
if (m_signed >= 0) { // interior angle less or equal than \pi
136-
if (prev[1] < curr[1] && next[1] < curr[1]) {
137-
node_category[i] = node_category_t::start;
138-
} else if (prev[1] > curr[1] && next[1] > curr[1]) {
139-
node_category[i] = node_category_t::end;
149+
if (below(prev, curr) && below(next, curr)) {
150+
node_category[nodes[it->id()]] = node_category_t::start;
151+
} else if (below(curr, prev) && below(curr, next)) {
152+
node_category[nodes[it->id()]] = node_category_t::end;
140153
} else {
141-
node_category[i] = node_category_t::regular;
154+
node_category[nodes[it->id()]] = node_category_t::regular;
142155
}
143-
} else if (m_signed < 0) { // interior angle grater than \pi
144-
if (prev[1] < curr[1] && next[1] < curr[1]) {
145-
node_category[i] = node_category_t::split;
146-
} else if (prev[1] > curr[1] && next[1] > curr[1]) {
147-
node_category[i] = node_category_t::merge;
156+
} else { // interior angle grater than \pi
157+
if (below(prev, curr) && below(next, curr)) {
158+
node_category[nodes[it->id()]] = node_category_t::split;
159+
} else if (below(curr, prev) && below(curr, next)) {
160+
node_category[nodes[it->id()]] = node_category_t::merge;
148161
} else {
149-
node_category[i] = node_category_t::regular;
162+
node_category[nodes[it->id()]] = node_category_t::regular;
150163
}
151164
}
152165
}
153166
// O(nlog(n)) y-coordinate sort (break tiles using x-coordinate)
154-
std::sort(node_ptr.begin(), node_ptr.end(), [&](node_ptr_t n, node_ptr_t m) {
155-
int i = n->id(), j = m->id();
156-
return nodes(i, 1) > nodes(j, 1) || (nodes(i, 1) == nodes(j, 1) && nodes(i, 0) < nodes(j, 0));
167+
std::sort(nodes.begin(), nodes.end(), [&](halfedge_ptr_t n, halfedge_ptr_t m) {
168+
int i = n->node()->id(), j = m->node()->id();
169+
return coords(i, 1) > coords(j, 1) || (coords(i, 1) == coords(j, 1) && coords(i, 0) > coords(j, 0));
157170
});
158171
// monotone partitioning algorithm (details in section 3.2 of De Berg, M. (2000). Computational geometry:
159172
// algorithms and applications. Springer Science & Business Media.)
160-
std::vector<edge_ptr_t> diagonals;
161-
std::vector<typename std::set<edge_t_>::iterator> sweep_line_it(n_edges);
162-
sweep_line_it.resize(n_edges);
163-
// given edge {p1, p2}, asserts true if the polygon interior is on the left of {p1, p2}
164-
auto interior_on_left = []<typename point_t>(const point_t& p1, const point_t& p2) {
165-
if (p1[1] < p2[1]) {
166-
return true;
167-
} else if (p1[1] == p2[1]) {
168-
if (p1[0] < p2[0]) { return true; }
169-
}
170-
return false;
171-
};
172-
173-
for (node_ptr_t v : node_ptr) { // loops in decreasing y-coordinate order
174-
int v_i = v->id();
175-
// process i-th node
176-
switch (node_category[v_i]) {
173+
std::vector<typename std::set<edge_t>::iterator> sweep_line_it(n_nodes);
174+
sweep_line_it.resize(n_nodes);
175+
for (halfedge_ptr_t v : nodes) { // loops in decreasing y-coordinate order
176+
int prev = v->node()->prev()->id();
177+
int curr = v->node()->id();
178+
int next = v->node()->next()->id();
179+
// process i-th node
180+
switch (node_category[v]) {
177181
case node_category_t::start: {
178-
auto ref = sweep_line.emplace(
179-
v_i, nodes(v_i, 0), nodes(v_i, 1), nodes(node_next(v_i), 0), nodes(node_next(v_i), 1));
180-
sweep_line_it[v_i] = ref.first;
181-
helper[v_i] = v;
182+
auto ref = sweep_line.emplace(curr, coords(curr, 0), coords(curr, 1), coords(next, 0), coords(next, 1));
183+
sweep_line_it[curr] = ref.first;
184+
helper[curr] = v;
182185
break;
183186
}
184187
case node_category_t::end: {
185-
if (node_category[helper[node_prev(v_i)]->id()] == node_category_t::merge) {
186-
diagonals.push_back(poly.insert_edge(v, helper[node_prev(v_i)]));
187-
}
188-
sweep_line.erase(sweep_line_it[node_prev(v_i)]);
188+
if (node_category[helper[prev]] == node_category_t::merge) { dcel.insert_edge(v, helper[prev]); }
189+
sweep_line.erase(sweep_line_it[prev]);
189190
break;
190191
}
191192
case node_category_t::split: {
192-
// search edge in sweep line directly left to v_i
193-
edge_t_ edge(nodes(v_i, 0), nodes(v_i, 1), nodes(v_i, 0), nodes(v_i, 1));
193+
// search edge in sweep line directly left to curr
194+
edge_t edge(coords(curr, 0), coords(curr, 1), coords(curr, 0), coords(curr, 1));
194195
auto it = sweep_line.lower_bound(edge);
195-
// insert diagonal
196-
diagonals.push_back(poly.insert_edge(v, helper[it->id]));
197-
helper[it->id] = v;
198-
auto ref = sweep_line.emplace(
199-
v_i, nodes(v_i, 0), nodes(v_i, 1), nodes(node_next(v_i), 0), nodes(node_next(v_i), 1));
200-
sweep_line_it[v_i] = ref.first;
201-
helper[v_i] = v;
196+
// insert diagonal (update v to point to the inserted diagonal)
197+
helper[it->id] = dcel.insert_edge(v, helper[it->id]);
198+
auto ref = sweep_line.emplace(curr, coords(curr, 0), coords(curr, 1), coords(next, 0), coords(next, 1));
199+
sweep_line_it[curr] = ref.first;
200+
helper[curr] = v;
202201
break;
203202
}
204203
case node_category_t::merge: {
205-
if (node_category[helper[node_prev(v_i)]->id()] == node_category_t::merge) { // insert diagonal
206-
diagonals.push_back(poly.insert_edge(v, helper[node_prev(v_i)]));
207-
}
208-
sweep_line.erase(sweep_line_it[node_prev(v_i)]);
209-
// search edge in sweep line directly left to v_i
210-
edge_t_ edge(nodes(v_i, 0), nodes(v_i, 1), nodes(v_i, 0), nodes(v_i, 1));
204+
if (node_category[helper[prev]] == node_category_t::merge) { dcel.insert_edge(v, helper[prev]); }
205+
sweep_line.erase(sweep_line_it[prev]);
206+
// search edge in sweep line directly left to curr
207+
edge_t edge(coords(curr, 0), coords(curr, 1), coords(curr, 0), coords(curr, 1));
211208
auto it = sweep_line.lower_bound(edge);
212-
if (node_category[it->id] == node_category_t::merge) {
213-
diagonals.push_back(poly.insert_edge(v, helper[it->id]));
209+
if (node_category[helper[it->id]] == node_category_t::merge) {
210+
v = dcel.insert_edge(v, helper[it->id]); // update v to point to the inserted diagonal
214211
}
215212
helper[it->id] = v;
216213
break;
217214
}
218215
case node_category_t::regular: {
219-
if (!interior_on_left(nodes.row(v_i), nodes.row(node_next(v_i)))) {
220-
if (node_category[helper[node_prev(v_i)]->id()] == node_category_t::merge) { // insert diagonal
221-
diagonals.push_back(poly.insert_edge(v, helper[node_prev(v_i)]));
222-
}
223-
sweep_line.erase(sweep_line_it[node_prev(v_i)]);
224-
auto ref = sweep_line.emplace(
225-
v_i, nodes(v_i, 0), nodes(v_i, 1), nodes(node_next(v_i), 0), nodes(node_next(v_i), 1));
226-
sweep_line_it[v_i] = ref.first;
227-
helper[v_i] = v;
216+
if (!below(coords.row(curr), coords.row(next))) {
217+
// polygon interior is on the right of this halfedge
218+
if (node_category[helper[prev]] == node_category_t::merge) { dcel.insert_edge(v, helper[prev]); }
219+
sweep_line.erase(sweep_line_it[prev]);
220+
auto ref =
221+
sweep_line.emplace(curr, coords(curr, 0), coords(curr, 1), coords(next, 0), coords(next, 1));
222+
sweep_line_it[curr] = ref.first;
223+
helper[curr] = v;
228224
} else {
229-
// search edge in sweep line directly left to v_i
230-
edge_t_ edge(nodes(v_i, 0), nodes(v_i, 1), nodes(v_i, 0), nodes(v_i, 1));
225+
// search edge in sweep line directly left to curr
226+
edge_t edge(coords(curr, 0), coords(curr, 1), coords(curr, 0), coords(curr, 1));
231227
auto it = sweep_line.lower_bound(edge);
232-
if (node_category[helper[it->id]->id()] == node_category_t::merge) {
233-
diagonals.push_back(poly.insert_edge(v, helper[it->id]));
228+
if (node_category[helper[it->id]] == node_category_t::merge) {
229+
v = dcel.insert_edge(v, helper[it->id]);
234230
}
235231
helper[it->id] = v;
236232
}
@@ -239,24 +235,24 @@ template <int LocalDim, int EmbedDim> class Polygon {
239235
}
240236
}
241237
// recover from the DCEL structure the node numbering of each monotone polygon
242-
std::vector<bool> visited(poly.n_halfedges(), false);
238+
std::vector<bool> visited(dcel.n_halfedges(), false);
243239
std::vector<std::vector<int>> monotone_partition;
244-
auto collect_polygon_node_id_from = [&]<typename halfedge_t>(halfedge_t* chain) mutable {
245-
auto& it = monotone_partition.emplace_back();
246-
for (auto jt = poly.halfedge_begin(chain); jt != poly.halfedge_end(chain); ++jt) {
247-
int node_id = jt->node()->id();
248-
it.push_back(node_id);
249-
visited[jt->id()] = true;
240+
if (dcel.n_cells() == 1) { // polygon was already monotone
241+
auto& it = monotone_partition.emplace_back();
242+
it.resize(n_nodes);
243+
std::iota(it.begin(), it.end(), 0);
244+
} else {
245+
for (auto cell = dcel.cells_begin(); cell != dcel.cells_end(); ++cell) {
246+
auto& it = monotone_partition.emplace_back();
247+
halfedge_t *chain = cell->halfedge(), *end = chain;
248+
do {
249+
it.push_back(chain->node()->id());
250+
chain = chain->next();
251+
} while (end != chain);
250252
}
251-
};
252-
for (edge_ptr_t diagonal : diagonals) {
253-
// follow the chains from the halfedges composing this edge, if halfedge not already visited
254-
if (!visited[diagonal->first() ->id()]) { collect_polygon_node_id_from(diagonal->first ()); }
255-
if (!visited[diagonal->second()->id()]) { collect_polygon_node_id_from(diagonal->second()); }
256253
}
257254
return monotone_partition;
258255
}
259-
260256
// triangulate monotone polygon (returns a RowMajor ordered matrix of cells).
261257
std::vector<int> triangulate_monotone_(const DMatrix<double>& nodes) {
262258
// every triangulation of a polygon of n points has n - 2 triangles (lemma 1.2.2 of (1))
@@ -317,7 +313,6 @@ template <int LocalDim, int EmbedDim> class Polygon {
317313
for (std::size_t j = 2, n = node.size(); j < n; ++j) {
318314
node_i = *(reflex_chain.end() - 1);
319315
node_j = node[j];
320-
321316
if (are_in_opposite_chains(node_i, node_j)) { // triangulate
322317
node_i = *reflex_chain.begin();
323318
on_left = !on_left;

0 commit comments

Comments
 (0)