Merge pull request #3 from saurabhb17/master

secondary files

1 files changed, 739 insertions, 0 deletions

diff --git a/common/geometry/hetriang.cpp b/common/geometry/hetriang.cpp
new file mode 100644
index 0000000..e5cf26c
--- /dev/null
+++ b/common/geometry/hetriang.cpp
@@ -0,0 +1,739 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * <http://www.gnu.org/licenses/>.
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#include <ttl/halfedge/hetriang.h>
+#include <ttl/halfedge/hetraits.h>
+#include <ttl/ttl.h>
+#include <algorithm>
+#include <fstream>
+#include <limits>
+#include <boost/foreach.hpp>
+#include <boost/make_shared.hpp>
+#include <class_board_connected_item.h>
+
+using namespace hed;
+
+#ifdef TTL_USE_NODE_ID
+  int NODE::id_count = 0;
+#endif
+
+
+void NODE::updateLayers()
+{
+    assert( m_layers.none() );
+
+    BOOST_FOREACH( const BOARD_CONNECTED_ITEM* item, m_parents )
+        m_layers |= item->GetLayerSet();
+}
+
+
+//#define DEBUG_HE
+#ifdef DEBUG_HE
+#include <iostream>
+static void errorAndExit( char* aMessage )
+{
+    cout << "\n!!! ERROR: "<< aMessage << " !!!\n" << endl;
+    exit( -1 );
+}
+#endif
+
+
+static EDGE_PTR getLeadingEdgeInTriangle( const EDGE_PTR& aEdge )
+{
+    EDGE_PTR edge = aEdge;
+
+    // Code: 3EF (assumes triangle)
+    if( !edge->IsLeadingEdge() )
+    {
+        edge = edge->GetNextEdgeInFace();
+
+        if( !edge->IsLeadingEdge() )
+            edge = edge->GetNextEdgeInFace();
+    }
+
+    if( !edge->IsLeadingEdge() )
+    {
+        return EDGE_PTR();
+    }
+
+    return edge;
+}
+
+
+static void getLimits( NODES_CONTAINER::iterator aFirst, NODES_CONTAINER::iterator aLast,
+                       int& aXmin, int& aYmin, int& aXmax, int& aYmax)
+{
+    aXmin = aYmin = std::numeric_limits<int>::min();
+    aXmax = aYmax = std::numeric_limits<int>::max();
+
+    NODES_CONTAINER::iterator it;
+
+    for( it = aFirst; it != aLast; ++it )
+    {
+        aXmin = std::min( aXmin, ( *it )->GetX() );
+        aYmin = std::min( aYmin, ( *it )->GetY() );
+        aXmax = std::max( aXmax, ( *it )->GetX() );
+        aYmax = std::max( aYmax, ( *it )->GetY() );
+    }
+}
+
+
+EDGE_PTR TRIANGULATION::InitTwoEnclosingTriangles( NODES_CONTAINER::iterator aFirst,
+                                                  NODES_CONTAINER::iterator aLast)
+{
+    int xmin, ymin, xmax, ymax;
+    getLimits( aFirst, aLast, xmin, ymin, xmax, ymax );
+
+    // Add 10% of range:
+    double fac = 10.0;
+    double dx = ( xmax - xmin ) / fac;
+    double dy = ( ymax - ymin ) / fac;
+
+    NODE_PTR n1 = boost::make_shared<NODE>( xmin - dx, ymin - dy );
+    NODE_PTR n2 = boost::make_shared<NODE>( xmax + dx, ymin - dy );
+    NODE_PTR n3 = boost::make_shared<NODE>( xmax + dx, ymax + dy );
+    NODE_PTR n4 = boost::make_shared<NODE>( xmin - dx, ymax + dy );
+
+    // diagonal
+    EDGE_PTR e1d = boost::make_shared<EDGE>();
+    EDGE_PTR e2d = boost::make_shared<EDGE>();
+
+    // lower triangle
+    EDGE_PTR e11 = boost::make_shared<EDGE>();
+    EDGE_PTR e12 = boost::make_shared<EDGE>();
+
+    // upper triangle
+    EDGE_PTR e21 = boost::make_shared<EDGE>();
+    EDGE_PTR e22 = boost::make_shared<EDGE>();
+
+    // lower triangle
+    e1d->SetSourceNode( n3 );
+    e1d->SetNextEdgeInFace( e11 );
+    e1d->SetTwinEdge( e2d );
+    addLeadingEdge( e1d );
+
+    e11->SetSourceNode( n1 );
+    e11->SetNextEdgeInFace( e12 );
+
+    e12->SetSourceNode( n2 );
+    e12->SetNextEdgeInFace( e1d );
+
+    // upper triangle
+    e2d->SetSourceNode( n1 );
+    e2d->SetNextEdgeInFace( e21 );
+    e2d->SetTwinEdge( e1d );
+    addLeadingEdge( e2d );
+
+    e21->SetSourceNode( n3 );
+    e21->SetNextEdgeInFace( e22 );
+
+    e22->SetSourceNode( n4 );
+    e22->SetNextEdgeInFace( e2d );
+
+    return e11;
+}
+
+
+TRIANGULATION::TRIANGULATION()
+{
+    m_helper = new ttl::TRIANGULATION_HELPER( *this );
+}
+
+
+TRIANGULATION::TRIANGULATION( const TRIANGULATION& aTriangulation )
+{
+    m_helper = 0;   // make coverity and static analysers quiet.
+    // Triangulation: Copy constructor not present
+    assert( false );
+}
+
+
+TRIANGULATION::~TRIANGULATION()
+{
+    cleanAll();
+    delete m_helper;
+}
+
+
+void TRIANGULATION::CreateDelaunay( NODES_CONTAINER::iterator aFirst,
+                                    NODES_CONTAINER::iterator aLast )
+{
+    cleanAll();
+
+    EDGE_PTR bedge = InitTwoEnclosingTriangles( aFirst, aLast );
+    DART dc( bedge );
+
+    DART d_iter = dc;
+
+    NODES_CONTAINER::iterator it;
+    for( it = aFirst; it != aLast; ++it )
+    {
+        m_helper->InsertNode<TTLtraits>( d_iter, *it );
+    }
+
+    // In general (e.g. for the triangle based data structure), the initial dart
+    // may have been changed.
+    // It is the users responsibility to get a valid boundary dart here.
+    // The half-edge data structure preserves the initial dart.
+    // (A dart at the boundary can also be found by trying to locate a
+    // triangle "outside" the triangulation.)
+
+    // Assumes rectangular domain
+    m_helper->RemoveRectangularBoundary<TTLtraits>( dc );
+}
+
+
+void TRIANGULATION::RemoveTriangle( EDGE_PTR& aEdge )
+{
+  EDGE_PTR e1 = getLeadingEdgeInTriangle( aEdge );
+
+#ifdef DEBUG_HE
+    if( !e1 )
+    errorAndExit( "Triangulation::removeTriangle: could not find leading aEdge" );
+#endif
+
+    removeLeadingEdgeFromList( e1 );
+    // cout << "No leading edges = " << leadingEdges_.size() << endl;
+    // Remove the triangle
+    EDGE_PTR e2( e1->GetNextEdgeInFace() );
+    EDGE_PTR e3( e2->GetNextEdgeInFace() );
+
+    e1->Clear();
+    e2->Clear();
+    e3->Clear();
+}
+
+
+void TRIANGULATION::ReverseSplitTriangle( EDGE_PTR& aEdge )
+{
+    // Reverse operation of splitTriangle
+    EDGE_PTR e1( aEdge->GetNextEdgeInFace() );
+    EDGE_PTR le( getLeadingEdgeInTriangle( e1 ) );
+#ifdef DEBUG_HE
+    if (!le)
+    errorAndExit("Triangulation::removeTriangle: could not find leading edge");
+#endif
+    removeLeadingEdgeFromList( le );
+
+    EDGE_PTR e2( e1->GetNextEdgeInFace()->GetTwinEdge()->GetNextEdgeInFace() );
+    le = getLeadingEdgeInTriangle( e2 );
+#ifdef DEBUG_HE
+    if (!le)
+    errorAndExit("Triangulation::removeTriangle: could not find leading edge");
+#endif
+    removeLeadingEdgeFromList( le );
+
+    EDGE_PTR e3( aEdge->GetTwinEdge()->GetNextEdgeInFace()->GetNextEdgeInFace() );
+    le = getLeadingEdgeInTriangle( e3 );
+#ifdef DEBUG_HE
+    if (!le)
+    errorAndExit("Triangulation::removeTriangle: could not find leading edge");
+#endif
+    removeLeadingEdgeFromList( le );
+
+    // The three triangles at the node have now been removed
+    // from the triangulation, but the arcs have not been deleted.
+    // Next delete the 6 half edges radiating from the node
+    // The node is maintained by handle and need not be deleted explicitly
+    EDGE_PTR estar = aEdge;
+    EDGE_PTR enext = estar->GetTwinEdge()->GetNextEdgeInFace();
+    estar->GetTwinEdge()->Clear();
+    estar->Clear();
+
+    estar = enext;
+    enext = estar->GetTwinEdge()->GetNextEdgeInFace();
+    estar->GetTwinEdge()->Clear();
+    estar->Clear();
+
+    enext->GetTwinEdge()->Clear();
+    enext->Clear();
+
+    // Create the new triangle
+    e1->SetNextEdgeInFace( e2 );
+    e2->SetNextEdgeInFace( e3 );
+    e3->SetNextEdgeInFace( e1 );
+    addLeadingEdge( e1 );
+}
+
+
+DART TRIANGULATION::CreateDart()
+{
+  // Return an arbitrary CCW dart
+    return DART( *m_leadingEdges.begin() );
+}
+
+
+bool TRIANGULATION::removeLeadingEdgeFromList( EDGE_PTR& aLeadingEdge )
+{
+    // Remove the edge from the list of leading edges,
+    // but don't delete it.
+    // Also set flag for leading edge to false.
+    // Must search from start of list. Since edges are added to the
+    // start of the list during triangulation, this operation will
+    // normally be fast (when used in the triangulation algorithm)
+    std::list<EDGE_PTR>::iterator it;
+    for( it = m_leadingEdges.begin(); it != m_leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        if( edge == aLeadingEdge )
+        {
+            edge->SetAsLeadingEdge( false );
+            it = m_leadingEdges.erase( it );
+
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+void TRIANGULATION::cleanAll()
+{
+    BOOST_FOREACH( EDGE_PTR& edge, m_leadingEdges )
+        edge->SetNextEdgeInFace( EDGE_PTR() );
+}
+
+
+void TRIANGULATION::swapEdge( DART& aDart )
+{
+    SwapEdge( aDart.GetEdge() );
+}
+
+
+void TRIANGULATION::splitTriangle( DART& aDart, const NODE_PTR& aPoint )
+{
+    EDGE_PTR edge = SplitTriangle( aDart.GetEdge(), aPoint );
+    aDart.Init( edge );
+}
+
+
+void TRIANGULATION::reverseSplitTriangle( DART& aDart )
+{
+    ReverseSplitTriangle( aDart.GetEdge() );
+}
+
+
+void TRIANGULATION::removeBoundaryTriangle( DART& aDart )
+{
+    RemoveTriangle( aDart.GetEdge() );
+}
+
+
+#ifdef TTL_USE_NODE_FLAG
+void TRIANGULATION::FlagNodes( bool aFlag ) const
+{
+    std::list<EDGE_PTR>::const_iterator it;
+    for( it = m_leadingEdges.begin(); it != m_leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        for( int i = 0; i < 3; ++i )
+        {
+            edge->GetSourceNode()->SetFlag( aFlag );
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+}
+
+
+std::list<NODE_PTR>* TRIANGULATION::GetNodes() const
+{
+    FlagNodes( false );
+    std::list<NODE_PTR>* nodeList = new std::list<NODE_PTR>;
+    std::list<EDGE_PTR>::const_iterator it;
+
+    for( it = m_leadingEdges.begin(); it != m_leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        for( int i = 0; i < 3; ++i )
+        {
+            const NODE_PTR& node = edge->GetSourceNode();
+
+            if( node->GetFlag() == false )
+            {
+                nodeList->push_back( node );
+                node->SetFlag( true );
+            }
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+    return nodeList;
+}
+#endif
+
+
+std::list<EDGE_PTR>* TRIANGULATION::GetEdges( bool aSkipBoundaryEdges ) const
+{
+    // collect all arcs (one half edge for each arc)
+    // (boundary edges are also collected).
+    std::list<EDGE_PTR>::const_iterator it;
+    std::list<EDGE_PTR>* elist = new std::list<EDGE_PTR>;
+
+    for( it = m_leadingEdges.begin(); it != m_leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+        for( int i = 0; i < 3; ++i )
+        {
+            EDGE_PTR twinedge = edge->GetTwinEdge();
+            // only one of the half-edges
+
+            if( ( !twinedge && !aSkipBoundaryEdges )
+                    || ( twinedge && ( (size_t) edge.get() > (size_t) twinedge.get() ) ) )
+                elist->push_front( edge );
+
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+
+    return elist;
+}
+
+
+EDGE_PTR TRIANGULATION::SplitTriangle( EDGE_PTR& aEdge, const NODE_PTR& aPoint )
+{
+    // Add a node by just splitting a triangle into three triangles
+    // Assumes the half aEdge is located in the triangle
+    // Returns a half aEdge with source node as the new node
+
+    // e#_n are new edges
+    // e# are existing edges
+    // e#_n and e##_n are new twin edges
+    // e##_n are edges incident to the new node
+
+    // Add the node to the structure
+    //NODE_PTR new_node(new Node(x,y,z));
+
+    NODE_PTR n1( aEdge->GetSourceNode() );
+    EDGE_PTR e1( aEdge );
+
+    EDGE_PTR e2( aEdge->GetNextEdgeInFace() );
+    NODE_PTR n2( e2->GetSourceNode() );
+
+    EDGE_PTR e3( e2->GetNextEdgeInFace() );
+    NODE_PTR n3( e3->GetSourceNode() );
+
+    EDGE_PTR e1_n = boost::make_shared<EDGE>();
+    EDGE_PTR e11_n = boost::make_shared<EDGE>();
+    EDGE_PTR e2_n = boost::make_shared<EDGE>();
+    EDGE_PTR e22_n = boost::make_shared<EDGE>();
+    EDGE_PTR e3_n = boost::make_shared<EDGE>();
+    EDGE_PTR e33_n = boost::make_shared<EDGE>();
+
+    e1_n->SetSourceNode( n1 );
+    e11_n->SetSourceNode( aPoint );
+    e2_n->SetSourceNode( n2 );
+    e22_n->SetSourceNode( aPoint );
+    e3_n->SetSourceNode( n3 );
+    e33_n->SetSourceNode( aPoint );
+
+    e1_n->SetTwinEdge( e11_n );
+    e11_n->SetTwinEdge( e1_n );
+    e2_n->SetTwinEdge( e22_n );
+    e22_n->SetTwinEdge( e2_n );
+    e3_n->SetTwinEdge( e33_n );
+    e33_n->SetTwinEdge( e3_n );
+
+    e1_n->SetNextEdgeInFace( e33_n );
+    e2_n->SetNextEdgeInFace( e11_n );
+    e3_n->SetNextEdgeInFace( e22_n );
+
+    e11_n->SetNextEdgeInFace( e1 );
+    e22_n->SetNextEdgeInFace( e2 );
+    e33_n->SetNextEdgeInFace( e3 );
+
+    // and update old's next aEdge
+    e1->SetNextEdgeInFace( e2_n );
+    e2->SetNextEdgeInFace( e3_n );
+    e3->SetNextEdgeInFace( e1_n );
+
+    // add the three new leading edges,
+    // Must remove the old leading aEdge from the list.
+    // Use the field telling if an aEdge is a leading aEdge
+    // NOTE: Must search in the list!!!
+
+    if( e1->IsLeadingEdge() )
+        removeLeadingEdgeFromList( e1 );
+    else if( e2->IsLeadingEdge() )
+        removeLeadingEdgeFromList( e2 );
+    else if( e3->IsLeadingEdge() )
+        removeLeadingEdgeFromList( e3 );
+    else
+        assert( false );        // one of the edges should be leading
+
+    addLeadingEdge( e1_n );
+    addLeadingEdge( e2_n );
+    addLeadingEdge( e3_n );
+
+    // Return a half aEdge incident to the new node (with the new node as source node)
+
+    return e11_n;
+}
+
+
+void TRIANGULATION::SwapEdge( EDGE_PTR& aDiagonal )
+{
+    // Note that diagonal is both input and output and it is always
+    // kept in counterclockwise direction (this is not required by all
+    // functions in TriangulationHelper now)
+
+    // Swap by rotating counterclockwise
+    // Use the same objects - no deletion or new objects
+    EDGE_PTR eL( aDiagonal );
+    EDGE_PTR eR( eL->GetTwinEdge() );
+    EDGE_PTR eL_1( eL->GetNextEdgeInFace() );
+    EDGE_PTR eL_2( eL_1->GetNextEdgeInFace() );
+    EDGE_PTR eR_1( eR->GetNextEdgeInFace() );
+    EDGE_PTR eR_2( eR_1->GetNextEdgeInFace() );
+
+    // avoid node to be dereferenced to zero and deleted
+    NODE_PTR nR( eR_2->GetSourceNode() );
+    NODE_PTR nL( eL_2->GetSourceNode() );
+
+    eL->SetSourceNode( nR );
+    eR->SetSourceNode( nL );
+
+    // and now 6 1-sewings
+    eL->SetNextEdgeInFace( eL_2 );
+    eL_2->SetNextEdgeInFace( eR_1 );
+    eR_1->SetNextEdgeInFace( eL );
+
+    eR->SetNextEdgeInFace( eR_2 );
+    eR_2->SetNextEdgeInFace( eL_1 );
+    eL_1->SetNextEdgeInFace( eR );
+
+    if( eL->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eL );
+    else if( eL_1->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eL_1 );
+    else if( eL_2->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eL_2 );
+
+    if( eR->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eR );
+    else if( eR_1->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eR_1 );
+    else if( eR_2->IsLeadingEdge() )
+        removeLeadingEdgeFromList( eR_2 );
+
+    addLeadingEdge( eL );
+    addLeadingEdge( eR );
+}
+
+
+bool TRIANGULATION::CheckDelaunay() const
+{
+    // ???? outputs !!!!
+    // ofstream os("qweND.dat");
+    const std::list<EDGE_PTR>& leadingEdges = GetLeadingEdges();
+
+    std::list<EDGE_PTR>::const_iterator it;
+    bool ok = true;
+    int noNotDelaunay = 0;
+
+    for( it = leadingEdges.begin(); it != leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        for( int i = 0; i < 3; ++i )
+        {
+            EDGE_PTR twinedge = edge->GetTwinEdge();
+
+            // only one of the half-edges
+            if( !twinedge || (size_t) edge.get() > (size_t) twinedge.get() )
+            {
+                DART dart( edge );
+                if( m_helper->SwapTestDelaunay<TTLtraits>( dart ) )
+                {
+                    noNotDelaunay++;
+
+                    //printEdge(dart,os); os << "\n";
+                    ok = false;
+                    //cout << "............. not Delaunay .... " << endl;
+                }
+            }
+
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+
+#ifdef DEBUG_HE
+    cout << "!!! Triangulation is NOT Delaunay: " << noNotDelaunay << " edges\n" << endl;
+#endif
+
+    return ok;
+}
+
+
+void TRIANGULATION::OptimizeDelaunay()
+{
+    // This function is also present in ttl where it is implemented
+    // generically.
+    // The implementation below is tailored for the half-edge data structure,
+    // and is thus more efficient
+
+    // Collect all interior edges (one half edge for each arc)
+    bool skip_boundary_edges = true;
+    std::list<EDGE_PTR>* elist = GetEdges( skip_boundary_edges );
+
+    // Assumes that elist has only one half-edge for each arc.
+    bool cycling_check = true;
+    bool optimal = false;
+    std::list<EDGE_PTR>::const_iterator it;
+
+    while( !optimal )
+    {
+        optimal = true;
+
+        for( it = elist->begin(); it != elist->end(); ++it )
+        {
+            EDGE_PTR edge = *it;
+
+            DART dart( edge );
+            // Constrained edges should not be swapped
+            if( m_helper->SwapTestDelaunay<TTLtraits>( dart, cycling_check ) )
+            {
+                optimal = false;
+                SwapEdge( edge );
+            }
+        }
+    }
+
+    delete elist;
+}
+
+
+EDGE_PTR TRIANGULATION::GetInteriorNode() const
+{
+    const std::list<EDGE_PTR>& leadingEdges = GetLeadingEdges();
+    std::list<EDGE_PTR>::const_iterator it;
+
+    for( it = leadingEdges.begin(); it != leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        // multiple checks, but only until found
+        for( int i = 0; i < 3; ++i )
+        {
+            if( edge->GetTwinEdge() )
+            {
+                if( !m_helper->IsBoundaryNode( DART( edge ) ) )
+                    return edge;
+            }
+
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+
+    return EDGE_PTR(); // no boundary nodes
+}
+
+
+EDGE_PTR TRIANGULATION::GetBoundaryEdgeInTriangle( const EDGE_PTR& aEdge ) const
+{
+    EDGE_PTR edge = aEdge;
+
+    if( m_helper->IsBoundaryEdge( DART( edge ) ) )
+        return edge;
+
+    edge = edge->GetNextEdgeInFace();
+    if( m_helper->IsBoundaryEdge( DART( edge ) ) )
+        return edge;
+
+    edge = edge->GetNextEdgeInFace();
+    if( m_helper->IsBoundaryEdge( DART( edge ) ) )
+        return edge;
+
+    return EDGE_PTR();
+}
+
+
+EDGE_PTR TRIANGULATION::GetBoundaryEdge() const
+{
+    // Get an arbitrary (CCW) boundary edge
+    // If the triangulation is closed, NULL is returned
+    const std::list<EDGE_PTR>& leadingEdges = GetLeadingEdges();
+    std::list<EDGE_PTR>::const_iterator it;
+    EDGE_PTR edge;
+
+    for( it = leadingEdges.begin(); it != leadingEdges.end(); ++it )
+    {
+        edge = GetBoundaryEdgeInTriangle( *it );
+
+        if( edge )
+            return edge;
+    }
+    return EDGE_PTR();
+}
+
+
+void TRIANGULATION::PrintEdges( std::ofstream& aOutput ) const
+{
+    // Print source node and target node for each edge face by face,
+    // but only one of the half-edges.
+    const std::list<EDGE_PTR>& leadingEdges = GetLeadingEdges();
+    std::list<EDGE_PTR>::const_iterator it;
+
+    for( it = leadingEdges.begin(); it != leadingEdges.end(); ++it )
+    {
+        EDGE_PTR edge = *it;
+
+        for( int i = 0; i < 3; ++i )
+        {
+            EDGE_PTR twinedge = edge->GetTwinEdge();
+
+            // Print only one edge (the highest value of the pointer)
+            if( !twinedge || (size_t) edge.get() > (size_t) twinedge.get() )
+            {
+                // Print source node and target node
+                NODE_PTR node = edge->GetSourceNode();
+                aOutput << node->GetX() << " " << node->GetY() << std::endl;
+                node = edge->GetTargetNode();
+                aOutput << node->GetX() << " " << node->GetY() << std::endl;
+                aOutput << '\n'; // blank line
+            }
+
+            edge = edge->GetNextEdgeInFace();
+        }
+    }
+}