1 |
5 |
maiden |
/*
|
2 |
|
|
* Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors
|
3 |
|
|
* http://code.google.com/p/poly2tri/
|
4 |
|
|
*
|
5 |
|
|
* All rights reserved.
|
6 |
|
|
*
|
7 |
|
|
* Redistribution and use in source and binary forms, with or without modification,
|
8 |
|
|
* are permitted provided that the following conditions are met:
|
9 |
|
|
*
|
10 |
|
|
* * Redistributions of source code must retain the above copyright notice,
|
11 |
|
|
* this list of conditions and the following disclaimer.
|
12 |
|
|
* * Redistributions in binary form must reproduce the above copyright notice,
|
13 |
|
|
* this list of conditions and the following disclaimer in the documentation
|
14 |
|
|
* and/or other materials provided with the distribution.
|
15 |
|
|
* * Neither the name of Poly2Tri nor the names of its contributors may be
|
16 |
|
|
* used to endorse or promote products derived from this software without specific
|
17 |
|
|
* prior written permission.
|
18 |
|
|
*
|
19 |
|
|
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
20 |
|
|
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
21 |
|
|
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
22 |
|
|
* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
|
23 |
|
|
* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
24 |
|
|
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
25 |
|
|
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
26 |
|
|
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
27 |
|
|
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
28 |
|
|
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
29 |
|
|
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
30 |
|
|
*/
|
31 |
|
|
#include <stdexcept>
|
32 |
|
|
#include "sweep.h"
|
33 |
|
|
#include "sweep_context.h"
|
34 |
|
|
#include "advancing_front.h"
|
35 |
|
|
#include "../common/utils.h"
|
36 |
|
|
|
37 |
|
|
namespace p2t {
|
38 |
|
|
|
39 |
|
|
// Triangulate simple polygon with holes
|
40 |
|
|
void Sweep::Triangulate(SweepContext& tcx)
|
41 |
|
|
{
|
42 |
|
|
tcx.InitTriangulation();
|
43 |
|
|
tcx.CreateAdvancingFront(nodes_);
|
44 |
|
|
// Sweep points; build mesh
|
45 |
|
|
SweepPoints(tcx);
|
46 |
|
|
// Clean up
|
47 |
|
|
FinalizationPolygon(tcx);
|
48 |
|
|
}
|
49 |
|
|
|
50 |
|
|
void Sweep::SweepPoints(SweepContext& tcx)
|
51 |
|
|
{
|
52 |
|
|
for (int i = 1; i < tcx.point_count(); i++) {
|
53 |
|
|
Point& point = *tcx.GetPoint(i);
|
54 |
|
|
Node* node = &PointEvent(tcx, point);
|
55 |
|
|
for (unsigned int i = 0; i < point.edge_list.size(); i++) {
|
56 |
|
|
EdgeEvent(tcx, point.edge_list[i], node);
|
57 |
|
|
}
|
58 |
|
|
}
|
59 |
|
|
}
|
60 |
|
|
|
61 |
|
|
void Sweep::FinalizationPolygon(SweepContext& tcx)
|
62 |
|
|
{
|
63 |
|
|
// Get an Internal triangle to start with
|
64 |
|
|
Triangle* t = tcx.front()->head()->next->triangle;
|
65 |
|
|
Point* p = tcx.front()->head()->next->point;
|
66 |
|
|
while (!t->GetConstrainedEdgeCW(*p)) {
|
67 |
|
|
t = t->NeighborCCW(*p);
|
68 |
|
|
}
|
69 |
|
|
|
70 |
|
|
// Collect interior triangles constrained by edges
|
71 |
|
|
tcx.MeshClean(*t);
|
72 |
|
|
}
|
73 |
|
|
|
74 |
|
|
Node& Sweep::PointEvent(SweepContext& tcx, Point& point)
|
75 |
|
|
{
|
76 |
|
|
Node& node = tcx.LocateNode(point);
|
77 |
|
|
Node& new_node = NewFrontTriangle(tcx, point, node);
|
78 |
|
|
|
79 |
|
|
// Only need to check +epsilon since point never have smaller
|
80 |
|
|
// x value than node due to how we fetch nodes from the front
|
81 |
|
|
if (point.x <= node.point->x + EPSILON) {
|
82 |
|
|
Fill(tcx, node);
|
83 |
|
|
}
|
84 |
|
|
|
85 |
|
|
//tcx.AddNode(new_node);
|
86 |
|
|
|
87 |
|
|
FillAdvancingFront(tcx, new_node);
|
88 |
|
|
return new_node;
|
89 |
|
|
}
|
90 |
|
|
|
91 |
|
|
void Sweep::EdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
|
92 |
|
|
{
|
93 |
|
|
tcx.edge_event.constrained_edge = edge;
|
94 |
|
|
tcx.edge_event.right = (edge->p->x > edge->q->x);
|
95 |
|
|
|
96 |
|
|
if (IsEdgeSideOfTriangle(*node->triangle, *edge->p, *edge->q)) {
|
97 |
|
|
return;
|
98 |
|
|
}
|
99 |
|
|
|
100 |
|
|
// For now we will do all needed filling
|
101 |
|
|
// TODO: integrate with flip process might give some better performance
|
102 |
|
|
// but for now this avoid the issue with cases that needs both flips and fills
|
103 |
|
|
FillEdgeEvent(tcx, edge, node);
|
104 |
|
|
EdgeEvent(tcx, *edge->p, *edge->q, node->triangle, *edge->q);
|
105 |
|
|
}
|
106 |
|
|
|
107 |
|
|
void Sweep::EdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point)
|
108 |
|
|
{
|
109 |
|
|
if (IsEdgeSideOfTriangle(*triangle, ep, eq)) {
|
110 |
|
|
return;
|
111 |
|
|
}
|
112 |
|
|
|
113 |
|
|
Point* p1 = triangle->PointCCW(point);
|
114 |
|
|
Orientation o1 = Orient2d(eq, *p1, ep);
|
115 |
|
|
if (o1 == COLLINEAR) {
|
116 |
|
|
if( triangle->Contains(&eq, p1)) {
|
117 |
|
|
triangle->MarkConstrainedEdge(&eq, p1 );
|
118 |
|
|
// We are modifying the constraint maybe it would be better to
|
119 |
|
|
// not change the given constraint and just keep a variable for the new constraint
|
120 |
|
|
tcx.edge_event.constrained_edge->q = p1;
|
121 |
|
|
triangle = &triangle->NeighborAcross(point);
|
122 |
|
|
EdgeEvent( tcx, ep, *p1, triangle, *p1 );
|
123 |
|
|
} else {
|
124 |
|
|
std::runtime_error("EdgeEvent - collinear points not supported");
|
125 |
|
|
assert(0);
|
126 |
|
|
}
|
127 |
|
|
return;
|
128 |
|
|
}
|
129 |
|
|
|
130 |
|
|
Point* p2 = triangle->PointCW(point);
|
131 |
|
|
Orientation o2 = Orient2d(eq, *p2, ep);
|
132 |
|
|
if (o2 == COLLINEAR) {
|
133 |
|
|
if( triangle->Contains(&eq, p2)) {
|
134 |
|
|
triangle->MarkConstrainedEdge(&eq, p2 );
|
135 |
|
|
// We are modifying the constraint maybe it would be better to
|
136 |
|
|
// not change the given constraint and just keep a variable for the new constraint
|
137 |
|
|
tcx.edge_event.constrained_edge->q = p2;
|
138 |
|
|
triangle = &triangle->NeighborAcross(point);
|
139 |
|
|
EdgeEvent( tcx, ep, *p2, triangle, *p2 );
|
140 |
|
|
} else {
|
141 |
|
|
std::runtime_error("EdgeEvent - collinear points not supported");
|
142 |
|
|
assert(0);
|
143 |
|
|
}
|
144 |
|
|
return;
|
145 |
|
|
}
|
146 |
|
|
|
147 |
|
|
if (o1 == o2) {
|
148 |
|
|
// Need to decide if we are rotating CW or CCW to get to a triangle
|
149 |
|
|
// that will cross edge
|
150 |
|
|
if (o1 == CW) {
|
151 |
|
|
triangle = triangle->NeighborCCW(point);
|
152 |
|
|
} else{
|
153 |
|
|
triangle = triangle->NeighborCW(point);
|
154 |
|
|
}
|
155 |
|
|
EdgeEvent(tcx, ep, eq, triangle, point);
|
156 |
|
|
} else {
|
157 |
|
|
// This triangle crosses constraint so lets flippin start!
|
158 |
|
|
FlipEdgeEvent(tcx, ep, eq, triangle, point);
|
159 |
|
|
}
|
160 |
|
|
}
|
161 |
|
|
|
162 |
|
|
bool Sweep::IsEdgeSideOfTriangle(Triangle& triangle, Point& ep, Point& eq)
|
163 |
|
|
{
|
164 |
|
|
int index = triangle.EdgeIndex(&ep, &eq);
|
165 |
|
|
|
166 |
|
|
if (index != -1) {
|
167 |
|
|
triangle.MarkConstrainedEdge(index);
|
168 |
|
|
Triangle* t = triangle.GetNeighbor(index);
|
169 |
|
|
if (t) {
|
170 |
|
|
t->MarkConstrainedEdge(&ep, &eq);
|
171 |
|
|
}
|
172 |
|
|
return true;
|
173 |
|
|
}
|
174 |
|
|
return false;
|
175 |
|
|
}
|
176 |
|
|
|
177 |
|
|
Node& Sweep::NewFrontTriangle(SweepContext& tcx, Point& point, Node& node)
|
178 |
|
|
{
|
179 |
|
|
Triangle* triangle = new Triangle(point, *node.point, *node.next->point);
|
180 |
|
|
|
181 |
|
|
triangle->MarkNeighbor(*node.triangle);
|
182 |
|
|
tcx.AddToMap(triangle);
|
183 |
|
|
|
184 |
|
|
Node* new_node = new Node(point);
|
185 |
|
|
nodes_.push_back(new_node);
|
186 |
|
|
|
187 |
|
|
new_node->next = node.next;
|
188 |
|
|
new_node->prev = &node;
|
189 |
|
|
node.next->prev = new_node;
|
190 |
|
|
node.next = new_node;
|
191 |
|
|
|
192 |
|
|
if (!Legalize(tcx, *triangle)) {
|
193 |
|
|
tcx.MapTriangleToNodes(*triangle);
|
194 |
|
|
}
|
195 |
|
|
|
196 |
|
|
return *new_node;
|
197 |
|
|
}
|
198 |
|
|
|
199 |
|
|
void Sweep::Fill(SweepContext& tcx, Node& node)
|
200 |
|
|
{
|
201 |
|
|
Triangle* triangle = new Triangle(*node.prev->point, *node.point, *node.next->point);
|
202 |
|
|
|
203 |
|
|
// TODO: should copy the constrained_edge value from neighbor triangles
|
204 |
|
|
// for now constrained_edge values are copied during the legalize
|
205 |
|
|
triangle->MarkNeighbor(*node.prev->triangle);
|
206 |
|
|
triangle->MarkNeighbor(*node.triangle);
|
207 |
|
|
|
208 |
|
|
tcx.AddToMap(triangle);
|
209 |
|
|
|
210 |
|
|
// Update the advancing front
|
211 |
|
|
node.prev->next = node.next;
|
212 |
|
|
node.next->prev = node.prev;
|
213 |
|
|
|
214 |
|
|
// If it was legalized the triangle has already been mapped
|
215 |
|
|
if (!Legalize(tcx, *triangle)) {
|
216 |
|
|
tcx.MapTriangleToNodes(*triangle);
|
217 |
|
|
}
|
218 |
|
|
|
219 |
|
|
}
|
220 |
|
|
|
221 |
|
|
void Sweep::FillAdvancingFront(SweepContext& tcx, Node& n)
|
222 |
|
|
{
|
223 |
|
|
|
224 |
|
|
// Fill right holes
|
225 |
|
|
Node* node = n.next;
|
226 |
|
|
|
227 |
|
|
while (node->next) {
|
228 |
|
|
// if HoleAngle exceeds 90 degrees then break.
|
229 |
|
|
if (LargeHole_DontFill(node)) break;
|
230 |
|
|
Fill(tcx, *node);
|
231 |
|
|
node = node->next;
|
232 |
|
|
}
|
233 |
|
|
|
234 |
|
|
// Fill left holes
|
235 |
|
|
node = n.prev;
|
236 |
|
|
|
237 |
|
|
while (node->prev) {
|
238 |
|
|
// if HoleAngle exceeds 90 degrees then break.
|
239 |
|
|
if (LargeHole_DontFill(node)) break;
|
240 |
|
|
Fill(tcx, *node);
|
241 |
|
|
node = node->prev;
|
242 |
|
|
}
|
243 |
|
|
|
244 |
|
|
// Fill right basins
|
245 |
|
|
if (n.next && n.next->next) {
|
246 |
|
|
double angle = BasinAngle(n);
|
247 |
|
|
if (angle < PI_3div4) {
|
248 |
|
|
FillBasin(tcx, n);
|
249 |
|
|
}
|
250 |
|
|
}
|
251 |
|
|
}
|
252 |
|
|
|
253 |
|
|
// True if HoleAngle exceeds 90 degrees.
|
254 |
|
|
bool Sweep::LargeHole_DontFill(Node* node) {
|
255 |
|
|
|
256 |
|
|
Node* nextNode = node->next;
|
257 |
|
|
Node* prevNode = node->prev;
|
258 |
|
|
if (!AngleExceeds90Degrees(node->point, nextNode->point, prevNode->point))
|
259 |
|
|
return false;
|
260 |
|
|
|
261 |
|
|
// Check additional points on front.
|
262 |
|
|
Node* next2Node = nextNode->next;
|
263 |
|
|
// "..Plus.." because only want angles on same side as point being added.
|
264 |
|
|
if ((next2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, next2Node->point, prevNode->point))
|
265 |
|
|
return false;
|
266 |
|
|
|
267 |
|
|
Node* prev2Node = prevNode->prev;
|
268 |
|
|
// "..Plus.." because only want angles on same side as point being added.
|
269 |
|
|
if ((prev2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, nextNode->point, prev2Node->point))
|
270 |
|
|
return false;
|
271 |
|
|
|
272 |
|
|
return true;
|
273 |
|
|
}
|
274 |
|
|
|
275 |
|
|
bool Sweep::AngleExceeds90Degrees(Point* origin, Point* pa, Point* pb) {
|
276 |
|
|
double angle = Angle(*origin, *pa, *pb);
|
277 |
|
|
bool exceeds90Degrees = ((angle > PI_div2) || (angle < -PI_div2));
|
278 |
|
|
return exceeds90Degrees;
|
279 |
|
|
}
|
280 |
|
|
|
281 |
|
|
bool Sweep::AngleExceedsPlus90DegreesOrIsNegative(Point* origin, Point* pa, Point* pb) {
|
282 |
|
|
double angle = Angle(*origin, *pa, *pb);
|
283 |
|
|
bool exceedsPlus90DegreesOrIsNegative = (angle > PI_div2) || (angle < 0);
|
284 |
|
|
return exceedsPlus90DegreesOrIsNegative;
|
285 |
|
|
}
|
286 |
|
|
|
287 |
|
|
double Sweep::Angle(Point& origin, Point& pa, Point& pb) {
|
288 |
|
|
/* Complex plane
|
289 |
|
|
* ab = cosA +i*sinA
|
290 |
|
|
* ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
|
291 |
|
|
* atan2(y,x) computes the principal value of the argument function
|
292 |
|
|
* applied to the complex number x+iy
|
293 |
|
|
* Where x = ax*bx + ay*by
|
294 |
|
|
* y = ax*by - ay*bx
|
295 |
|
|
*/
|
296 |
|
|
double px = origin.x;
|
297 |
|
|
double py = origin.y;
|
298 |
|
|
double ax = pa.x- px;
|
299 |
|
|
double ay = pa.y - py;
|
300 |
|
|
double bx = pb.x - px;
|
301 |
|
|
double by = pb.y - py;
|
302 |
|
|
double x = ax * by - ay * bx;
|
303 |
|
|
double y = ax * bx + ay * by;
|
304 |
|
|
double angle = atan2(x, y);
|
305 |
|
|
return angle;
|
306 |
|
|
}
|
307 |
|
|
|
308 |
|
|
double Sweep::BasinAngle(Node& node)
|
309 |
|
|
{
|
310 |
|
|
double ax = node.point->x - node.next->next->point->x;
|
311 |
|
|
double ay = node.point->y - node.next->next->point->y;
|
312 |
|
|
return atan2(ay, ax);
|
313 |
|
|
}
|
314 |
|
|
|
315 |
|
|
double Sweep::HoleAngle(Node& node)
|
316 |
|
|
{
|
317 |
|
|
/* Complex plane
|
318 |
|
|
* ab = cosA +i*sinA
|
319 |
|
|
* ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
|
320 |
|
|
* atan2(y,x) computes the principal value of the argument function
|
321 |
|
|
* applied to the complex number x+iy
|
322 |
|
|
* Where x = ax*bx + ay*by
|
323 |
|
|
* y = ax*by - ay*bx
|
324 |
|
|
*/
|
325 |
|
|
double ax = node.next->point->x - node.point->x;
|
326 |
|
|
double ay = node.next->point->y - node.point->y;
|
327 |
|
|
double bx = node.prev->point->x - node.point->x;
|
328 |
|
|
double by = node.prev->point->y - node.point->y;
|
329 |
|
|
return atan2(ax * by - ay * bx, ax * bx + ay * by);
|
330 |
|
|
}
|
331 |
|
|
|
332 |
|
|
bool Sweep::Legalize(SweepContext& tcx, Triangle& t)
|
333 |
|
|
{
|
334 |
|
|
// To legalize a triangle we start by finding if any of the three edges
|
335 |
|
|
// violate the Delaunay condition
|
336 |
|
|
for (int i = 0; i < 3; i++) {
|
337 |
|
|
if (t.delaunay_edge[i])
|
338 |
|
|
continue;
|
339 |
|
|
|
340 |
|
|
Triangle* ot = t.GetNeighbor(i);
|
341 |
|
|
|
342 |
|
|
if (ot) {
|
343 |
|
|
Point* p = t.GetPoint(i);
|
344 |
|
|
Point* op = ot->OppositePoint(t, *p);
|
345 |
|
|
int oi = ot->Index(op);
|
346 |
|
|
|
347 |
|
|
// If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization)
|
348 |
|
|
// then we should not try to legalize
|
349 |
|
|
if (ot->constrained_edge[oi] || ot->delaunay_edge[oi]) {
|
350 |
|
|
t.constrained_edge[i] = ot->constrained_edge[oi];
|
351 |
|
|
continue;
|
352 |
|
|
}
|
353 |
|
|
|
354 |
|
|
bool inside = Incircle(*p, *t.PointCCW(*p), *t.PointCW(*p), *op);
|
355 |
|
|
|
356 |
|
|
if (inside) {
|
357 |
|
|
// Lets mark this shared edge as Delaunay
|
358 |
|
|
t.delaunay_edge[i] = true;
|
359 |
|
|
ot->delaunay_edge[oi] = true;
|
360 |
|
|
|
361 |
|
|
// Lets rotate shared edge one vertex CW to legalize it
|
362 |
|
|
RotateTrianglePair(t, *p, *ot, *op);
|
363 |
|
|
|
364 |
|
|
// We now got one valid Delaunay Edge shared by two triangles
|
365 |
|
|
// This gives us 4 new edges to check for Delaunay
|
366 |
|
|
|
367 |
|
|
// Make sure that triangle to node mapping is done only one time for a specific triangle
|
368 |
|
|
bool not_legalized = !Legalize(tcx, t);
|
369 |
|
|
if (not_legalized) {
|
370 |
|
|
tcx.MapTriangleToNodes(t);
|
371 |
|
|
}
|
372 |
|
|
|
373 |
|
|
not_legalized = !Legalize(tcx, *ot);
|
374 |
|
|
if (not_legalized)
|
375 |
|
|
tcx.MapTriangleToNodes(*ot);
|
376 |
|
|
|
377 |
|
|
// Reset the Delaunay edges, since they only are valid Delaunay edges
|
378 |
|
|
// until we add a new triangle or point.
|
379 |
|
|
// XXX: need to think about this. Can these edges be tried after we
|
380 |
|
|
// return to previous recursive level?
|
381 |
|
|
t.delaunay_edge[i] = false;
|
382 |
|
|
ot->delaunay_edge[oi] = false;
|
383 |
|
|
|
384 |
|
|
// If triangle have been legalized no need to check the other edges since
|
385 |
|
|
// the recursive legalization will handles those so we can end here.
|
386 |
|
|
return true;
|
387 |
|
|
}
|
388 |
|
|
}
|
389 |
|
|
}
|
390 |
|
|
return false;
|
391 |
|
|
}
|
392 |
|
|
|
393 |
|
|
bool Sweep::Incircle(Point& pa, Point& pb, Point& pc, Point& pd)
|
394 |
|
|
{
|
395 |
|
|
double adx = pa.x - pd.x;
|
396 |
|
|
double ady = pa.y - pd.y;
|
397 |
|
|
double bdx = pb.x - pd.x;
|
398 |
|
|
double bdy = pb.y - pd.y;
|
399 |
|
|
|
400 |
|
|
double adxbdy = adx * bdy;
|
401 |
|
|
double bdxady = bdx * ady;
|
402 |
|
|
double oabd = adxbdy - bdxady;
|
403 |
|
|
|
404 |
|
|
if (oabd <= 0)
|
405 |
|
|
return false;
|
406 |
|
|
|
407 |
|
|
double cdx = pc.x - pd.x;
|
408 |
|
|
double cdy = pc.y - pd.y;
|
409 |
|
|
|
410 |
|
|
double cdxady = cdx * ady;
|
411 |
|
|
double adxcdy = adx * cdy;
|
412 |
|
|
double ocad = cdxady - adxcdy;
|
413 |
|
|
|
414 |
|
|
if (ocad <= 0)
|
415 |
|
|
return false;
|
416 |
|
|
|
417 |
|
|
double bdxcdy = bdx * cdy;
|
418 |
|
|
double cdxbdy = cdx * bdy;
|
419 |
|
|
|
420 |
|
|
double alift = adx * adx + ady * ady;
|
421 |
|
|
double blift = bdx * bdx + bdy * bdy;
|
422 |
|
|
double clift = cdx * cdx + cdy * cdy;
|
423 |
|
|
|
424 |
|
|
double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd;
|
425 |
|
|
|
426 |
|
|
return det > 0;
|
427 |
|
|
}
|
428 |
|
|
|
429 |
|
|
void Sweep::RotateTrianglePair(Triangle& t, Point& p, Triangle& ot, Point& op)
|
430 |
|
|
{
|
431 |
|
|
Triangle* n1, *n2, *n3, *n4;
|
432 |
|
|
n1 = t.NeighborCCW(p);
|
433 |
|
|
n2 = t.NeighborCW(p);
|
434 |
|
|
n3 = ot.NeighborCCW(op);
|
435 |
|
|
n4 = ot.NeighborCW(op);
|
436 |
|
|
|
437 |
|
|
bool ce1, ce2, ce3, ce4;
|
438 |
|
|
ce1 = t.GetConstrainedEdgeCCW(p);
|
439 |
|
|
ce2 = t.GetConstrainedEdgeCW(p);
|
440 |
|
|
ce3 = ot.GetConstrainedEdgeCCW(op);
|
441 |
|
|
ce4 = ot.GetConstrainedEdgeCW(op);
|
442 |
|
|
|
443 |
|
|
bool de1, de2, de3, de4;
|
444 |
|
|
de1 = t.GetDelunayEdgeCCW(p);
|
445 |
|
|
de2 = t.GetDelunayEdgeCW(p);
|
446 |
|
|
de3 = ot.GetDelunayEdgeCCW(op);
|
447 |
|
|
de4 = ot.GetDelunayEdgeCW(op);
|
448 |
|
|
|
449 |
|
|
t.Legalize(p, op);
|
450 |
|
|
ot.Legalize(op, p);
|
451 |
|
|
|
452 |
|
|
// Remap delaunay_edge
|
453 |
|
|
ot.SetDelunayEdgeCCW(p, de1);
|
454 |
|
|
t.SetDelunayEdgeCW(p, de2);
|
455 |
|
|
t.SetDelunayEdgeCCW(op, de3);
|
456 |
|
|
ot.SetDelunayEdgeCW(op, de4);
|
457 |
|
|
|
458 |
|
|
// Remap constrained_edge
|
459 |
|
|
ot.SetConstrainedEdgeCCW(p, ce1);
|
460 |
|
|
t.SetConstrainedEdgeCW(p, ce2);
|
461 |
|
|
t.SetConstrainedEdgeCCW(op, ce3);
|
462 |
|
|
ot.SetConstrainedEdgeCW(op, ce4);
|
463 |
|
|
|
464 |
|
|
// Remap neighbors
|
465 |
|
|
// XXX: might optimize the markNeighbor by keeping track of
|
466 |
|
|
// what side should be assigned to what neighbor after the
|
467 |
|
|
// rotation. Now mark neighbor does lots of testing to find
|
468 |
|
|
// the right side.
|
469 |
|
|
t.ClearNeighbors();
|
470 |
|
|
ot.ClearNeighbors();
|
471 |
|
|
if (n1) ot.MarkNeighbor(*n1);
|
472 |
|
|
if (n2) t.MarkNeighbor(*n2);
|
473 |
|
|
if (n3) t.MarkNeighbor(*n3);
|
474 |
|
|
if (n4) ot.MarkNeighbor(*n4);
|
475 |
|
|
t.MarkNeighbor(ot);
|
476 |
|
|
}
|
477 |
|
|
|
478 |
|
|
void Sweep::FillBasin(SweepContext& tcx, Node& node)
|
479 |
|
|
{
|
480 |
|
|
if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
|
481 |
|
|
tcx.basin.left_node = node.next->next;
|
482 |
|
|
} else {
|
483 |
|
|
tcx.basin.left_node = node.next;
|
484 |
|
|
}
|
485 |
|
|
|
486 |
|
|
// Find the bottom and right node
|
487 |
|
|
tcx.basin.bottom_node = tcx.basin.left_node;
|
488 |
|
|
while (tcx.basin.bottom_node->next
|
489 |
|
|
&& tcx.basin.bottom_node->point->y >= tcx.basin.bottom_node->next->point->y) {
|
490 |
|
|
tcx.basin.bottom_node = tcx.basin.bottom_node->next;
|
491 |
|
|
}
|
492 |
|
|
if (tcx.basin.bottom_node == tcx.basin.left_node) {
|
493 |
|
|
// No valid basin
|
494 |
|
|
return;
|
495 |
|
|
}
|
496 |
|
|
|
497 |
|
|
tcx.basin.right_node = tcx.basin.bottom_node;
|
498 |
|
|
while (tcx.basin.right_node->next
|
499 |
|
|
&& tcx.basin.right_node->point->y < tcx.basin.right_node->next->point->y) {
|
500 |
|
|
tcx.basin.right_node = tcx.basin.right_node->next;
|
501 |
|
|
}
|
502 |
|
|
if (tcx.basin.right_node == tcx.basin.bottom_node) {
|
503 |
|
|
// No valid basins
|
504 |
|
|
return;
|
505 |
|
|
}
|
506 |
|
|
|
507 |
|
|
tcx.basin.width = tcx.basin.right_node->point->x - tcx.basin.left_node->point->x;
|
508 |
|
|
tcx.basin.left_highest = tcx.basin.left_node->point->y > tcx.basin.right_node->point->y;
|
509 |
|
|
|
510 |
|
|
FillBasinReq(tcx, tcx.basin.bottom_node);
|
511 |
|
|
}
|
512 |
|
|
|
513 |
|
|
void Sweep::FillBasinReq(SweepContext& tcx, Node* node)
|
514 |
|
|
{
|
515 |
|
|
// if shallow stop filling
|
516 |
|
|
if (IsShallow(tcx, *node)) {
|
517 |
|
|
return;
|
518 |
|
|
}
|
519 |
|
|
|
520 |
|
|
Fill(tcx, *node);
|
521 |
|
|
|
522 |
|
|
if (node->prev == tcx.basin.left_node && node->next == tcx.basin.right_node) {
|
523 |
|
|
return;
|
524 |
|
|
} else if (node->prev == tcx.basin.left_node) {
|
525 |
|
|
Orientation o = Orient2d(*node->point, *node->next->point, *node->next->next->point);
|
526 |
|
|
if (o == CW) {
|
527 |
|
|
return;
|
528 |
|
|
}
|
529 |
|
|
node = node->next;
|
530 |
|
|
} else if (node->next == tcx.basin.right_node) {
|
531 |
|
|
Orientation o = Orient2d(*node->point, *node->prev->point, *node->prev->prev->point);
|
532 |
|
|
if (o == CCW) {
|
533 |
|
|
return;
|
534 |
|
|
}
|
535 |
|
|
node = node->prev;
|
536 |
|
|
} else {
|
537 |
|
|
// Continue with the neighbor node with lowest Y value
|
538 |
|
|
if (node->prev->point->y < node->next->point->y) {
|
539 |
|
|
node = node->prev;
|
540 |
|
|
} else {
|
541 |
|
|
node = node->next;
|
542 |
|
|
}
|
543 |
|
|
}
|
544 |
|
|
|
545 |
|
|
FillBasinReq(tcx, node);
|
546 |
|
|
}
|
547 |
|
|
|
548 |
|
|
bool Sweep::IsShallow(SweepContext& tcx, Node& node)
|
549 |
|
|
{
|
550 |
|
|
double height;
|
551 |
|
|
|
552 |
|
|
if (tcx.basin.left_highest) {
|
553 |
|
|
height = tcx.basin.left_node->point->y - node.point->y;
|
554 |
|
|
} else {
|
555 |
|
|
height = tcx.basin.right_node->point->y - node.point->y;
|
556 |
|
|
}
|
557 |
|
|
|
558 |
|
|
// if shallow stop filling
|
559 |
|
|
if (tcx.basin.width > height) {
|
560 |
|
|
return true;
|
561 |
|
|
}
|
562 |
|
|
return false;
|
563 |
|
|
}
|
564 |
|
|
|
565 |
|
|
void Sweep::FillEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
|
566 |
|
|
{
|
567 |
|
|
if (tcx.edge_event.right) {
|
568 |
|
|
FillRightAboveEdgeEvent(tcx, edge, node);
|
569 |
|
|
} else {
|
570 |
|
|
FillLeftAboveEdgeEvent(tcx, edge, node);
|
571 |
|
|
}
|
572 |
|
|
}
|
573 |
|
|
|
574 |
|
|
void Sweep::FillRightAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
|
575 |
|
|
{
|
576 |
|
|
while (node->next->point->x < edge->p->x) {
|
577 |
|
|
// Check if next node is below the edge
|
578 |
|
|
if (Orient2d(*edge->q, *node->next->point, *edge->p) == CCW) {
|
579 |
|
|
FillRightBelowEdgeEvent(tcx, edge, *node);
|
580 |
|
|
} else {
|
581 |
|
|
node = node->next;
|
582 |
|
|
}
|
583 |
|
|
}
|
584 |
|
|
}
|
585 |
|
|
|
586 |
|
|
void Sweep::FillRightBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
587 |
|
|
{
|
588 |
|
|
if (node.point->x < edge->p->x) {
|
589 |
|
|
if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
|
590 |
|
|
// Concave
|
591 |
|
|
FillRightConcaveEdgeEvent(tcx, edge, node);
|
592 |
|
|
} else{
|
593 |
|
|
// Convex
|
594 |
|
|
FillRightConvexEdgeEvent(tcx, edge, node);
|
595 |
|
|
// Retry this one
|
596 |
|
|
FillRightBelowEdgeEvent(tcx, edge, node);
|
597 |
|
|
}
|
598 |
|
|
}
|
599 |
|
|
}
|
600 |
|
|
|
601 |
|
|
void Sweep::FillRightConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
602 |
|
|
{
|
603 |
|
|
Fill(tcx, *node.next);
|
604 |
|
|
if (node.next->point != edge->p) {
|
605 |
|
|
// Next above or below edge?
|
606 |
|
|
if (Orient2d(*edge->q, *node.next->point, *edge->p) == CCW) {
|
607 |
|
|
// Below
|
608 |
|
|
if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
|
609 |
|
|
// Next is concave
|
610 |
|
|
FillRightConcaveEdgeEvent(tcx, edge, node);
|
611 |
|
|
} else {
|
612 |
|
|
// Next is convex
|
613 |
|
|
}
|
614 |
|
|
}
|
615 |
|
|
}
|
616 |
|
|
|
617 |
|
|
}
|
618 |
|
|
|
619 |
|
|
void Sweep::FillRightConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
620 |
|
|
{
|
621 |
|
|
// Next concave or convex?
|
622 |
|
|
if (Orient2d(*node.next->point, *node.next->next->point, *node.next->next->next->point) == CCW) {
|
623 |
|
|
// Concave
|
624 |
|
|
FillRightConcaveEdgeEvent(tcx, edge, *node.next);
|
625 |
|
|
} else{
|
626 |
|
|
// Convex
|
627 |
|
|
// Next above or below edge?
|
628 |
|
|
if (Orient2d(*edge->q, *node.next->next->point, *edge->p) == CCW) {
|
629 |
|
|
// Below
|
630 |
|
|
FillRightConvexEdgeEvent(tcx, edge, *node.next);
|
631 |
|
|
} else{
|
632 |
|
|
// Above
|
633 |
|
|
}
|
634 |
|
|
}
|
635 |
|
|
}
|
636 |
|
|
|
637 |
|
|
void Sweep::FillLeftAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
|
638 |
|
|
{
|
639 |
|
|
while (node->prev->point->x > edge->p->x) {
|
640 |
|
|
// Check if next node is below the edge
|
641 |
|
|
if (Orient2d(*edge->q, *node->prev->point, *edge->p) == CW) {
|
642 |
|
|
FillLeftBelowEdgeEvent(tcx, edge, *node);
|
643 |
|
|
} else {
|
644 |
|
|
node = node->prev;
|
645 |
|
|
}
|
646 |
|
|
}
|
647 |
|
|
}
|
648 |
|
|
|
649 |
|
|
void Sweep::FillLeftBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
650 |
|
|
{
|
651 |
|
|
if (node.point->x > edge->p->x) {
|
652 |
|
|
if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) {
|
653 |
|
|
// Concave
|
654 |
|
|
FillLeftConcaveEdgeEvent(tcx, edge, node);
|
655 |
|
|
} else {
|
656 |
|
|
// Convex
|
657 |
|
|
FillLeftConvexEdgeEvent(tcx, edge, node);
|
658 |
|
|
// Retry this one
|
659 |
|
|
FillLeftBelowEdgeEvent(tcx, edge, node);
|
660 |
|
|
}
|
661 |
|
|
}
|
662 |
|
|
}
|
663 |
|
|
|
664 |
|
|
void Sweep::FillLeftConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
665 |
|
|
{
|
666 |
|
|
// Next concave or convex?
|
667 |
|
|
if (Orient2d(*node.prev->point, *node.prev->prev->point, *node.prev->prev->prev->point) == CW) {
|
668 |
|
|
// Concave
|
669 |
|
|
FillLeftConcaveEdgeEvent(tcx, edge, *node.prev);
|
670 |
|
|
} else{
|
671 |
|
|
// Convex
|
672 |
|
|
// Next above or below edge?
|
673 |
|
|
if (Orient2d(*edge->q, *node.prev->prev->point, *edge->p) == CW) {
|
674 |
|
|
// Below
|
675 |
|
|
FillLeftConvexEdgeEvent(tcx, edge, *node.prev);
|
676 |
|
|
} else{
|
677 |
|
|
// Above
|
678 |
|
|
}
|
679 |
|
|
}
|
680 |
|
|
}
|
681 |
|
|
|
682 |
|
|
void Sweep::FillLeftConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
|
683 |
|
|
{
|
684 |
|
|
Fill(tcx, *node.prev);
|
685 |
|
|
if (node.prev->point != edge->p) {
|
686 |
|
|
// Next above or below edge?
|
687 |
|
|
if (Orient2d(*edge->q, *node.prev->point, *edge->p) == CW) {
|
688 |
|
|
// Below
|
689 |
|
|
if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) {
|
690 |
|
|
// Next is concave
|
691 |
|
|
FillLeftConcaveEdgeEvent(tcx, edge, node);
|
692 |
|
|
} else{
|
693 |
|
|
// Next is convex
|
694 |
|
|
}
|
695 |
|
|
}
|
696 |
|
|
}
|
697 |
|
|
|
698 |
|
|
}
|
699 |
|
|
|
700 |
|
|
void Sweep::FlipEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p)
|
701 |
|
|
{
|
702 |
|
|
Triangle& ot = t->NeighborAcross(p);
|
703 |
|
|
Point& op = *ot.OppositePoint(*t, p);
|
704 |
|
|
|
705 |
|
|
if (&ot == NULL) {
|
706 |
|
|
// If we want to integrate the fillEdgeEvent do it here
|
707 |
|
|
// With current implementation we should never get here
|
708 |
|
|
//throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle");
|
709 |
|
|
assert(0);
|
710 |
|
|
}
|
711 |
|
|
|
712 |
|
|
if (InScanArea(p, *t->PointCCW(p), *t->PointCW(p), op)) {
|
713 |
|
|
// Lets rotate shared edge one vertex CW
|
714 |
|
|
RotateTrianglePair(*t, p, ot, op);
|
715 |
|
|
tcx.MapTriangleToNodes(*t);
|
716 |
|
|
tcx.MapTriangleToNodes(ot);
|
717 |
|
|
|
718 |
|
|
if (p == eq && op == ep) {
|
719 |
|
|
if (eq == *tcx.edge_event.constrained_edge->q && ep == *tcx.edge_event.constrained_edge->p) {
|
720 |
|
|
t->MarkConstrainedEdge(&ep, &eq);
|
721 |
|
|
ot.MarkConstrainedEdge(&ep, &eq);
|
722 |
|
|
Legalize(tcx, *t);
|
723 |
|
|
Legalize(tcx, ot);
|
724 |
|
|
} else {
|
725 |
|
|
// XXX: I think one of the triangles should be legalized here?
|
726 |
|
|
}
|
727 |
|
|
} else {
|
728 |
|
|
Orientation o = Orient2d(eq, op, ep);
|
729 |
|
|
t = &NextFlipTriangle(tcx, (int)o, *t, ot, p, op);
|
730 |
|
|
FlipEdgeEvent(tcx, ep, eq, t, p);
|
731 |
|
|
}
|
732 |
|
|
} else {
|
733 |
|
|
Point& newP = NextFlipPoint(ep, eq, ot, op);
|
734 |
|
|
FlipScanEdgeEvent(tcx, ep, eq, *t, ot, newP);
|
735 |
|
|
EdgeEvent(tcx, ep, eq, t, p);
|
736 |
|
|
}
|
737 |
|
|
}
|
738 |
|
|
|
739 |
|
|
Triangle& Sweep::NextFlipTriangle(SweepContext& tcx, int o, Triangle& t, Triangle& ot, Point& p, Point& op)
|
740 |
|
|
{
|
741 |
|
|
if (o == CCW) {
|
742 |
|
|
// ot is not crossing edge after flip
|
743 |
|
|
int edge_index = ot.EdgeIndex(&p, &op);
|
744 |
|
|
ot.delaunay_edge[edge_index] = true;
|
745 |
|
|
Legalize(tcx, ot);
|
746 |
|
|
ot.ClearDelunayEdges();
|
747 |
|
|
return t;
|
748 |
|
|
}
|
749 |
|
|
|
750 |
|
|
// t is not crossing edge after flip
|
751 |
|
|
int edge_index = t.EdgeIndex(&p, &op);
|
752 |
|
|
|
753 |
|
|
t.delaunay_edge[edge_index] = true;
|
754 |
|
|
Legalize(tcx, t);
|
755 |
|
|
t.ClearDelunayEdges();
|
756 |
|
|
return ot;
|
757 |
|
|
}
|
758 |
|
|
|
759 |
|
|
Point& Sweep::NextFlipPoint(Point& ep, Point& eq, Triangle& ot, Point& op)
|
760 |
|
|
{
|
761 |
|
|
Orientation o2d = Orient2d(eq, op, ep);
|
762 |
|
|
if (o2d == CW) {
|
763 |
|
|
// Right
|
764 |
|
|
return *ot.PointCCW(op);
|
765 |
|
|
} else if (o2d == CCW) {
|
766 |
|
|
// Left
|
767 |
|
|
return *ot.PointCW(op);
|
768 |
|
|
} else{
|
769 |
|
|
//throw new RuntimeException("[Unsupported] Opposing point on constrained edge");
|
770 |
|
|
assert(0);
|
771 |
|
|
}
|
772 |
|
|
}
|
773 |
|
|
|
774 |
|
|
void Sweep::FlipScanEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle,
|
775 |
|
|
Triangle& t, Point& p)
|
776 |
|
|
{
|
777 |
|
|
Triangle& ot = t.NeighborAcross(p);
|
778 |
|
|
Point& op = *ot.OppositePoint(t, p);
|
779 |
|
|
|
780 |
|
|
if (&t.NeighborAcross(p) == NULL) {
|
781 |
|
|
// If we want to integrate the fillEdgeEvent do it here
|
782 |
|
|
// With current implementation we should never get here
|
783 |
|
|
//throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle");
|
784 |
|
|
assert(0);
|
785 |
|
|
}
|
786 |
|
|
|
787 |
|
|
if (InScanArea(eq, *flip_triangle.PointCCW(eq), *flip_triangle.PointCW(eq), op)) {
|
788 |
|
|
// flip with new edge op->eq
|
789 |
|
|
FlipEdgeEvent(tcx, eq, op, &ot, op);
|
790 |
|
|
// TODO: Actually I just figured out that it should be possible to
|
791 |
|
|
// improve this by getting the next ot and op before the the above
|
792 |
|
|
// flip and continue the flipScanEdgeEvent here
|
793 |
|
|
// set new ot and op here and loop back to inScanArea test
|
794 |
|
|
// also need to set a new flip_triangle first
|
795 |
|
|
// Turns out at first glance that this is somewhat complicated
|
796 |
|
|
// so it will have to wait.
|
797 |
|
|
} else{
|
798 |
|
|
Point& newP = NextFlipPoint(ep, eq, ot, op);
|
799 |
|
|
FlipScanEdgeEvent(tcx, ep, eq, flip_triangle, ot, newP);
|
800 |
|
|
}
|
801 |
|
|
}
|
802 |
|
|
|
803 |
|
|
Sweep::~Sweep() {
|
804 |
|
|
|
805 |
|
|
// Clean up memory
|
806 |
|
|
for(unsigned int i = 0; i < nodes_.size(); i++) {
|
807 |
|
|
delete nodes_[i];
|
808 |
|
|
}
|
809 |
|
|
|
810 |
|
|
}
|
811 |
|
|
|
812 |
|
|
}
|
813 |
|
|
|