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