diff --git a/polygon/poly2tri/common/shapes.h b/polygon/poly2tri/common/shapes.h index c65f485eaf..e8dd811c23 100644 --- a/polygon/poly2tri/common/shapes.h +++ b/polygon/poly2tri/common/shapes.h @@ -56,7 +56,7 @@ struct Point std::vector edge_list; /// Construct using coordinates. - Point( double x, double y ) : x( x ), y( y ) {} + Point( double ax, double ay ) : x( ax ), y( ay ) {} /// Set this point to all zeros. void set_zero() diff --git a/polygon/poly2tri/sweep/advancing_front.cc b/polygon/poly2tri/sweep/advancing_front.cc index 019df4a6eb..b87dfc5451 100644 --- a/polygon/poly2tri/sweep/advancing_front.cc +++ b/polygon/poly2tri/sweep/advancing_front.cc @@ -31,79 +31,105 @@ #include "advancing_front.h" namespace p2t { - -AdvancingFront::AdvancingFront(Node& head, Node& tail) +AdvancingFront::AdvancingFront( Node& head, Node& tail ) { - head_ = &head; - tail_ = &tail; - search_node_ = &head; + head_ = &head; + tail_ = &tail; + search_node_ = &head; } -Node* AdvancingFront::LocateNode(const double& x) -{ - Node* node = search_node_; - if (x < node->value) { - while ((node = node->prev) != NULL) { - if (x >= node->value) { +Node* AdvancingFront::LocateNode( const double& x ) +{ + Node* node = search_node_; + + if( x < node->value ) + { + while( (node = node->prev) != NULL ) + { + if( x >= node->value ) + { + search_node_ = node; + return node; + } + } + } + else + { + while( (node = node->next) != NULL ) + { + if( x < node->value ) + { + search_node_ = node->prev; + return node->prev; + } + } + } + + return NULL; +} + + +Node* AdvancingFront::FindSearchNode( const double& x ) +{ + (void) x; // suppress compiler warnings "unused parameter 'x'" + // TODO: implement BST index + return search_node_; +} + + +Node* AdvancingFront::LocatePoint( const Point* point ) +{ + const double px = point->x; + Node* node = FindSearchNode( px ); + const double nx = node->point->x; + + if( px == nx ) + { + if( point != node->point ) + { + // We might have two nodes with same x value for a short time + if( point == node->prev->point ) + { + node = node->prev; + } + else if( point == node->next->point ) + { + node = node->next; + } + else + { + assert( 0 ); + } + } + } + else if( px < nx ) + { + while( (node = node->prev) != NULL ) + { + if( point == node->point ) + { + break; + } + } + } + else + { + while( (node = node->next) != NULL ) + { + if( point == node->point ) + break; + } + } + + if( node ) search_node_ = node; - return node; - } - } - } else { - while ((node = node->next) != NULL) { - if (x < node->value) { - search_node_ = node->prev; - return node->prev; - } - } - } - return NULL; + + return node; } -Node* AdvancingFront::FindSearchNode(const double& x) -{ - (void)x; // suppress compiler warnings "unused parameter 'x'" - // TODO: implement BST index - return search_node_; -} - -Node* AdvancingFront::LocatePoint(const Point* point) -{ - const double px = point->x; - Node* node = FindSearchNode(px); - const double nx = node->point->x; - - if (px == nx) { - if (point != node->point) { - // We might have two nodes with same x value for a short time - if (point == node->prev->point) { - node = node->prev; - } else if (point == node->next->point) { - node = node->next; - } else { - assert(0); - } - } - } else if (px < nx) { - while ((node = node->prev) != NULL) { - if (point == node->point) { - break; - } - } - } else { - while ((node = node->next) != NULL) { - if (point == node->point) - break; - } - } - if(node) search_node_ = node; - return node; -} AdvancingFront::~AdvancingFront() { } - } - diff --git a/polygon/poly2tri/sweep/advancing_front.h b/polygon/poly2tri/sweep/advancing_front.h index bab73d449c..179d22c82f 100644 --- a/polygon/poly2tri/sweep/advancing_front.h +++ b/polygon/poly2tri/sweep/advancing_front.h @@ -35,84 +35,91 @@ #include "../common/shapes.h" namespace p2t { - struct Node; // Advancing front node -struct Node { - Point* point; - Triangle* triangle; +struct Node +{ + Point* point; + Triangle* triangle; - Node* next; - Node* prev; + Node* next; + Node* prev; - double value; + double value; - Node(Point& p) : point(&p), triangle(NULL), next(NULL), prev(NULL), value(p.x) - { - } - - Node(Point& p, Triangle& t) : point(&p), triangle(&t), next(NULL), prev(NULL), value(p.x) - { - } + Node( Point& p ) : point( &p ), triangle( NULL ), next( NULL ), prev( NULL ), value( p.x ) + { + } + Node( Point& p, Triangle& t ) : point( &p ), triangle( &t ), next( NULL ), prev( NULL ), value( + p.x ) + { + } }; // Advancing front -class AdvancingFront { +class AdvancingFront +{ public: -AdvancingFront(Node& head, Node& tail); + AdvancingFront( Node& head, Node& tail ); // Destructor -~AdvancingFront(); + ~AdvancingFront(); -Node* head(); -void set_head(Node* node); -Node* tail(); -void set_tail(Node* node); -Node* search(); -void set_search(Node* node); + Node* head(); + void set_head( Node* node ); + Node* tail(); + void set_tail( Node* node ); + Node* search(); + void set_search( Node* node ); /// Locate insertion point along advancing front -Node* LocateNode(const double& x); + Node* LocateNode( const double& x ); -Node* LocatePoint(const Point* point); + Node* LocatePoint( const Point* point ); private: -Node* head_, *tail_, *search_node_; + Node* head_, * tail_, * search_node_; -Node* FindSearchNode(const double& x); + Node* FindSearchNode( const double& x ); }; inline Node* AdvancingFront::head() { - return head_; + return head_; } -inline void AdvancingFront::set_head(Node* node) + + +inline void AdvancingFront::set_head( Node* node ) { - head_ = node; + head_ = node; } + inline Node* AdvancingFront::tail() { - return tail_; + return tail_; } -inline void AdvancingFront::set_tail(Node* node) + + +inline void AdvancingFront::set_tail( Node* node ) { - tail_ = node; + tail_ = node; } + inline Node* AdvancingFront::search() { - return search_node_; + return search_node_; } -inline void AdvancingFront::set_search(Node* node) + +inline void AdvancingFront::set_search( Node* node ) { - search_node_ = node; + search_node_ = node; } - } #endif diff --git a/polygon/poly2tri/sweep/cdt.cc b/polygon/poly2tri/sweep/cdt.cc index d7838257c7..0cb910b6c2 100644 --- a/polygon/poly2tri/sweep/cdt.cc +++ b/polygon/poly2tri/sweep/cdt.cc @@ -1,4 +1,4 @@ -/* +/* * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors * http://code.google.com/p/poly2tri/ * @@ -31,42 +31,46 @@ #include "cdt.h" namespace p2t { - -CDT::CDT(std::vector polyline) +CDT::CDT( std::vector polyline ) { - sweep_context_ = new SweepContext(polyline); - sweep_ = new Sweep; + sweep_context_ = new SweepContext( polyline ); + sweep_ = new Sweep; } -void CDT::AddHole(std::vector polyline) + +void CDT::AddHole( std::vector polyline ) { - sweep_context_->AddHole(polyline); + sweep_context_->AddHole( polyline ); } -void CDT::AddPoint(Point* point) { - sweep_context_->AddPoint(point); + +void CDT::AddPoint( Point* point ) +{ + sweep_context_->AddPoint( point ); } + void CDT::Triangulate() { - sweep_->Triangulate(*sweep_context_); + sweep_->Triangulate( *sweep_context_ ); } + std::vector CDT::GetTriangles() { - return sweep_context_->GetTriangles(); + return sweep_context_->GetTriangles(); } + std::list CDT::GetMap() { - return sweep_context_->GetMap(); + return sweep_context_->GetMap(); } + CDT::~CDT() { - delete sweep_context_; - delete sweep_; + delete sweep_context_; + delete sweep_; } - } - diff --git a/polygon/poly2tri/sweep/sweep.cc b/polygon/poly2tri/sweep/sweep.cc index 75e7adf9ee..095437b7f6 100644 --- a/polygon/poly2tri/sweep/sweep.cc +++ b/polygon/poly2tri/sweep/sweep.cc @@ -35,783 +35,989 @@ #include "../common/utils.h" namespace p2t { - // Triangulate simple polygon with holes -void Sweep::Triangulate(SweepContext& tcx) +void Sweep::Triangulate( SweepContext& tcx ) { - tcx.InitTriangulation(); - tcx.CreateAdvancingFront(nodes_); - // Sweep points; build mesh - SweepPoints(tcx); - // Clean up - FinalizationPolygon(tcx); + tcx.InitTriangulation(); + tcx.CreateAdvancingFront( nodes_ ); + // Sweep points; build mesh + SweepPoints( tcx ); + // Clean up + FinalizationPolygon( tcx ); } -void Sweep::SweepPoints(SweepContext& tcx) + +void Sweep::SweepPoints( SweepContext& tcx ) { - for (int i = 1; i < tcx.point_count(); i++) { - Point& point = *tcx.GetPoint(i); - Node* node = &PointEvent(tcx, point); - for (unsigned int i = 0; i < point.edge_list.size(); i++) { - EdgeEvent(tcx, point.edge_list[i], node); + for( int jj = 1; jj < tcx.point_count(); jj++ ) + { + Point& point = *tcx.GetPoint( jj ); + Node* node = &PointEvent( tcx, point ); + + for( unsigned int i = 0; i < point.edge_list.size(); i++ ) + { + EdgeEvent( tcx, point.edge_list[i], node ); + } } - } } -void Sweep::FinalizationPolygon(SweepContext& tcx) + +void Sweep::FinalizationPolygon( SweepContext& tcx ) { - // Get an Internal triangle to start with - Triangle* t = tcx.front()->head()->next->triangle; - Point* p = tcx.front()->head()->next->point; - while (!t->GetConstrainedEdgeCW(*p)) { - t = t->NeighborCCW(*p); - } + // Get an Internal triangle to start with + Triangle* t = tcx.front()->head()->next->triangle; + Point* p = tcx.front()->head()->next->point; - // Collect interior triangles constrained by edges - tcx.MeshClean(*t); -} - -Node& Sweep::PointEvent(SweepContext& tcx, Point& point) -{ - Node& node = tcx.LocateNode(point); - Node& new_node = NewFrontTriangle(tcx, point, node); - - // Only need to check +epsilon since point never have smaller - // x value than node due to how we fetch nodes from the front - if (point.x <= node.point->x + EPSILON) { - Fill(tcx, node); - } - - //tcx.AddNode(new_node); - - FillAdvancingFront(tcx, new_node); - return new_node; -} - -void Sweep::EdgeEvent(SweepContext& tcx, Edge* edge, Node* node) -{ - tcx.edge_event.constrained_edge = edge; - tcx.edge_event.right = (edge->p->x > edge->q->x); - - if (IsEdgeSideOfTriangle(*node->triangle, *edge->p, *edge->q)) { - return; - } - - // For now we will do all needed filling - // TODO: integrate with flip process might give some better performance - // but for now this avoid the issue with cases that needs both flips and fills - FillEdgeEvent(tcx, edge, node); - EdgeEvent(tcx, *edge->p, *edge->q, node->triangle, *edge->q); -} - -void Sweep::EdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point) -{ - if (IsEdgeSideOfTriangle(*triangle, ep, eq)) { - return; - } - - Point* p1 = triangle->PointCCW(point); - Orientation o1 = Orient2d(eq, *p1, ep); - if (o1 == COLLINEAR) { - if( triangle->Contains(&eq, p1)) { - triangle->MarkConstrainedEdge(&eq, p1 ); - // We are modifying the constraint maybe it would be better to - // not change the given constraint and just keep a variable for the new constraint - tcx.edge_event.constrained_edge->q = p1; - triangle = &triangle->NeighborAcross(point); - EdgeEvent( tcx, ep, *p1, triangle, *p1 ); - } else { - std::runtime_error("EdgeEvent - collinear points not supported"); - assert(0); + while( !t->GetConstrainedEdgeCW( *p ) ) + { + t = t->NeighborCCW( *p ); } - return; - } - Point* p2 = triangle->PointCW(point); - Orientation o2 = Orient2d(eq, *p2, ep); - if (o2 == COLLINEAR) { - if( triangle->Contains(&eq, p2)) { - triangle->MarkConstrainedEdge(&eq, p2 ); - // We are modifying the constraint maybe it would be better to - // not change the given constraint and just keep a variable for the new constraint - tcx.edge_event.constrained_edge->q = p2; - triangle = &triangle->NeighborAcross(point); - EdgeEvent( tcx, ep, *p2, triangle, *p2 ); - } else { - std::runtime_error("EdgeEvent - collinear points not supported"); - assert(0); + // Collect interior triangles constrained by edges + tcx.MeshClean( *t ); +} + + +Node& Sweep::PointEvent( SweepContext& tcx, Point& point ) +{ + Node& node = tcx.LocateNode( point ); + Node& new_node = NewFrontTriangle( tcx, point, node ); + + // Only need to check +epsilon since point never have smaller + // x value than node due to how we fetch nodes from the front + if( point.x <= node.point->x + EPSILON ) + { + Fill( tcx, node ); } - return; - } - if (o1 == o2) { - // Need to decide if we are rotating CW or CCW to get to a triangle - // that will cross edge - if (o1 == CW) { - triangle = triangle->NeighborCCW(point); - } else{ - triangle = triangle->NeighborCW(point); + // tcx.AddNode(new_node); + + FillAdvancingFront( tcx, new_node ); + return new_node; +} + + +void Sweep::EdgeEvent( SweepContext& tcx, Edge* edge, Node* node ) +{ + tcx.edge_event.constrained_edge = edge; + tcx.edge_event.right = (edge->p->x > edge->q->x); + + if( IsEdgeSideOfTriangle( *node->triangle, *edge->p, *edge->q ) ) + { + return; } - EdgeEvent(tcx, ep, eq, triangle, point); - } else { - // This triangle crosses constraint so lets flippin start! - FlipEdgeEvent(tcx, ep, eq, triangle, point); - } + + // For now we will do all needed filling + // TODO: integrate with flip process might give some better performance + // but for now this avoid the issue with cases that needs both flips and fills + FillEdgeEvent( tcx, edge, node ); + EdgeEvent( tcx, *edge->p, *edge->q, node->triangle, *edge->q ); } -bool Sweep::IsEdgeSideOfTriangle(Triangle& triangle, Point& ep, Point& eq) -{ - int index = triangle.EdgeIndex(&ep, &eq); - if (index != -1) { - triangle.MarkConstrainedEdge(index); - Triangle* t = triangle.GetNeighbor(index); - if (t) { - t->MarkConstrainedEdge(&ep, &eq); +void Sweep::EdgeEvent( SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point ) +{ + if( IsEdgeSideOfTriangle( *triangle, ep, eq ) ) + { + return; } - return true; - } - return false; -} -Node& Sweep::NewFrontTriangle(SweepContext& tcx, Point& point, Node& node) -{ - Triangle* triangle = new Triangle(point, *node.point, *node.next->point); + Point* p1 = triangle->PointCCW( point ); + Orientation o1 = Orient2d( eq, *p1, ep ); - triangle->MarkNeighbor(*node.triangle); - tcx.AddToMap(triangle); - - Node* new_node = new Node(point); - nodes_.push_back(new_node); - - new_node->next = node.next; - new_node->prev = &node; - node.next->prev = new_node; - node.next = new_node; - - if (!Legalize(tcx, *triangle)) { - tcx.MapTriangleToNodes(*triangle); - } - - return *new_node; -} - -void Sweep::Fill(SweepContext& tcx, Node& node) -{ - Triangle* triangle = new Triangle(*node.prev->point, *node.point, *node.next->point); - - // TODO: should copy the constrained_edge value from neighbor triangles - // for now constrained_edge values are copied during the legalize - triangle->MarkNeighbor(*node.prev->triangle); - triangle->MarkNeighbor(*node.triangle); - - tcx.AddToMap(triangle); - - // Update the advancing front - node.prev->next = node.next; - node.next->prev = node.prev; - - // If it was legalized the triangle has already been mapped - if (!Legalize(tcx, *triangle)) { - tcx.MapTriangleToNodes(*triangle); - } - -} - -void Sweep::FillAdvancingFront(SweepContext& tcx, Node& n) -{ - - // Fill right holes - Node* node = n.next; - - while (node->next) { - // if HoleAngle exceeds 90 degrees then break. - if (LargeHole_DontFill(node)) break; - Fill(tcx, *node); - node = node->next; - } - - // Fill left holes - node = n.prev; - - while (node->prev) { - // if HoleAngle exceeds 90 degrees then break. - if (LargeHole_DontFill(node)) break; - Fill(tcx, *node); - node = node->prev; - } - - // Fill right basins - if (n.next && n.next->next) { - double angle = BasinAngle(n); - if (angle < PI_3div4) { - FillBasin(tcx, n); - } - } -} - -// True if HoleAngle exceeds 90 degrees. -bool Sweep::LargeHole_DontFill(Node* node) { - - Node* nextNode = node->next; - Node* prevNode = node->prev; - if (!AngleExceeds90Degrees(node->point, nextNode->point, prevNode->point)) - return false; - - // Check additional points on front. - Node* next2Node = nextNode->next; - // "..Plus.." because only want angles on same side as point being added. - if ((next2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, next2Node->point, prevNode->point)) - return false; - - Node* prev2Node = prevNode->prev; - // "..Plus.." because only want angles on same side as point being added. - if ((prev2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, nextNode->point, prev2Node->point)) - return false; - - return true; -} - -bool Sweep::AngleExceeds90Degrees(Point* origin, Point* pa, Point* pb) { - double angle = Angle(*origin, *pa, *pb); - bool exceeds90Degrees = ((angle > PI_div2) || (angle < -PI_div2)); - return exceeds90Degrees; -} - -bool Sweep::AngleExceedsPlus90DegreesOrIsNegative(Point* origin, Point* pa, Point* pb) { - double angle = Angle(*origin, *pa, *pb); - bool exceedsPlus90DegreesOrIsNegative = (angle > PI_div2) || (angle < 0); - return exceedsPlus90DegreesOrIsNegative; -} - -double Sweep::Angle(Point& origin, Point& pa, Point& pb) { - /* Complex plane - * ab = cosA +i*sinA - * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) - * atan2(y,x) computes the principal value of the argument function - * applied to the complex number x+iy - * Where x = ax*bx + ay*by - * y = ax*by - ay*bx - */ - double px = origin.x; - double py = origin.y; - double ax = pa.x- px; - double ay = pa.y - py; - double bx = pb.x - px; - double by = pb.y - py; - double x = ax * by - ay * bx; - double y = ax * bx + ay * by; - double angle = atan2(x, y); - return angle; -} - -double Sweep::BasinAngle(Node& node) -{ - double ax = node.point->x - node.next->next->point->x; - double ay = node.point->y - node.next->next->point->y; - return atan2(ay, ax); -} - -double Sweep::HoleAngle(Node& node) -{ - /* Complex plane - * ab = cosA +i*sinA - * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) - * atan2(y,x) computes the principal value of the argument function - * applied to the complex number x+iy - * Where x = ax*bx + ay*by - * y = ax*by - ay*bx - */ - double ax = node.next->point->x - node.point->x; - double ay = node.next->point->y - node.point->y; - double bx = node.prev->point->x - node.point->x; - double by = node.prev->point->y - node.point->y; - return atan2(ax * by - ay * bx, ax * bx + ay * by); -} - -bool Sweep::Legalize(SweepContext& tcx, Triangle& t) -{ - // To legalize a triangle we start by finding if any of the three edges - // violate the Delaunay condition - for (int i = 0; i < 3; i++) { - if (t.delaunay_edge[i]) - continue; - - Triangle* ot = t.GetNeighbor(i); - - if (ot) { - Point* p = t.GetPoint(i); - Point* op = ot->OppositePoint(t, *p); - int oi = ot->Index(op); - - // If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization) - // then we should not try to legalize - if (ot->constrained_edge[oi] || ot->delaunay_edge[oi]) { - t.constrained_edge[i] = ot->constrained_edge[oi]; - continue; - } - - bool inside = Incircle(*p, *t.PointCCW(*p), *t.PointCW(*p), *op); - - if (inside) { - // Lets mark this shared edge as Delaunay - t.delaunay_edge[i] = true; - ot->delaunay_edge[oi] = true; - - // Lets rotate shared edge one vertex CW to legalize it - RotateTrianglePair(t, *p, *ot, *op); - - // We now got one valid Delaunay Edge shared by two triangles - // This gives us 4 new edges to check for Delaunay - - // Make sure that triangle to node mapping is done only one time for a specific triangle - bool not_legalized = !Legalize(tcx, t); - if (not_legalized) { - tcx.MapTriangleToNodes(t); + if( o1 == COLLINEAR ) + { + if( triangle->Contains( &eq, p1 ) ) + { + triangle->MarkConstrainedEdge( &eq, p1 ); + // We are modifying the constraint maybe it would be better to + // not change the given constraint and just keep a variable for the new constraint + tcx.edge_event.constrained_edge->q = p1; + triangle = &triangle->NeighborAcross( point ); + EdgeEvent( tcx, ep, *p1, triangle, *p1 ); + } + else + { + std::runtime_error( "EdgeEvent - collinear points not supported" ); + assert( 0 ); } - not_legalized = !Legalize(tcx, *ot); - if (not_legalized) - tcx.MapTriangleToNodes(*ot); + return; + } - // Reset the Delaunay edges, since they only are valid Delaunay edges - // until we add a new triangle or point. - // XXX: need to think about this. Can these edges be tried after we - // return to previous recursive level? - t.delaunay_edge[i] = false; - ot->delaunay_edge[oi] = false; + Point* p2 = triangle->PointCW( point ); + Orientation o2 = Orient2d( eq, *p2, ep ); + + if( o2 == COLLINEAR ) + { + if( triangle->Contains( &eq, p2 ) ) + { + triangle->MarkConstrainedEdge( &eq, p2 ); + // We are modifying the constraint maybe it would be better to + // not change the given constraint and just keep a variable for the new constraint + tcx.edge_event.constrained_edge->q = p2; + triangle = &triangle->NeighborAcross( point ); + EdgeEvent( tcx, ep, *p2, triangle, *p2 ); + } + else + { + std::runtime_error( "EdgeEvent - collinear points not supported" ); + assert( 0 ); + } + + return; + } + + if( o1 == o2 ) + { + // Need to decide if we are rotating CW or CCW to get to a triangle + // that will cross edge + if( o1 == CW ) + { + triangle = triangle->NeighborCCW( point ); + } + else + { + triangle = triangle->NeighborCW( point ); + } + + EdgeEvent( tcx, ep, eq, triangle, point ); + } + else + { + // This triangle crosses constraint so lets flippin start! + FlipEdgeEvent( tcx, ep, eq, triangle, point ); + } +} + + +bool Sweep::IsEdgeSideOfTriangle( Triangle& triangle, Point& ep, Point& eq ) +{ + int index = triangle.EdgeIndex( &ep, &eq ); + + if( index != -1 ) + { + triangle.MarkConstrainedEdge( index ); + Triangle* t = triangle.GetNeighbor( index ); + + if( t ) + { + t->MarkConstrainedEdge( &ep, &eq ); + } - // If triangle have been legalized no need to check the other edges since - // the recursive legalization will handles those so we can end here. return true; - } } - } - return false; -} -bool Sweep::Incircle(Point& pa, Point& pb, Point& pc, Point& pd) -{ - double adx = pa.x - pd.x; - double ady = pa.y - pd.y; - double bdx = pb.x - pd.x; - double bdy = pb.y - pd.y; - - double adxbdy = adx * bdy; - double bdxady = bdx * ady; - double oabd = adxbdy - bdxady; - - if (oabd <= 0) return false; - - double cdx = pc.x - pd.x; - double cdy = pc.y - pd.y; - - double cdxady = cdx * ady; - double adxcdy = adx * cdy; - double ocad = cdxady - adxcdy; - - if (ocad <= 0) - return false; - - double bdxcdy = bdx * cdy; - double cdxbdy = cdx * bdy; - - double alift = adx * adx + ady * ady; - double blift = bdx * bdx + bdy * bdy; - double clift = cdx * cdx + cdy * cdy; - - double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd; - - return det > 0; } -void Sweep::RotateTrianglePair(Triangle& t, Point& p, Triangle& ot, Point& op) + +Node& Sweep::NewFrontTriangle( SweepContext& tcx, Point& point, Node& node ) { - Triangle* n1, *n2, *n3, *n4; - n1 = t.NeighborCCW(p); - n2 = t.NeighborCW(p); - n3 = ot.NeighborCCW(op); - n4 = ot.NeighborCW(op); + Triangle* triangle = new Triangle( point, *node.point, *node.next->point ); - bool ce1, ce2, ce3, ce4; - ce1 = t.GetConstrainedEdgeCCW(p); - ce2 = t.GetConstrainedEdgeCW(p); - ce3 = ot.GetConstrainedEdgeCCW(op); - ce4 = ot.GetConstrainedEdgeCW(op); + triangle->MarkNeighbor( *node.triangle ); + tcx.AddToMap( triangle ); - bool de1, de2, de3, de4; - de1 = t.GetDelunayEdgeCCW(p); - de2 = t.GetDelunayEdgeCW(p); - de3 = ot.GetDelunayEdgeCCW(op); - de4 = ot.GetDelunayEdgeCW(op); + Node* new_node = new Node( point ); + nodes_.push_back( new_node ); - t.Legalize(p, op); - ot.Legalize(op, p); + new_node->next = node.next; + new_node->prev = &node; + node.next->prev = new_node; + node.next = new_node; - // Remap delaunay_edge - ot.SetDelunayEdgeCCW(p, de1); - t.SetDelunayEdgeCW(p, de2); - t.SetDelunayEdgeCCW(op, de3); - ot.SetDelunayEdgeCW(op, de4); - - // Remap constrained_edge - ot.SetConstrainedEdgeCCW(p, ce1); - t.SetConstrainedEdgeCW(p, ce2); - t.SetConstrainedEdgeCCW(op, ce3); - ot.SetConstrainedEdgeCW(op, ce4); - - // Remap neighbors - // XXX: might optimize the markNeighbor by keeping track of - // what side should be assigned to what neighbor after the - // rotation. Now mark neighbor does lots of testing to find - // the right side. - t.ClearNeighbors(); - ot.ClearNeighbors(); - if (n1) ot.MarkNeighbor(*n1); - if (n2) t.MarkNeighbor(*n2); - if (n3) t.MarkNeighbor(*n3); - if (n4) ot.MarkNeighbor(*n4); - t.MarkNeighbor(ot); -} - -void Sweep::FillBasin(SweepContext& tcx, Node& node) -{ - if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { - tcx.basin.left_node = node.next->next; - } else { - tcx.basin.left_node = node.next; - } - - // Find the bottom and right node - tcx.basin.bottom_node = tcx.basin.left_node; - while (tcx.basin.bottom_node->next - && tcx.basin.bottom_node->point->y >= tcx.basin.bottom_node->next->point->y) { - tcx.basin.bottom_node = tcx.basin.bottom_node->next; - } - if (tcx.basin.bottom_node == tcx.basin.left_node) { - // No valid basin - return; - } - - tcx.basin.right_node = tcx.basin.bottom_node; - while (tcx.basin.right_node->next - && tcx.basin.right_node->point->y < tcx.basin.right_node->next->point->y) { - tcx.basin.right_node = tcx.basin.right_node->next; - } - if (tcx.basin.right_node == tcx.basin.bottom_node) { - // No valid basins - return; - } - - tcx.basin.width = tcx.basin.right_node->point->x - tcx.basin.left_node->point->x; - tcx.basin.left_highest = tcx.basin.left_node->point->y > tcx.basin.right_node->point->y; - - FillBasinReq(tcx, tcx.basin.bottom_node); -} - -void Sweep::FillBasinReq(SweepContext& tcx, Node* node) -{ - // if shallow stop filling - if (IsShallow(tcx, *node)) { - return; - } - - Fill(tcx, *node); - - if (node->prev == tcx.basin.left_node && node->next == tcx.basin.right_node) { - return; - } else if (node->prev == tcx.basin.left_node) { - Orientation o = Orient2d(*node->point, *node->next->point, *node->next->next->point); - if (o == CW) { - return; + if( !Legalize( tcx, *triangle ) ) + { + tcx.MapTriangleToNodes( *triangle ); } - node = node->next; - } else if (node->next == tcx.basin.right_node) { - Orientation o = Orient2d(*node->point, *node->prev->point, *node->prev->prev->point); - if (o == CCW) { - return; - } - node = node->prev; - } else { - // Continue with the neighbor node with lowest Y value - if (node->prev->point->y < node->next->point->y) { - node = node->prev; - } else { - node = node->next; - } - } - FillBasinReq(tcx, node); + return *new_node; } -bool Sweep::IsShallow(SweepContext& tcx, Node& node) + +void Sweep::Fill( SweepContext& tcx, Node& node ) { - double height; + Triangle* triangle = new Triangle( *node.prev->point, *node.point, *node.next->point ); - if (tcx.basin.left_highest) { - height = tcx.basin.left_node->point->y - node.point->y; - } else { - height = tcx.basin.right_node->point->y - node.point->y; - } + // TODO: should copy the constrained_edge value from neighbor triangles + // for now constrained_edge values are copied during the legalize + triangle->MarkNeighbor( *node.prev->triangle ); + triangle->MarkNeighbor( *node.triangle ); + + tcx.AddToMap( triangle ); + + // Update the advancing front + node.prev->next = node.next; + node.next->prev = node.prev; + + // If it was legalized the triangle has already been mapped + if( !Legalize( tcx, *triangle ) ) + { + tcx.MapTriangleToNodes( *triangle ); + } +} + + +void Sweep::FillAdvancingFront( SweepContext& tcx, Node& n ) +{ + // Fill right holes + Node* node = n.next; + + while( node->next ) + { + // if HoleAngle exceeds 90 degrees then break. + if( LargeHole_DontFill( node ) ) + break; + + Fill( tcx, *node ); + node = node->next; + } + + // Fill left holes + node = n.prev; + + while( node->prev ) + { + // if HoleAngle exceeds 90 degrees then break. + if( LargeHole_DontFill( node ) ) + break; + + Fill( tcx, *node ); + node = node->prev; + } + + // Fill right basins + if( n.next && n.next->next ) + { + double angle = BasinAngle( n ); + + if( angle < PI_3div4 ) + { + FillBasin( tcx, n ); + } + } +} + + +// True if HoleAngle exceeds 90 degrees. +bool Sweep::LargeHole_DontFill( Node* node ) +{ + Node* nextNode = node->next; + Node* prevNode = node->prev; + + if( !AngleExceeds90Degrees( node->point, nextNode->point, prevNode->point ) ) + return false; + + // Check additional points on front. + Node* next2Node = nextNode->next; + + // "..Plus.." because only want angles on same side as point being added. + if( (next2Node != NULL) + && !AngleExceedsPlus90DegreesOrIsNegative( node->point, next2Node->point, + prevNode->point ) ) + return false; + + Node* prev2Node = prevNode->prev; + + // "..Plus.." because only want angles on same side as point being added. + if( (prev2Node != NULL) + && !AngleExceedsPlus90DegreesOrIsNegative( node->point, nextNode->point, + prev2Node->point ) ) + return false; - // if shallow stop filling - if (tcx.basin.width > height) { return true; - } - return false; } -void Sweep::FillEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) + +bool Sweep::AngleExceeds90Degrees( Point* origin, Point* pa, Point* pb ) { - if (tcx.edge_event.right) { - FillRightAboveEdgeEvent(tcx, edge, node); - } else { - FillLeftAboveEdgeEvent(tcx, edge, node); - } + double angle = Angle( *origin, *pa, *pb ); + bool exceeds90Degrees = ( (angle > PI_div2) || (angle < -PI_div2) ); + + return exceeds90Degrees; } -void Sweep::FillRightAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) + +bool Sweep::AngleExceedsPlus90DegreesOrIsNegative( Point* origin, Point* pa, Point* pb ) { - while (node->next->point->x < edge->p->x) { - // Check if next node is below the edge - if (Orient2d(*edge->q, *node->next->point, *edge->p) == CCW) { - FillRightBelowEdgeEvent(tcx, edge, *node); - } else { - node = node->next; + double angle = Angle( *origin, *pa, *pb ); + bool exceedsPlus90DegreesOrIsNegative = (angle > PI_div2) || (angle < 0); + + return exceedsPlus90DegreesOrIsNegative; +} + + +double Sweep::Angle( Point& origin, Point& pa, Point& pb ) +{ + /* Complex plane + * ab = cosA +i*sinA + * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) + * atan2(y,x) computes the principal value of the argument function + * applied to the complex number x+iy + * Where x = ax*bx + ay*by + * y = ax*by - ay*bx + */ + double px = origin.x; + double py = origin.y; + double ax = pa.x - px; + double ay = pa.y - py; + double bx = pb.x - px; + double by = pb.y - py; + double x = ax * by - ay * bx; + double y = ax * bx + ay * by; + double angle = atan2( x, y ); + + return angle; +} + + +double Sweep::BasinAngle( Node& node ) +{ + double ax = node.point->x - node.next->next->point->x; + double ay = node.point->y - node.next->next->point->y; + + return atan2( ay, ax ); +} + + +double Sweep::HoleAngle( Node& node ) +{ + /* Complex plane + * ab = cosA +i*sinA + * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) + * atan2(y,x) computes the principal value of the argument function + * applied to the complex number x+iy + * Where x = ax*bx + ay*by + * y = ax*by - ay*bx + */ + double ax = node.next->point->x - node.point->x; + double ay = node.next->point->y - node.point->y; + double bx = node.prev->point->x - node.point->x; + double by = node.prev->point->y - node.point->y; + + return atan2( ax * by - ay * bx, ax * bx + ay * by ); +} + + +bool Sweep::Legalize( SweepContext& tcx, Triangle& t ) +{ + // To legalize a triangle we start by finding if any of the three edges + // violate the Delaunay condition + for( int i = 0; i < 3; i++ ) + { + if( t.delaunay_edge[i] ) + continue; + + Triangle* ot = t.GetNeighbor( i ); + + if( ot ) + { + Point* p = t.GetPoint( i ); + Point* op = ot->OppositePoint( t, *p ); + int oi = ot->Index( op ); + + // If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization) + // then we should not try to legalize + if( ot->constrained_edge[oi] || ot->delaunay_edge[oi] ) + { + t.constrained_edge[i] = ot->constrained_edge[oi]; + continue; + } + + bool inside = Incircle( *p, *t.PointCCW( *p ), *t.PointCW( *p ), *op ); + + if( inside ) + { + // Lets mark this shared edge as Delaunay + t.delaunay_edge[i] = true; + ot->delaunay_edge[oi] = true; + + // Lets rotate shared edge one vertex CW to legalize it + RotateTrianglePair( t, *p, *ot, *op ); + + // We now got one valid Delaunay Edge shared by two triangles + // This gives us 4 new edges to check for Delaunay + + // Make sure that triangle to node mapping is done only one time for a specific triangle + bool not_legalized = !Legalize( tcx, t ); + + if( not_legalized ) + { + tcx.MapTriangleToNodes( t ); + } + + not_legalized = !Legalize( tcx, *ot ); + + if( not_legalized ) + tcx.MapTriangleToNodes( *ot ); + + // Reset the Delaunay edges, since they only are valid Delaunay edges + // until we add a new triangle or point. + // XXX: need to think about this. Can these edges be tried after we + // return to previous recursive level? + t.delaunay_edge[i] = false; + ot->delaunay_edge[oi] = false; + + // If triangle have been legalized no need to check the other edges since + // the recursive legalization will handles those so we can end here. + return true; + } + } } - } + + return false; } -void Sweep::FillRightBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) + +bool Sweep::Incircle( Point& pa, Point& pb, Point& pc, Point& pd ) { - if (node.point->x < edge->p->x) { - if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { - // Concave - FillRightConcaveEdgeEvent(tcx, edge, node); - } else{ - // Convex - FillRightConvexEdgeEvent(tcx, edge, node); - // Retry this one - FillRightBelowEdgeEvent(tcx, edge, node); + double adx = pa.x - pd.x; + double ady = pa.y - pd.y; + double bdx = pb.x - pd.x; + double bdy = pb.y - pd.y; + + double adxbdy = adx * bdy; + double bdxady = bdx * ady; + double oabd = adxbdy - bdxady; + + if( oabd <= 0 ) + return false; + + double cdx = pc.x - pd.x; + double cdy = pc.y - pd.y; + + double cdxady = cdx * ady; + double adxcdy = adx * cdy; + double ocad = cdxady - adxcdy; + + if( ocad <= 0 ) + return false; + + double bdxcdy = bdx * cdy; + double cdxbdy = cdx * bdy; + + double alift = adx * adx + ady * ady; + double blift = bdx * bdx + bdy * bdy; + double clift = cdx * cdx + cdy * cdy; + + double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd; + + return det > 0; +} + + +void Sweep::RotateTrianglePair( Triangle& t, Point& p, Triangle& ot, Point& op ) +{ + Triangle* n1, * n2, * n3, * n4; + + n1 = t.NeighborCCW( p ); + n2 = t.NeighborCW( p ); + n3 = ot.NeighborCCW( op ); + n4 = ot.NeighborCW( op ); + + bool ce1, ce2, ce3, ce4; + ce1 = t.GetConstrainedEdgeCCW( p ); + ce2 = t.GetConstrainedEdgeCW( p ); + ce3 = ot.GetConstrainedEdgeCCW( op ); + ce4 = ot.GetConstrainedEdgeCW( op ); + + bool de1, de2, de3, de4; + de1 = t.GetDelunayEdgeCCW( p ); + de2 = t.GetDelunayEdgeCW( p ); + de3 = ot.GetDelunayEdgeCCW( op ); + de4 = ot.GetDelunayEdgeCW( op ); + + t.Legalize( p, op ); + ot.Legalize( op, p ); + + // Remap delaunay_edge + ot.SetDelunayEdgeCCW( p, de1 ); + t.SetDelunayEdgeCW( p, de2 ); + t.SetDelunayEdgeCCW( op, de3 ); + ot.SetDelunayEdgeCW( op, de4 ); + + // Remap constrained_edge + ot.SetConstrainedEdgeCCW( p, ce1 ); + t.SetConstrainedEdgeCW( p, ce2 ); + t.SetConstrainedEdgeCCW( op, ce3 ); + ot.SetConstrainedEdgeCW( op, ce4 ); + + // Remap neighbors + // XXX: might optimize the markNeighbor by keeping track of + // what side should be assigned to what neighbor after the + // rotation. Now mark neighbor does lots of testing to find + // the right side. + t.ClearNeighbors(); + ot.ClearNeighbors(); + + if( n1 ) + ot.MarkNeighbor( *n1 ); + + if( n2 ) + t.MarkNeighbor( *n2 ); + + if( n3 ) + t.MarkNeighbor( *n3 ); + + if( n4 ) + ot.MarkNeighbor( *n4 ); + + t.MarkNeighbor( ot ); +} + + +void Sweep::FillBasin( SweepContext& tcx, Node& node ) +{ + if( Orient2d( *node.point, *node.next->point, *node.next->next->point ) == CCW ) + { + tcx.basin.left_node = node.next->next; } - } -} - -void Sweep::FillRightConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) -{ - Fill(tcx, *node.next); - if (node.next->point != edge->p) { - // Next above or below edge? - if (Orient2d(*edge->q, *node.next->point, *edge->p) == CCW) { - // Below - if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { - // Next is concave - FillRightConcaveEdgeEvent(tcx, edge, node); - } else { - // Next is convex - } + else + { + tcx.basin.left_node = node.next; } - } -} + // Find the bottom and right node + tcx.basin.bottom_node = tcx.basin.left_node; -void Sweep::FillRightConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) -{ - // Next concave or convex? - if (Orient2d(*node.next->point, *node.next->next->point, *node.next->next->next->point) == CCW) { - // Concave - FillRightConcaveEdgeEvent(tcx, edge, *node.next); - } else{ - // Convex - // Next above or below edge? - if (Orient2d(*edge->q, *node.next->next->point, *edge->p) == CCW) { - // Below - FillRightConvexEdgeEvent(tcx, edge, *node.next); - } else{ - // Above + while( tcx.basin.bottom_node->next + && tcx.basin.bottom_node->point->y >= tcx.basin.bottom_node->next->point->y ) + { + tcx.basin.bottom_node = tcx.basin.bottom_node->next; } - } -} -void Sweep::FillLeftAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) -{ - while (node->prev->point->x > edge->p->x) { - // Check if next node is below the edge - if (Orient2d(*edge->q, *node->prev->point, *edge->p) == CW) { - FillLeftBelowEdgeEvent(tcx, edge, *node); - } else { - node = node->prev; + if( tcx.basin.bottom_node == tcx.basin.left_node ) + { + // No valid basin + return; } - } -} -void Sweep::FillLeftBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) -{ - if (node.point->x > edge->p->x) { - if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) { - // Concave - FillLeftConcaveEdgeEvent(tcx, edge, node); - } else { - // Convex - FillLeftConvexEdgeEvent(tcx, edge, node); - // Retry this one - FillLeftBelowEdgeEvent(tcx, edge, node); + tcx.basin.right_node = tcx.basin.bottom_node; + + while( tcx.basin.right_node->next + && tcx.basin.right_node->point->y < tcx.basin.right_node->next->point->y ) + { + tcx.basin.right_node = tcx.basin.right_node->next; } - } -} -void Sweep::FillLeftConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) -{ - // Next concave or convex? - if (Orient2d(*node.prev->point, *node.prev->prev->point, *node.prev->prev->prev->point) == CW) { - // Concave - FillLeftConcaveEdgeEvent(tcx, edge, *node.prev); - } else{ - // Convex - // Next above or below edge? - if (Orient2d(*edge->q, *node.prev->prev->point, *edge->p) == CW) { - // Below - FillLeftConvexEdgeEvent(tcx, edge, *node.prev); - } else{ - // Above + if( tcx.basin.right_node == tcx.basin.bottom_node ) + { + // No valid basins + return; } - } + + tcx.basin.width = tcx.basin.right_node->point->x - tcx.basin.left_node->point->x; + tcx.basin.left_highest = tcx.basin.left_node->point->y > tcx.basin.right_node->point->y; + + FillBasinReq( tcx, tcx.basin.bottom_node ); } -void Sweep::FillLeftConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) + +void Sweep::FillBasinReq( SweepContext& tcx, Node* node ) { - Fill(tcx, *node.prev); - if (node.prev->point != edge->p) { - // Next above or below edge? - if (Orient2d(*edge->q, *node.prev->point, *edge->p) == CW) { - // Below - if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) { - // Next is concave - FillLeftConcaveEdgeEvent(tcx, edge, node); - } else{ - // Next is convex - } + // if shallow stop filling + if( IsShallow( tcx, *node ) ) + { + return; } - } -} + Fill( tcx, *node ); -void Sweep::FlipEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p) -{ - Triangle& ot = t->NeighborAcross(p); - Point& op = *ot.OppositePoint(*t, p); - - if (&ot == NULL) { - // If we want to integrate the fillEdgeEvent do it here - // With current implementation we should never get here - //throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); - assert(0); - } - - if (InScanArea(p, *t->PointCCW(p), *t->PointCW(p), op)) { - // Lets rotate shared edge one vertex CW - RotateTrianglePair(*t, p, ot, op); - tcx.MapTriangleToNodes(*t); - tcx.MapTriangleToNodes(ot); - - if (p == eq && op == ep) { - if (eq == *tcx.edge_event.constrained_edge->q && ep == *tcx.edge_event.constrained_edge->p) { - t->MarkConstrainedEdge(&ep, &eq); - ot.MarkConstrainedEdge(&ep, &eq); - Legalize(tcx, *t); - Legalize(tcx, ot); - } else { - // XXX: I think one of the triangles should be legalized here? - } - } else { - Orientation o = Orient2d(eq, op, ep); - t = &NextFlipTriangle(tcx, (int)o, *t, ot, p, op); - FlipEdgeEvent(tcx, ep, eq, t, p); + if( node->prev == tcx.basin.left_node && node->next == tcx.basin.right_node ) + { + return; } - } else { - Point& newP = NextFlipPoint(ep, eq, ot, op); - FlipScanEdgeEvent(tcx, ep, eq, *t, ot, newP); - EdgeEvent(tcx, ep, eq, t, p); - } + else if( node->prev == tcx.basin.left_node ) + { + Orientation o = Orient2d( *node->point, *node->next->point, *node->next->next->point ); + + if( o == CW ) + { + return; + } + + node = node->next; + } + else if( node->next == tcx.basin.right_node ) + { + Orientation o = Orient2d( *node->point, *node->prev->point, *node->prev->prev->point ); + + if( o == CCW ) + { + return; + } + + node = node->prev; + } + else + { + // Continue with the neighbor node with lowest Y value + if( node->prev->point->y < node->next->point->y ) + { + node = node->prev; + } + else + { + node = node->next; + } + } + + FillBasinReq( tcx, node ); } -Triangle& Sweep::NextFlipTriangle(SweepContext& tcx, int o, Triangle& t, Triangle& ot, Point& p, Point& op) + +bool Sweep::IsShallow( SweepContext& tcx, Node& node ) { - if (o == CCW) { - // ot is not crossing edge after flip - int edge_index = ot.EdgeIndex(&p, &op); - ot.delaunay_edge[edge_index] = true; - Legalize(tcx, ot); - ot.ClearDelunayEdges(); - return t; - } + double height; - // t is not crossing edge after flip - int edge_index = t.EdgeIndex(&p, &op); + if( tcx.basin.left_highest ) + { + height = tcx.basin.left_node->point->y - node.point->y; + } + else + { + height = tcx.basin.right_node->point->y - node.point->y; + } - t.delaunay_edge[edge_index] = true; - Legalize(tcx, t); - t.ClearDelunayEdges(); - return ot; + // if shallow stop filling + if( tcx.basin.width > height ) + { + return true; + } + + return false; } -Point& Sweep::NextFlipPoint(Point& ep, Point& eq, Triangle& ot, Point& op) -{ - Orientation o2d = Orient2d(eq, op, ep); - if (o2d == CW) { - // Right - return *ot.PointCCW(op); - } else if (o2d == CCW) { - // Left - return *ot.PointCW(op); - } - //throw new RuntimeException("[Unsupported] Opposing point on constrained edge"); - assert(0); +void Sweep::FillEdgeEvent( SweepContext& tcx, Edge* edge, Node* node ) +{ + if( tcx.edge_event.right ) + { + FillRightAboveEdgeEvent( tcx, edge, node ); + } + else + { + FillLeftAboveEdgeEvent( tcx, edge, node ); + } +} + + +void Sweep::FillRightAboveEdgeEvent( SweepContext& tcx, Edge* edge, Node* node ) +{ + while( node->next->point->x < edge->p->x ) + { + // Check if next node is below the edge + if( Orient2d( *edge->q, *node->next->point, *edge->p ) == CCW ) + { + FillRightBelowEdgeEvent( tcx, edge, *node ); + } + else + { + node = node->next; + } + } +} + + +void Sweep::FillRightBelowEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + if( node.point->x < edge->p->x ) + { + if( Orient2d( *node.point, *node.next->point, *node.next->next->point ) == CCW ) + { + // Concave + FillRightConcaveEdgeEvent( tcx, edge, node ); + } + else + { + // Convex + FillRightConvexEdgeEvent( tcx, edge, node ); + // Retry this one + FillRightBelowEdgeEvent( tcx, edge, node ); + } + } +} + + +void Sweep::FillRightConcaveEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + Fill( tcx, *node.next ); + + if( node.next->point != edge->p ) + { + // Next above or below edge? + if( Orient2d( *edge->q, *node.next->point, *edge->p ) == CCW ) + { + // Below + if( Orient2d( *node.point, *node.next->point, *node.next->next->point ) == CCW ) + { + // Next is concave + FillRightConcaveEdgeEvent( tcx, edge, node ); + } + else + { + // Next is convex + } + } + } +} + + +void Sweep::FillRightConvexEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + // Next concave or convex? + if( Orient2d( *node.next->point, *node.next->next->point, + *node.next->next->next->point ) == CCW ) + { + // Concave + FillRightConcaveEdgeEvent( tcx, edge, *node.next ); + } + else + { + // Convex + // Next above or below edge? + if( Orient2d( *edge->q, *node.next->next->point, *edge->p ) == CCW ) + { + // Below + FillRightConvexEdgeEvent( tcx, edge, *node.next ); + } + else + { + // Above + } + } +} + + +void Sweep::FillLeftAboveEdgeEvent( SweepContext& tcx, Edge* edge, Node* node ) +{ + while( node->prev->point->x > edge->p->x ) + { + // Check if next node is below the edge + if( Orient2d( *edge->q, *node->prev->point, *edge->p ) == CW ) + { + FillLeftBelowEdgeEvent( tcx, edge, *node ); + } + else + { + node = node->prev; + } + } +} + + +void Sweep::FillLeftBelowEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + if( node.point->x > edge->p->x ) + { + if( Orient2d( *node.point, *node.prev->point, *node.prev->prev->point ) == CW ) + { + // Concave + FillLeftConcaveEdgeEvent( tcx, edge, node ); + } + else + { + // Convex + FillLeftConvexEdgeEvent( tcx, edge, node ); + // Retry this one + FillLeftBelowEdgeEvent( tcx, edge, node ); + } + } +} + + +void Sweep::FillLeftConvexEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + // Next concave or convex? + if( Orient2d( *node.prev->point, *node.prev->prev->point, + *node.prev->prev->prev->point ) == CW ) + { + // Concave + FillLeftConcaveEdgeEvent( tcx, edge, *node.prev ); + } + else + { + // Convex + // Next above or below edge? + if( Orient2d( *edge->q, *node.prev->prev->point, *edge->p ) == CW ) + { + // Below + FillLeftConvexEdgeEvent( tcx, edge, *node.prev ); + } + else + { + // Above + } + } +} + + +void Sweep::FillLeftConcaveEdgeEvent( SweepContext& tcx, Edge* edge, Node& node ) +{ + Fill( tcx, *node.prev ); + + if( node.prev->point != edge->p ) + { + // Next above or below edge? + if( Orient2d( *edge->q, *node.prev->point, *edge->p ) == CW ) + { + // Below + if( Orient2d( *node.point, *node.prev->point, *node.prev->prev->point ) == CW ) + { + // Next is concave + FillLeftConcaveEdgeEvent( tcx, edge, node ); + } + else + { + // Next is convex + } + } + } +} + + +void Sweep::FlipEdgeEvent( SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p ) +{ + Triangle& ot = t->NeighborAcross( p ); + Point& op = *ot.OppositePoint( *t, p ); + + if( &ot == NULL ) + { + // If we want to integrate the fillEdgeEvent do it here + // With current implementation we should never get here + // throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); + assert( 0 ); + } + + if( InScanArea( p, *t->PointCCW( p ), *t->PointCW( p ), op ) ) + { + // Lets rotate shared edge one vertex CW + RotateTrianglePair( *t, p, ot, op ); + tcx.MapTriangleToNodes( *t ); + tcx.MapTriangleToNodes( ot ); + + if( p == eq && op == ep ) + { + if( eq == *tcx.edge_event.constrained_edge->q + && ep == *tcx.edge_event.constrained_edge->p ) + { + t->MarkConstrainedEdge( &ep, &eq ); + ot.MarkConstrainedEdge( &ep, &eq ); + Legalize( tcx, *t ); + Legalize( tcx, ot ); + } + else + { + // XXX: I think one of the triangles should be legalized here? + } + } + else + { + Orientation o = Orient2d( eq, op, ep ); + t = &NextFlipTriangle( tcx, (int) o, *t, ot, p, op ); + FlipEdgeEvent( tcx, ep, eq, t, p ); + } + } + else + { + Point& newP = NextFlipPoint( ep, eq, ot, op ); + FlipScanEdgeEvent( tcx, ep, eq, *t, ot, newP ); + EdgeEvent( tcx, ep, eq, t, p ); + } +} + + +Triangle& Sweep::NextFlipTriangle( SweepContext& tcx, + int o, + Triangle& t, + Triangle& ot, + Point& p, + Point& op ) +{ + if( o == CCW ) + { + // ot is not crossing edge after flip + int edge_index = ot.EdgeIndex( &p, &op ); + ot.delaunay_edge[edge_index] = true; + Legalize( tcx, ot ); + ot.ClearDelunayEdges(); + return t; + } + + // t is not crossing edge after flip + int edge_index = t.EdgeIndex( &p, &op ); + + t.delaunay_edge[edge_index] = true; + Legalize( tcx, t ); + t.ClearDelunayEdges(); + return ot; +} + + +Point& Sweep::NextFlipPoint( Point& ep, Point& eq, Triangle& ot, Point& op ) +{ + Orientation o2d = Orient2d( eq, op, ep ); + + if( o2d == CW ) + { + // Right + return *ot.PointCCW( op ); + } + else if( o2d == CCW ) + { + // Left + return *ot.PointCW( op ); + } + + // throw new RuntimeException("[Unsupported] Opposing point on constrained edge"); + assert( 0 ); // Never executed, due tu assert( 0 ). Just to avoid compil warning return ep; } -void Sweep::FlipScanEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle, - Triangle& t, Point& p) + +void Sweep::FlipScanEdgeEvent( SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle, + Triangle& t, Point& p ) { - Triangle& ot = t.NeighborAcross(p); - Point& op = *ot.OppositePoint(t, p); + Triangle& ot = t.NeighborAcross( p ); + Point& op = *ot.OppositePoint( t, p ); - if (&t.NeighborAcross(p) == NULL) { - // If we want to integrate the fillEdgeEvent do it here - // With current implementation we should never get here - //throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); - assert(0); - } + if( &t.NeighborAcross( p ) == NULL ) + { + // If we want to integrate the fillEdgeEvent do it here + // With current implementation we should never get here + // throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); + assert( 0 ); + } - if (InScanArea(eq, *flip_triangle.PointCCW(eq), *flip_triangle.PointCW(eq), op)) { - // flip with new edge op->eq - FlipEdgeEvent(tcx, eq, op, &ot, op); - // TODO: Actually I just figured out that it should be possible to - // improve this by getting the next ot and op before the the above - // flip and continue the flipScanEdgeEvent here - // set new ot and op here and loop back to inScanArea test - // also need to set a new flip_triangle first - // Turns out at first glance that this is somewhat complicated - // so it will have to wait. - } else{ - Point& newP = NextFlipPoint(ep, eq, ot, op); - FlipScanEdgeEvent(tcx, ep, eq, flip_triangle, ot, newP); - } + if( InScanArea( eq, *flip_triangle.PointCCW( eq ), *flip_triangle.PointCW( eq ), op ) ) + { + // flip with new edge op->eq + FlipEdgeEvent( tcx, eq, op, &ot, op ); + // TODO: Actually I just figured out that it should be possible to + // improve this by getting the next ot and op before the the above + // flip and continue the flipScanEdgeEvent here + // set new ot and op here and loop back to inScanArea test + // also need to set a new flip_triangle first + // Turns out at first glance that this is somewhat complicated + // so it will have to wait. + } + else + { + Point& newP = NextFlipPoint( ep, eq, ot, op ); + FlipScanEdgeEvent( tcx, ep, eq, flip_triangle, ot, newP ); + } } -Sweep::~Sweep() { +Sweep::~Sweep() +{ // Clean up memory for( unsigned i = 0; i < nodes_.size(); i++ ) { delete nodes_[i]; } - } - } - diff --git a/polygon/poly2tri/sweep/sweep_context.cc b/polygon/poly2tri/sweep/sweep_context.cc index 6c0b0447dd..012421560d 100644 --- a/polygon/poly2tri/sweep/sweep_context.cc +++ b/polygon/poly2tri/sweep/sweep_context.cc @@ -33,164 +33,193 @@ #include "advancing_front.h" namespace p2t { - -SweepContext::SweepContext(std::vector polyline) : - front_(0), - head_(0), - tail_(0), - af_head_(0), - af_middle_(0), - af_tail_(0) +SweepContext::SweepContext( std::vector polyline ) : + front_( 0 ), + head_( 0 ), + tail_( 0 ), + af_head_( 0 ), + af_middle_( 0 ), + af_tail_( 0 ) { - basin = Basin(); - edge_event = EdgeEvent(); + basin = Basin(); + edge_event = EdgeEvent(); - points_ = polyline; + points_ = polyline; - InitEdges(points_); + InitEdges( points_ ); } -void SweepContext::AddHole(std::vector polyline) + +void SweepContext::AddHole( std::vector polyline ) { - InitEdges(polyline); - for(unsigned int i = 0; i < polyline.size(); i++) { - points_.push_back(polyline[i]); - } + InitEdges( polyline ); + + for( unsigned int i = 0; i < polyline.size(); i++ ) + { + points_.push_back( polyline[i] ); + } } -void SweepContext::AddPoint(Point* point) { - points_.push_back(point); + +void SweepContext::AddPoint( Point* point ) +{ + points_.push_back( point ); } + std::vector SweepContext::GetTriangles() { - return triangles_; + return triangles_; } + std::list SweepContext::GetMap() { - return map_; + return map_; } + void SweepContext::InitTriangulation() { - double xmax(points_[0]->x), xmin(points_[0]->x); - double ymax(points_[0]->y), ymin(points_[0]->y); + double xmax( points_[0]->x ), xmin( points_[0]->x ); + double ymax( points_[0]->y ), ymin( points_[0]->y ); - // Calculate bounds. - for (unsigned int i = 0; i < points_.size(); i++) { - Point& p = *points_[i]; - if (p.x > xmax) - xmax = p.x; - if (p.x < xmin) - xmin = p.x; - if (p.y > ymax) - ymax = p.y; - if (p.y < ymin) - ymin = p.y; - } + // Calculate bounds. + for( unsigned int i = 0; i < points_.size(); i++ ) + { + Point& p = *points_[i]; - double dx = kAlpha * (xmax - xmin); - double dy = kAlpha * (ymax - ymin); - head_ = new Point(xmax + dx, ymin - dy); - tail_ = new Point(xmin - dx, ymin - dy); + if( p.x > xmax ) + xmax = p.x; - // Sort points along y-axis - std::sort(points_.begin(), points_.end(), cmp); + if( p.x < xmin ) + xmin = p.x; -} + if( p.y > ymax ) + ymax = p.y; -void SweepContext::InitEdges(std::vector polyline) -{ - int num_points = polyline.size(); - for (int i = 0; i < num_points; i++) { - int j = i < num_points - 1 ? i + 1 : 0; - edge_list.push_back(new Edge(*polyline[i], *polyline[j])); - } -} - -Point* SweepContext::GetPoint(const int& index) -{ - return points_[index]; -} - -void SweepContext::AddToMap(Triangle* triangle) -{ - map_.push_back(triangle); -} - -Node& SweepContext::LocateNode(Point& point) -{ - // TODO implement search tree - return *front_->LocateNode(point.x); -} - -void SweepContext::CreateAdvancingFront(std::vector nodes) -{ - - (void) nodes; - // Initial triangle - Triangle* triangle = new Triangle(*points_[0], *tail_, *head_); - - map_.push_back(triangle); - - af_head_ = new Node(*triangle->GetPoint(1), *triangle); - af_middle_ = new Node(*triangle->GetPoint(0), *triangle); - af_tail_ = new Node(*triangle->GetPoint(2)); - front_ = new AdvancingFront(*af_head_, *af_tail_); - - // TODO: More intuitive if head is middles next and not previous? - // so swap head and tail - af_head_->next = af_middle_; - af_middle_->next = af_tail_; - af_middle_->prev = af_head_; - af_tail_->prev = af_middle_; -} - -void SweepContext::RemoveNode(Node* node) -{ - delete node; -} - -void SweepContext::MapTriangleToNodes(Triangle& t) -{ - for (int i = 0; i < 3; i++) { - if (!t.GetNeighbor(i)) { - Node* n = front_->LocatePoint(t.PointCW(*t.GetPoint(i))); - if (n) - n->triangle = &t; + if( p.y < ymin ) + ymin = p.y; } - } + + double dx = kAlpha * (xmax - xmin); + double dy = kAlpha * (ymax - ymin); + head_ = new Point( xmax + dx, ymin - dy ); + tail_ = new Point( xmin - dx, ymin - dy ); + + // Sort points along y-axis + std::sort( points_.begin(), points_.end(), cmp ); } -void SweepContext::RemoveFromMap(Triangle* triangle) + +void SweepContext::InitEdges( std::vector polyline ) { - map_.remove(triangle); -} + int num_points = polyline.size(); -void SweepContext::MeshClean(Triangle& triangle) -{ - std::vector triangles; - triangles.push_back(&triangle); - - while(!triangles.empty()){ - Triangle *t = triangles.back(); - triangles.pop_back(); - - if (t != NULL && !t->IsInterior()) { - t->IsInterior(true); - triangles_.push_back(t); - for (int i = 0; i < 3; i++) { - if (!t->constrained_edge[i]) - triangles.push_back(t->GetNeighbor(i)); - } + for( int i = 0; i < num_points; i++ ) + { + int j = i < num_points - 1 ? i + 1 : 0; + edge_list.push_back( new Edge( *polyline[i], *polyline[j] ) ); } - } } + +Point* SweepContext::GetPoint( const int& index ) +{ + return points_[index]; +} + + +void SweepContext::AddToMap( Triangle* triangle ) +{ + map_.push_back( triangle ); +} + + +Node& SweepContext::LocateNode( Point& point ) +{ + // TODO implement search tree + return *front_->LocateNode( point.x ); +} + + +void SweepContext::CreateAdvancingFront( std::vector nodes ) +{ + (void) nodes; + // Initial triangle + Triangle* triangle = new Triangle( *points_[0], *tail_, *head_ ); + + map_.push_back( triangle ); + + af_head_ = new Node( *triangle->GetPoint( 1 ), *triangle ); + af_middle_ = new Node( *triangle->GetPoint( 0 ), *triangle ); + af_tail_ = new Node( *triangle->GetPoint( 2 ) ); + front_ = new AdvancingFront( *af_head_, *af_tail_ ); + + // TODO: More intuitive if head is middles next and not previous? + // so swap head and tail + af_head_->next = af_middle_; + af_middle_->next = af_tail_; + af_middle_->prev = af_head_; + af_tail_->prev = af_middle_; +} + + +void SweepContext::RemoveNode( Node* node ) +{ + delete node; +} + + +void SweepContext::MapTriangleToNodes( Triangle& t ) +{ + for( int i = 0; i < 3; i++ ) + { + if( !t.GetNeighbor( i ) ) + { + Node* n = front_->LocatePoint( t.PointCW( *t.GetPoint( i ) ) ); + + if( n ) + n->triangle = &t; + } + } +} + + +void SweepContext::RemoveFromMap( Triangle* triangle ) +{ + map_.remove( triangle ); +} + + +void SweepContext::MeshClean( Triangle& triangle ) +{ + std::vector triangles; + + triangles.push_back( &triangle ); + + while( !triangles.empty() ) + { + Triangle* t = triangles.back(); + triangles.pop_back(); + + if( t != NULL && !t->IsInterior() ) + { + t->IsInterior( true ); + triangles_.push_back( t ); + + for( int i = 0; i < 3; i++ ) + { + if( !t->constrained_edge[i] ) + triangles.push_back( t->GetNeighbor( i ) ); + } + } + } +} + + SweepContext::~SweepContext() { - // Clean up memory delete head_; @@ -202,15 +231,15 @@ SweepContext::~SweepContext() typedef std::list type_list; - for(type_list::iterator iter = map_.begin(); iter != map_.end(); ++iter) { + for( type_list::iterator iter = map_.begin(); iter != map_.end(); ++iter ) + { Triangle* ptr = *iter; delete ptr; } - for(unsigned int i = 0; i < edge_list.size(); i++) { + for( unsigned int i = 0; i < edge_list.size(); i++ ) + { delete edge_list[i]; } - } - } diff --git a/polygon/poly2tri/sweep/sweep_context.h b/polygon/poly2tri/sweep/sweep_context.h index 1010c0e8a8..241d18a729 100644 --- a/polygon/poly2tri/sweep/sweep_context.h +++ b/polygon/poly2tri/sweep/sweep_context.h @@ -37,7 +37,6 @@ #include namespace p2t { - // Inital triangle factor, seed triangle will extend 30% of // PointSet width to both left and right. const double kAlpha = 0.3; @@ -48,139 +47,147 @@ struct Node; struct Edge; class AdvancingFront; -class SweepContext { +class SweepContext +{ public: /// Constructor -SweepContext(std::vector polyline); + SweepContext( std::vector polyline ); /// Destructor -~SweepContext(); + ~SweepContext(); -void set_head(Point* p1); + void set_head( Point* p1 ); -Point* head(); + Point* head(); -void set_tail(Point* p1); + void set_tail( Point* p1 ); -Point* tail(); + Point* tail(); -int point_count(); + int point_count(); -Node& LocateNode(Point& point); + Node& LocateNode( Point& point ); -void RemoveNode(Node* node); + void RemoveNode( Node* node ); -void CreateAdvancingFront(std::vector nodes); + void CreateAdvancingFront( std::vector nodes ); /// Try to map a node to all sides of this triangle that don't have a neighbor -void MapTriangleToNodes(Triangle& t); + void MapTriangleToNodes( Triangle& t ); -void AddToMap(Triangle* triangle); + void AddToMap( Triangle* triangle ); -Point* GetPoint(const int& index); + Point* GetPoint( const int& index ); -Point* GetPoints(); + Point* GetPoints(); -void RemoveFromMap(Triangle* triangle); + void RemoveFromMap( Triangle* triangle ); -void AddHole(std::vector polyline); + void AddHole( std::vector polyline ); -void AddPoint(Point* point); + void AddPoint( Point* point ); -AdvancingFront* front(); + AdvancingFront* front(); -void MeshClean(Triangle& triangle); + void MeshClean( Triangle& triangle ); -std::vector GetTriangles(); -std::list GetMap(); + std::vector GetTriangles(); -std::vector edge_list; + std::list GetMap(); -struct Basin { - Node* left_node; - Node* bottom_node; - Node* right_node; - double width; - bool left_highest; + std::vector edge_list; - Basin() : left_node(NULL), bottom_node(NULL), right_node(NULL), width(0.0), left_highest(false) - { - } + struct Basin + { + Node* left_node; + Node* bottom_node; + Node* right_node; + double width; + bool left_highest; - void Clear() - { - left_node = NULL; - bottom_node = NULL; - right_node = NULL; - width = 0.0; - left_highest = false; - } -}; + Basin() : left_node( NULL ), bottom_node( NULL ), right_node( NULL ), width( 0.0 ), + left_highest( false ) + { + } -struct EdgeEvent { - Edge* constrained_edge; - bool right; + void Clear() + { + left_node = NULL; + bottom_node = NULL; + right_node = NULL; + width = 0.0; + left_highest = false; + } + }; - EdgeEvent() : constrained_edge(NULL), right(false) - { - } -}; + struct EdgeEvent + { + Edge* constrained_edge; + bool right; -Basin basin; -EdgeEvent edge_event; + EdgeEvent() : constrained_edge( NULL ), right( false ) + { + } + }; + + Basin basin; + EdgeEvent edge_event; private: -friend class Sweep; + friend class Sweep; -std::vector triangles_; -std::list map_; -std::vector points_; + std::vector triangles_; + std::list map_; + std::vector points_; // Advancing front -AdvancingFront* front_; + AdvancingFront* front_; // head point used with advancing front -Point* head_; + Point* head_; // tail point used with advancing front -Point* tail_; + Point* tail_; -Node *af_head_, *af_middle_, *af_tail_; - -void InitTriangulation(); -void InitEdges(std::vector polyline); + Node* af_head_, * af_middle_, * af_tail_; + void InitTriangulation(); + void InitEdges( std::vector polyline ); }; inline AdvancingFront* SweepContext::front() { - return front_; + return front_; } + inline int SweepContext::point_count() { - return points_.size(); + return points_.size(); } -inline void SweepContext::set_head(Point* p1) + +inline void SweepContext::set_head( Point* p1 ) { - head_ = p1; + head_ = p1; } + inline Point* SweepContext::head() { - return head_; + return head_; } -inline void SweepContext::set_tail(Point* p1) + +inline void SweepContext::set_tail( Point* p1 ) { - tail_ = p1; + tail_ = p1; } + inline Point* SweepContext::tail() { - return tail_; + return tail_; } - } #endif