0023105: Exception during Meshing / Missing triangles
[occt.git] / src / BRepMesh / BRepMesh_Delaun.cxx
index 4761f64..36ac549 100755 (executable)
 // and conditions governing the rights and limitations under the License.
 
 
-#include <BRepMesh_Delaun.ixx>
+#include <BRepMesh_Delaun.hxx>
+
+#include <gp.hxx>
 #include <gp_XY.hxx>
 #include <gp_Pnt2d.hxx>
 #include <gp_Vec2d.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
-#include <TColStd_ListOfInteger.hxx>
+
 #include <Precision.hxx>
 #include <Bnd_Box2d.hxx>
-#include <gp.hxx>
+#include <Bnd_B2d.hxx>
+
 #include <TColStd_MapOfInteger.hxx>
 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
 #include <TColStd_Array1OfBoolean.hxx>
+#include <TColStd_ListIteratorOfListOfInteger.hxx>
+#include <TColStd_ListOfInteger.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_Array1OfInteger.hxx>
+#include <TColStd_SequenceOfInteger.hxx>
+
 #include <BRepMesh_MapOfIntegerInteger.hxx>
 #include <BRepMesh_HeapSortIndexedVertexOfDelaun.hxx>
 #include <BRepMesh_ComparatorOfIndexedVertexOfDelaun.hxx>
 #include <BRepMesh_SelectorOfDataStructureOfDelaun.hxx>
 #include <BRepMesh_HeapSortVertexOfDelaun.hxx>
 #include <BRepMesh_ComparatorOfVertexOfDelaun.hxx>
-#include <TColgp_Array1OfXY.hxx>
-#include <TColStd_Array1OfReal.hxx>
+#include <BRepMesh_Array1OfVertexOfDelaun.hxx>
+
+#include <BRepMesh_Edge.hxx>
+#include <BRepMesh_Vertex.hxx>
+#include <BRepMesh_Triangle.hxx>
+
+#include <NCollection_Vector.hxx>
 
 typedef TColStd_ListIteratorOfListOfInteger  IteratorOnListOfInteger;
 typedef TColStd_ListOfInteger                ListOfInteger;
 
-const Standard_Real EPSEPS = Precision::PConfusion() * Precision::PConfusion();
-const gp_XY SortingDirection(M_SQRT1_2, M_SQRT1_2);
+const Standard_Real AngDeviation1Deg  = M_PI/180.;
+const Standard_Real AngDeviation90Deg = 90 * AngDeviation1Deg;
+const Standard_Real Angle2PI          = 2 * M_PI;
 
-//=======================================================================
-//function : fillBndBox
-//purpose  : Add boundig box for edge defined by start & end point to
-//           the given vector of bounding boxes for triangulation edges
-//=======================================================================
-static void fillBndBox( NCollection_Vector <Bnd_Box2d>& theVB,
-                        const BRepMesh_Vertex& theV1,
-                        const BRepMesh_Vertex& theV2 )
-{
-  Bnd_Box2d aBox;      
-  aBox.Add( gp_Pnt2d( theV1.Coord() ) );
-  aBox.Add( gp_Pnt2d( theV2.Coord() ) );
-  theVB.Append( aBox );
-}
+const Standard_Real Precision    = Precision::PConfusion();
+const Standard_Real EndPrecision = 1 - Precision;
+const Standard_Real Precision2   = Precision * Precision;
+const gp_XY SortingDirection(M_SQRT1_2, M_SQRT1_2);
 
 //=======================================================================
 //function : BRepMesh_Delaun
 //purpose  : Creates the triangulation with an empty Mesh data structure
 //=======================================================================
 BRepMesh_Delaun::BRepMesh_Delaun( BRepMesh_Array1OfVertexOfDelaun& theVertices,
-                                  const Standard_Boolean isPositive )
-: myPositiveOrientation( isPositive ),
+                                  const Standard_Boolean           isPositive )
+: myIsPositiveOrientation( isPositive ),
   myCircles( theVertices.Length(), new NCollection_IncAllocator() )
 {
   if ( theVertices.Length() > 2 )
@@ -83,10 +87,10 @@ BRepMesh_Delaun::BRepMesh_Delaun( BRepMesh_Array1OfVertexOfDelaun& theVertices,
 //function : BRepMesh_Delaun
 //purpose  : Creates the triangulation with and existent Mesh data structure
 //=======================================================================
-BRepMesh_Delaun::BRepMesh_Delaun ( const Handle( BRepMesh_DataStructureOfDelaun )& theOldMesh, 
-                                   BRepMesh_Array1OfVertexOfDelaun& theVertices,
-                                   const Standard_Boolean isPositive )
- : myPositiveOrientation( isPositive ),
+BRepMesh_Delaun::BRepMesh_Delaun( const Handle( BRepMesh_DataStructureOfDelaun )& theOldMesh,
+                                  BRepMesh_Array1OfVertexOfDelaun&                theVertices,
+                                  const Standard_Boolean                          isPositive )
+ : myIsPositiveOrientation( isPositive ),
    myCircles( theVertices.Length(), theOldMesh->Allocator() )
 {
   myMeshData = theOldMesh;
@@ -99,21 +103,21 @@ BRepMesh_Delaun::BRepMesh_Delaun ( const Handle( BRepMesh_DataStructureOfDelaun
 //purpose  : Creates the triangulation with and existent Mesh data structure
 //=======================================================================
 BRepMesh_Delaun::BRepMesh_Delaun( const Handle( BRepMesh_DataStructureOfDelaun )& theOldMesh, 
-                                  TColStd_Array1OfInteger& theVertexIndexes,
-                                  const Standard_Boolean isPositive )
- : myPositiveOrientation( isPositive ),
-   myCircles( theVertexIndexes.Length(), theOldMesh->Allocator() )
+                                  TColStd_Array1OfInteger&                        theVertexIndices,
+                                  const Standard_Boolean                          isPositive )
+ : myIsPositiveOrientation( isPositive ),
+   myCircles( theVertexIndices.Length(), theOldMesh->Allocator() )
 {
   myMeshData = theOldMesh;
-  if ( theVertexIndexes.Length() > 2 )
+  if ( theVertexIndices.Length() > 2 )
   {
     Bnd_Box2d aBox;
-    Standard_Integer anIndex = theVertexIndexes.Lower();
-    Standard_Integer anUpper = theVertexIndexes.Upper();
+    Standard_Integer anIndex = theVertexIndices.Lower();
+    Standard_Integer anUpper = theVertexIndices.Upper();
     for ( ; anIndex <= anUpper; ++anIndex )
-      aBox.Add( gp_Pnt2d( GetVertex( theVertexIndexes( anIndex) ).Coord() ) );
+      aBox.Add( gp_Pnt2d( GetVertex( theVertexIndices( anIndex) ).Coord() ) );
 
-    Perform( aBox, theVertexIndexes );
+    perform( aBox, theVertexIndices );
   }
 }
 
@@ -135,31 +139,31 @@ void BRepMesh_Delaun::Init( BRepMesh_Array1OfVertexOfDelaun& theVertices )
     aVertexIndexes( anIndex ) = myMeshData->AddNode( theVertices( anIndex ) );
   }
 
-  Perform( aBox, aVertexIndexes );
+  perform( aBox, aVertexIndexes );
 }
 
 //=======================================================================
-//function : Perform
+//function : perform
 //purpose  : Create super mesh and run triangulation procedure
 //=======================================================================
-void BRepMesh_Delaun::Perform( Bnd_Box2d& theBndBox,
+void BRepMesh_Delaun::perform( Bnd_Box2d&               theBndBox,
                                TColStd_Array1OfInteger& theVertexIndexes )
 {
-  theBndBox.Enlarge( Precision::PConfusion() );
-  SuperMesh( theBndBox );
+  theBndBox.Enlarge( Precision );
+  superMesh( theBndBox );
 
   BRepMesh_HeapSortIndexedVertexOfDelaun::Sort( theVertexIndexes, 
       BRepMesh_ComparatorOfIndexedVertexOfDelaun( SortingDirection,
-        Precision::PConfusion(), myMeshData ) );
+        Precision, myMeshData ) );
 
-  Compute( theVertexIndexes );
+  compute( theVertexIndexes );
 }
 
 //=======================================================================
-//function : SuperMesh
+//function : superMesh
 //purpose  : Build the super mesh
 //=======================================================================
-void BRepMesh_Delaun::SuperMesh( const Bnd_Box2d& theBox )
+void BRepMesh_Delaun::superMesh( const Bnd_Box2d& theBox )
 {
   Standard_Real aMinX, aMinY, aMaxX, aMaxY;
   theBox.Get( aMinX, aMinY, aMaxX, aMaxY );
@@ -181,41 +185,46 @@ void BRepMesh_Delaun::SuperMesh( const Bnd_Box2d& theBox )
   myCircles.SetCellSize( aDeltaX / aScaler,
                          aDeltaY / aScaler );
 
-  mySupVert1 = myMeshData->AddNode(
+  mySupVert[0] = myMeshData->AddNode(
     BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) );
     
-  mySupVert2 = myMeshData->AddNode(
+  mySupVert[1] = myMeshData->AddNode(
     BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
     
-  mySupVert3 = myMeshData->AddNode(
+  mySupVert[2] = myMeshData->AddNode(
     BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
 
-  if ( !myPositiveOrientation )
+  if ( !myIsPositiveOrientation )
   {
     Standard_Integer aTmp;
-    aTmp       = mySupVert2;
-    mySupVert2 = mySupVert3;
-    mySupVert3 = aTmp;
+    aTmp         = mySupVert[1];
+    mySupVert[1] = mySupVert[2];
+    mySupVert[2] = aTmp;
   }
 
-  Standard_Integer anEdge1, anEdge2, anEdge3;
+  Standard_Integer anEdgeId[3];
   
-  anEdge1 = myMeshData->AddLink( BRepMesh_Edge( mySupVert1, mySupVert2, BRepMesh_Free ) );
-  anEdge2 = myMeshData->AddLink( BRepMesh_Edge( mySupVert2, mySupVert3, BRepMesh_Free ) );
-  anEdge3 = myMeshData->AddLink( BRepMesh_Edge( mySupVert3, mySupVert1, BRepMesh_Free ) );
+  for (Standard_Integer aNodeId = 0; aNodeId < 3; ++aNodeId)
+  {
+    Standard_Integer aFirstNode = aNodeId;
+    Standard_Integer aLastNode  = (aNodeId + 1) % 3;
+    anEdgeId[aNodeId] = myMeshData->AddLink( BRepMesh_Edge( 
+      mySupVert[aFirstNode], mySupVert[aLastNode], BRepMesh_Free ) );
+  }
   
-  mySupTrian = BRepMesh_Triangle( Abs( anEdge1 ), Abs( anEdge2 ), Abs( anEdge3 ), 
-    ( anEdge1 > 0 ), ( anEdge2 > 0 ), ( anEdge1 > 0), BRepMesh_Free);
+  mySupTrian = BRepMesh_Triangle( 
+    Abs( anEdgeId[0] ),  Abs( anEdgeId[1] ),  Abs( anEdgeId[2] ), 
+    ( anEdgeId[0] > 0 ), ( anEdgeId[1] > 0 ), ( anEdgeId[2] > 0),
+    BRepMesh_Free);
 }
 
 //=======================================================================
-//function : DeleteTriangle
-//purpose : The concerned triangles are deleted and the freed edges are added in 
-//           <theLoopEdges>. If an edge is added twice, it does not exist and
-//          it is necessary to destroy it. This corresponds to the destruction of two
-//          connected triangles.
+//function : deleteTriangle
+//purpose  : Deletes the triangle with the given index and adds the free
+//           edges into the map.
+//           When an edge is suppressed more than one time it is destroyed.
 //=======================================================================
-void  BRepMesh_Delaun::DeleteTriangle( const Standard_Integer theIndex, 
+void  BRepMesh_Delaun::deleteTriangle( const Standard_Integer        theIndex, 
                                        BRepMesh_MapOfIntegerInteger& theLoopEdges )
 {
   myCircles.Delete( theIndex );
@@ -224,7 +233,7 @@ void  BRepMesh_Delaun::DeleteTriangle( const Standard_Integer theIndex,
   Standard_Boolean o[3];
   GetTriangle( theIndex ).Edges( e[0], e[1], e[2],
                                  o[0], o[1], o[2] );
-                                 
+  
   myMeshData->RemoveElement( theIndex );
 
   for ( Standard_Integer i = 0; i < 3; ++i )
@@ -238,11 +247,11 @@ void  BRepMesh_Delaun::DeleteTriangle( const Standard_Integer theIndex,
 }
 
 //=======================================================================
-//function : Compute
+//function : compute
 //purpose  : Computes the triangulation and add the vertices edges and 
 //           triangles to the Mesh data structure
 //=======================================================================
-void BRepMesh_Delaun::Compute( TColStd_Array1OfInteger& theVertexIndexes )
+void BRepMesh_Delaun::compute( TColStd_Array1OfInteger& theVertexIndexes )
 {
   // Insertion of edges of super triangles in the list of free edges: 
   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
@@ -259,22 +268,21 @@ void BRepMesh_Delaun::Compute( TColStd_Array1OfInteger& theVertexIndexes )
   {
     // Creation of 3 trianglers with the first node and the edges of the super triangle:
     Standard_Integer anVertexIdx = theVertexIndexes.Lower();
-    CreateTriangles( theVertexIndexes( anVertexIdx ), aLoopEdges );
+    createTriangles( theVertexIndexes( anVertexIdx ), aLoopEdges );
 
     // Add other nodes to the mesh
-    CreateTrianglesOnNewVertices( theVertexIndexes );
+    createTrianglesOnNewVertices( theVertexIndexes );
   }
 
   // Destruction of triangles containing a top of the super triangle
   BRepMesh_SelectorOfDataStructureOfDelaun aSelector( myMeshData );
-  aSelector.NeighboursOfNode( mySupVert1 );
-  aSelector.NeighboursOfNode( mySupVert2 );
-  aSelector.NeighboursOfNode( mySupVert3 );
+  for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
+    aSelector.NeighboursOfNode( mySupVert[aSupVertId] );
   
   aLoopEdges.Clear();
   BRepMesh_MapOfInteger::Iterator aFreeTriangles( aSelector.Elements() );
   for ( ; aFreeTriangles.More(); aFreeTriangles.Next() )
-    DeleteTriangle( aFreeTriangles.Key(), aLoopEdges );
+    deleteTriangle( aFreeTriangles.Key(), aLoopEdges );
 
   // All edges that remain free are removed from aLoopEdges;
   // only the boundary edges of the triangulation remain there
@@ -286,16 +294,15 @@ void BRepMesh_Delaun::Compute( TColStd_Array1OfInteger& theVertexIndexes )
   }
 
   // The tops of the super triangle are destroyed
-  myMeshData->RemoveNode( mySupVert1 );
-  myMeshData->RemoveNode( mySupVert2 );
-  myMeshData->RemoveNode( mySupVert3 );
+  for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
+    myMeshData->RemoveNode( mySupVert[aSupVertId] );
 }
 
 //=======================================================================
-//function : CreateTriangles
+//function : createTriangles
 //purpose  : Creates the triangles beetween the node and the polyline.
 //=======================================================================
-void BRepMesh_Delaun::CreateTriangles ( const Standard_Integer theVertexIndex,  
+void BRepMesh_Delaun::createTriangles ( const Standard_Integer        theVertexIndex,  
                                         BRepMesh_MapOfIntegerInteger& thePoly )
 {
   ListOfInteger aLoopEdges, anExternalEdges;
@@ -304,90 +311,85 @@ void BRepMesh_Delaun::CreateTriangles ( const Standard_Integer theVertexIndex,
   BRepMesh_MapOfIntegerInteger::Iterator anEdges( thePoly );
   for ( ; anEdges.More(); anEdges.Next() )
   {
-    const BRepMesh_Edge& anEdge = GetEdge( anEdges.Key() );
-    Standard_Integer aFirstNode = anEdge.FirstNode();
-    Standard_Integer aLastNode  = anEdge.LastNode();
-    Standard_Boolean isPositive = (Standard_Boolean)thePoly( anEdges.Key() );
-    if ( !isPositive )
+    Standard_Integer     anEdgeId = anEdges.Key();
+    const BRepMesh_Edge& anEdge   = GetEdge( anEdgeId );
+
+    Standard_Boolean isPositive = (Standard_Boolean)thePoly( anEdgeId );
+
+    Standard_Integer aNodes[3];
+    if ( isPositive )
     {
-      Standard_Integer aTmp;
-      aTmp       = aFirstNode;
-      aFirstNode = aLastNode;
-      aLastNode  = aTmp;
+      aNodes[0] = anEdge.FirstNode();
+      aNodes[2] = anEdge.LastNode();
     }
+    else
+    {
+      aNodes[0] = anEdge.LastNode();
+      aNodes[2] = anEdge.FirstNode();
+    }
+    aNodes[1] = theVertexIndex;
 
-    const BRepMesh_Vertex& aFirstVertex = GetVertex( aFirstNode );
-    const BRepMesh_Vertex& aLastVertex  = GetVertex( aLastNode  );
+    const BRepMesh_Vertex& aFirstVertex = GetVertex( aNodes[0] );
+    const BRepMesh_Vertex& aLastVertex  = GetVertex( aNodes[2] );
 
-    gp_XY aVEdge1, aVEdge2, aVEdge3;
-    aVEdge1 = aFirstVertex.Coord();
-    aVEdge1.Subtract( aVertexCoord );
-    
-    aVEdge2 = aLastVertex.Coord();
-    aVEdge2.Subtract( aFirstVertex.Coord() );
-    
-    aVEdge3 = aVertexCoord;
-    aVEdge3.Subtract( aLastVertex.Coord() );
+    gp_XY anEdgeDir( aLastVertex.Coord() - aFirstVertex.Coord() );
+    Standard_Real anEdgeLen = anEdgeDir.Modulus();
+    if ( anEdgeLen < Precision )
+      continue;
 
-    Standard_Real aDist12 = 0., aDist23 = 0.;
-    Standard_Real aV2Len  = aVEdge2.Modulus();
-    if ( aV2Len > Precision::PConfusion() )
+    anEdgeDir.SetCoord( anEdgeDir.X() / anEdgeLen,
+                        anEdgeDir.Y() / anEdgeLen );
+
+    gp_XY aFirstLinkDir( aFirstVertex.Coord() - aVertexCoord );
+    gp_XY aLastLinkDir ( aVertexCoord         - aLastVertex.Coord() );
+                      
+    Standard_Real aDist12 = aFirstLinkDir ^ anEdgeDir;
+    Standard_Real aDist23 = anEdgeDir     ^ aLastLinkDir;
+    if ( Abs( aDist12 ) < Precision || 
+         Abs( aDist23 ) < Precision )
     {
-      aVEdge2.SetCoord( aVEdge2.X() / aV2Len,
-                        aVEdge2.Y() / aV2Len );
-                        
-      aDist12 = aVEdge1 ^ aVEdge2;
-      aDist23 = aVEdge2 ^ aVEdge3;
+      continue;
     }
 
-    if ( Abs( aDist12 ) >= Precision::PConfusion() && 
-         Abs( aDist23 ) >= Precision::PConfusion() )
+    Standard_Boolean isSensOK;
+    if ( myIsPositiveOrientation )
+      isSensOK = ( aDist12 > 0. && aDist23 > 0.);
+    else
+      isSensOK = ( aDist12 < 0. && aDist23 < 0.);
+    
+
+    BRepMesh_Edge aFirstLink( aNodes[1], aNodes[0], BRepMesh_Free );
+    BRepMesh_Edge aLastLink ( aNodes[2], aNodes[1], BRepMesh_Free );
+
+    Standard_Integer anEdgesInfo[3] = {
+      myMeshData->AddLink( aFirstLink ),
+      isPositive ? anEdgeId : -anEdgeId,
+      myMeshData->AddLink( aLastLink ) };
+
+    if ( isSensOK )
     {
-      Standard_Boolean isSensOK;
-      if ( myPositiveOrientation )
-        isSensOK = ( aDist12 > 0. && aDist23 > 0.);
-      else
-        isSensOK = ( aDist12 < 0. && aDist23 < 0.);
-        
-      if ( isSensOK )
+      Standard_Integer anEdges[3];
+      Standard_Boolean anEdgesOri[3];
+      for ( Standard_Integer aTriLinkIt = 0; aTriLinkIt < 3; ++aTriLinkIt )
       {
-        BRepMesh_Edge aNewEdge1 (theVertexIndex, aFirstNode, BRepMesh_Free);
-        BRepMesh_Edge aNewEdge3 (aLastNode, theVertexIndex, BRepMesh_Free);
-
-        Standard_Integer iNewEdge1 = myMeshData->AddLink (aNewEdge1);
-        Standard_Integer iNewEdge3 = myMeshData->AddLink (aNewEdge3);
-          
-        BRepMesh_Triangle aNewTri (Abs (iNewEdge1), anEdges.Key(), Abs (iNewEdge3),
-                                   (iNewEdge1 > 0), isPositive, (iNewEdge3 > 0),
-                                   BRepMesh_Free);
-        Standard_Integer iNewTri = myMeshData->AddElement(aNewTri);
-
-        Standard_Boolean isCircleCreated = 
-          myCircles.Add (aVertexCoord, aFirstVertex.Coord(), aLastVertex.Coord(), iNewTri);
-                         
-        if ( !isCircleCreated )
-          myMeshData->RemoveElement (iNewTri);
+        const Standard_Integer& anEdgeInfo = anEdgesInfo[aTriLinkIt];
+        anEdges[aTriLinkIt]    = Abs( anEdgeInfo );
+        anEdgesOri[aTriLinkIt] = anEdgeInfo > 0;
       }
+
+      addTriangle( anEdges, anEdgesOri, aNodes );
+    }
+    else
+    {
+      if ( isPositive )
+        aLoopEdges.Append(  anEdges.Key() );
       else
-      {
-        if ( isPositive )
-          aLoopEdges.Append(  anEdges.Key() );
-        else
-          aLoopEdges.Append( -anEdges.Key() );
-          
-        if ( aVEdge1.SquareModulus() > aVEdge3.SquareModulus() )
-        {
-          BRepMesh_Edge aNewEdge (theVertexIndex, aFirstNode, BRepMesh_Free);
-          Standard_Integer iNewEdge = myMeshData->AddLink(aNewEdge);
-          anExternalEdges.Append (Abs (iNewEdge));
-        }
-        else
-        {
-          BRepMesh_Edge aNewEdge (aLastNode, theVertexIndex, BRepMesh_Free);
-          Standard_Integer iNewEdge = myMeshData->AddLink (aNewEdge);
-          anExternalEdges.Append (Abs (iNewEdge));
-        }
-      }
+        aLoopEdges.Append( -anEdges.Key() );
+        
+      if ( aFirstLinkDir.SquareModulus() > aLastLinkDir.SquareModulus() )
+        anExternalEdges.Append( Abs( anEdgesInfo[0] ) );
+      else
+        anExternalEdges.Append( Abs( anEdgesInfo[2] ) );
     }
   }
   
@@ -399,7 +401,7 @@ void BRepMesh_Delaun::CreateTriangles ( const Standard_Integer theVertexIndex,
     
     
     if ( !aPair.IsEmpty() )
-      DeleteTriangle( aPair.FirstIndex(), thePoly );
+      deleteTriangle( aPair.FirstIndex(), thePoly );
       
     anExternalEdges.RemoveFirst();
   }
@@ -416,7 +418,7 @@ void BRepMesh_Delaun::CreateTriangles ( const Standard_Integer theVertexIndex,
     if ( anEdge.Movability() != BRepMesh_Deleted )
     {
       Standard_Integer anEdgeIdx = aLoopEdges.First();
-      MeshLeftPolygonOf( Abs( anEdgeIdx ), ( anEdgeIdx > 0 ) );
+      meshLeftPolygonOf( Abs( anEdgeIdx ), ( anEdgeIdx > 0 ) );
     }
       
     aLoopEdges.RemoveFirst();
@@ -424,10 +426,10 @@ void BRepMesh_Delaun::CreateTriangles ( const Standard_Integer theVertexIndex,
 }
 
 //=======================================================================
-//function : CreateTrianglesOnNewVertices
+//function : createTrianglesOnNewVertices
 //purpose  : Creation of triangles from the new nodes
 //=======================================================================
-void BRepMesh_Delaun::CreateTrianglesOnNewVertices( TColStd_Array1OfInteger& theVertexIndexes )
+void BRepMesh_Delaun::createTrianglesOnNewVertices( TColStd_Array1OfInteger& theVertexIndexes )
 {
   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
 
@@ -472,7 +474,7 @@ void BRepMesh_Delaun::CreateTrianglesOnNewVertices( TColStd_Array1OfInteger& the
 
     if ( aTriangleId > 0 )
     {
-      DeleteTriangle( aTriangleId, aLoopEdges );
+      deleteTriangle( aTriangleId, aLoopEdges );
 
       isModify = Standard_True;    
       while ( isModify && !aCirclesList.IsEmpty() )
@@ -491,7 +493,7 @@ void BRepMesh_Delaun::CreateTrianglesOnNewVertices( TColStd_Array1OfInteger& the
                aLoopEdges.IsBound( e[2] ) )
           {
             isModify = Standard_True;
-            DeleteTriangle( aCircleIt1.Value(), aLoopEdges );
+            deleteTriangle( aCircleIt1.Value(), aLoopEdges );
             aCirclesList.Remove( aCircleIt1 );
             break;
           }
@@ -500,105 +502,166 @@ void BRepMesh_Delaun::CreateTrianglesOnNewVertices( TColStd_Array1OfInteger& the
 
       // Creation of triangles with the current node and free edges
       // and removal of these edges from the list of free edges
-      CreateTriangles( aVertexIdx, aLoopEdges );
+      createTriangles( aVertexIdx, aLoopEdges );
     }
   }
   // Check that internal edges are not crossed by triangles
-  BRepMesh_MapOfInteger::Iterator anInernalEdgesIt( InternalEdges() );
+  Handle(BRepMesh_MapOfInteger) anInternalEdges = InternalEdges();
 
   // Destruction of triancles intersecting internal edges 
   // and their replacement by makeshift triangles
-  anInernalEdgesIt.Reset();
+  BRepMesh_MapOfInteger::Iterator anInernalEdgesIt( *anInternalEdges );
   for ( ; anInernalEdgesIt.More(); anInernalEdgesIt.Next() )
   {
     Standard_Integer aNbC;
     aNbC = myMeshData->ElemConnectedTo( anInernalEdgesIt.Key() ).Extent();
     if ( aNbC == 0 )
     {
-      MeshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_True  ); 
-      MeshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_False ); 
+      meshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_True  ); 
+      meshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_False ); 
     }
   }
 
   // Adjustment of meshes to boundary edges
-  FrontierAdjust();
+  frontierAdjust();
 }
 
 //=======================================================================
-//function : CleanupMesh
-//purpose  : Cleanup mesh from the free triangles
+//function : isBoundToFrontier
+//purpose  : Goes through the neighbour triangles around the given node 
+//           started from the given link, returns TRUE if some triangle 
+//           has a bounding frontier edge or FALSE elsewhere.
+//           Stop link is used to prevent cycles.
+//           Previous element Id is used to identify next neighbor element.
 //=======================================================================
-void BRepMesh_Delaun::CleanupMesh()
+Standard_Boolean BRepMesh_Delaun::isBoundToFrontier(
+  const Standard_Integer theRefNodeId,
+  const Standard_Integer theRefLinkId,
+  const Standard_Integer theStopLinkId,
+  const Standard_Integer thePrevElementId)
 {
-  BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
-  ListOfInteger aTrianglesList;
+  const BRepMesh_PairOfIndex& aPair = 
+    myMeshData->ElemConnectedTo( theRefLinkId );
+  if ( aPair.IsEmpty() )
+    return Standard_False;
 
-  for(;;)
+  Standard_Integer aNbElements = aPair.Extent();
+  for ( Standard_Integer anElemIt = 1; anElemIt <= aNbElements; ++anElemIt )
   {
-    aTrianglesList.Clear();
-    aLoopEdges.Clear();
+    const Standard_Integer aTriId = aPair.Index(anElemIt);
+    if ( aTriId < 0 || aTriId == thePrevElementId )
+      continue;
+
+    Standard_Integer anEdges[3];
+    Standard_Boolean anEdgesOri[3];
+    GetTriangle( aTriId ).Edges( anEdges[0], anEdges[1], anEdges[2],
+      anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
+
+    for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
+    {
+      const Standard_Integer anEdgeId = anEdges[anEdgeIt];
+      if ( anEdgeId == theRefLinkId )
+        continue;
+
+      if ( anEdgeId == theStopLinkId )
+        return Standard_False;
+
+      const BRepMesh_Edge& anEdge = GetEdge( anEdgeId );
+      if ( anEdge.FirstNode() != theRefNodeId &&
+           anEdge.LastNode()  != theRefNodeId )
+      {
+        continue;
+      }
+
+      if ( anEdge.Movability() != BRepMesh_Free )
+        return Standard_True;
 
-    BRepMesh_MapOfInteger::Iterator aFreeEdgesIt( FreeEdges() );
+      return isBoundToFrontier( theRefNodeId, anEdgeId, theStopLinkId, aTriId );
+    }
+  }
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : cleanupMesh
+//purpose  : Cleanup mesh from the free triangles
+//=======================================================================
+void BRepMesh_Delaun::cleanupMesh()
+{
+  while ( Standard_True )
+  {
+    BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
+    NCollection_Map<Standard_Integer> aDelTriangles;
+
+    Handle(BRepMesh_MapOfInteger) aFreeEdges = FreeEdges();
+    BRepMesh_MapOfInteger::Iterator aFreeEdgesIt( *aFreeEdges );
     for ( ; aFreeEdgesIt.More(); aFreeEdgesIt.Next() )
     {
-      const BRepMesh_Edge& anEdge = GetEdge( aFreeEdgesIt.Key() );
-      if ( anEdge.Movability() != BRepMesh_Frontier )
+      const Standard_Integer& aFreeEdgeId = aFreeEdgesIt.Key();
+      const BRepMesh_Edge&    anEdge      = GetEdge( aFreeEdgeId );
+      if ( anEdge.Movability() == BRepMesh_Frontier )
+        continue;
+
+      const BRepMesh_PairOfIndex& aPair = 
+        myMeshData->ElemConnectedTo( aFreeEdgeId );
+      if ( aPair.IsEmpty() )
       {
-        Standard_Integer aFrontierNb = 0;
-        if ( myMeshData->ElemConnectedTo( aFreeEdgesIt.Key() ).IsEmpty() ) 
-          aLoopEdges.Bind( aFreeEdgesIt.Key(), Standard_True );
-        else
+        aLoopEdges.Bind( aFreeEdgeId, Standard_True );
+        continue;
+      }
+
+      Standard_Integer aTriId = aPair.FirstIndex();
+
+      // Check that the connected triangle is not surrounded by another triangles
+      Standard_Integer anEdges[3];
+      Standard_Boolean anEdgesOri[3];
+      GetTriangle( aTriId ).Edges( anEdges[0], anEdges[1], anEdges[2],
+        anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
+
+      Standard_Boolean isCanNotBeRemoved = Standard_True;
+      for ( Standard_Integer aCurEdgeIdx = 0; aCurEdgeIdx < 3; ++aCurEdgeIdx )
+      {
+        if ( anEdges[aCurEdgeIdx] != aFreeEdgeId )
+          continue;
+
+        for ( Standard_Integer anOtherEdgeIt = 1; anOtherEdgeIt <= 2; ++anOtherEdgeIt )
         {
-          BRepMesh_ListOfInteger::Iterator aConnectedLinksIt( 
-            myMeshData->LinkNeighboursOf( anEdge.FirstNode() ) );
-            
-          for ( ; aConnectedLinksIt.More(); aConnectedLinksIt.Next() )
-          {
-            if ( GetEdge( aConnectedLinksIt.Value() ).Movability() == BRepMesh_Frontier )
-            {
-              aFrontierNb++;
-              if ( aFrontierNb > 1 )
-                break;
-            }
-          }
-          
-          if ( aFrontierNb == 2 )
+          Standard_Integer anOtherEdgeId = ( aCurEdgeIdx + anOtherEdgeIt ) % 3;
+          const BRepMesh_PairOfIndex& anOtherEdgePair = 
+            myMeshData->ElemConnectedTo( anEdges[anOtherEdgeId] );
+
+          if ( anOtherEdgePair.Extent() < 2 )
           {
-            const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( aFreeEdgesIt.Key() );
-            for ( Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; ++j )
-            {
-              const Standard_Integer anElemId = aPair.Index(j);
-              if ( anElemId < 0 )
-                continue;
-                
-              Standard_Integer e[3];
-              Standard_Boolean o[3];
-              GetTriangle( anElemId ).Edges( e[0], e[1], e[2],
-                                             o[0], o[1], o[2] );
-              
-              Standard_Boolean isTriangleToDelete = Standard_True;
-              for ( Standard_Integer k = 0; k < 3; ++k )
-              {
-                if ( GetEdge( e[k] ).Movability() != BRepMesh_Free )
-                {
-                  isTriangleToDelete = Standard_False;
-                  break;
-                }
-              }
-
-              if ( isTriangleToDelete )
-                aTrianglesList.Append( anElemId );
-            }
+            isCanNotBeRemoved = Standard_False;
+            break;
           }
         }
+
+        break;
+      }
+
+      if ( isCanNotBeRemoved )
+        continue;
+
+      Standard_Boolean isConnected[2] = { Standard_False, Standard_False };
+      for ( Standard_Integer aLinkNodeIt = 0; aLinkNodeIt < 2; ++aLinkNodeIt )
+      {
+        isConnected[aLinkNodeIt] = isBoundToFrontier( ( aLinkNodeIt == 0 ) ? 
+          anEdge.FirstNode() : anEdge.LastNode(), 
+          aFreeEdgeId, aFreeEdgeId, -1);
       }
+
+      if ( !isConnected[0] || !isConnected[1] )
+        aDelTriangles.Add( aTriId );
     }
 
     // Destruction of triangles :
     Standard_Integer aDeletedTrianglesNb = 0;
-    for ( ; !aTrianglesList.IsEmpty(); aTrianglesList.RemoveFirst() )
+    NCollection_Map<Standard_Integer>::Iterator aDelTrianglesIt( aDelTriangles );
+    for ( ; aDelTrianglesIt.More(); aDelTrianglesIt.Next() )
     {
-      DeleteTriangle( aTrianglesList.First(), aLoopEdges );
+      deleteTriangle( aDelTrianglesIt.Key(), aLoopEdges );
       aDeletedTrianglesNb++;
     }
 
@@ -616,736 +679,1261 @@ void BRepMesh_Delaun::CleanupMesh()
 }
 
 //=======================================================================
-//function : FrontierAdjust
+//function : frontierAdjust
 //purpose  : Adjust the mesh on the frontier
 //=======================================================================
-void BRepMesh_Delaun::FrontierAdjust()
+void BRepMesh_Delaun::frontierAdjust()
 {
-  BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
-  ListOfInteger aTrianglesList;
-
+  Handle(BRepMesh_MapOfInteger)         aFrontier = Frontier();
+  NCollection_Vector<Standard_Integer>  aFailedFrontiers;
+  BRepMesh_MapOfIntegerInteger          aLoopEdges( 10, myMeshData->Allocator() );
+  BRepMesh_Delaun::HandleOfMapOfInteger aIntFrontierEdges = new NCollection_Map<Standard_Integer>();
   for ( Standard_Integer aPass = 1; aPass <= 2; ++aPass )
   {      
-    // 1 pass): find external triangles on boundary edges
-    // 2 pass): find external triangles on boundary edges
-    // because in comlex curved boundaries external triangles can appear  
-    // after "mesh left aPolygonon"
+    // 1 pass): find external triangles on boundary edges;
+    // 2 pass): find external triangles on boundary edges appeared 
+    //          during triangles replacement.
     
-    BRepMesh_MapOfInteger::Iterator aFrontierIt( Frontier() );
+    BRepMesh_MapOfInteger::Iterator aFrontierIt( *aFrontier );
     for ( ; aFrontierIt.More(); aFrontierIt.Next() )
     {
-      const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( aFrontierIt.Key() );
-      for( Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; ++j )
+      Standard_Integer aFrontierId = aFrontierIt.Key();
+      const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( aFrontierId );
+      Standard_Integer aNbElem = aPair.Extent();
+      for( Standard_Integer aElemIt = 1; aElemIt <= aNbElem; ++aElemIt )
       {
-        const Standard_Integer aPriorElemId = aPair.Index(j);
+        const Standard_Integer aPriorElemId = aPair.Index( aElemIt );
         if( aPriorElemId < 0 )
-            continue;
+          continue;
             
         Standard_Integer e[3];
         Standard_Boolean o[3];
         GetTriangle( aPriorElemId ).Edges( e[0], e[1], e[2],
                                            o[0], o[1], o[2] );
 
+        Standard_Boolean isTriangleFound = Standard_False;
         for ( Standard_Integer n = 0; n < 3; ++n )
         {
-          if ( aFrontierIt.Key() == e[n] && !o[n] )
+          if ( aFrontierId == e[n] && !o[n] )
           {
-            aTrianglesList.Append( aPriorElemId );
+            // Destruction  of external triangles on boundary edges
+            isTriangleFound = Standard_True;
+            deleteTriangle( aPriorElemId, aLoopEdges );
             break;
           }
         }
+
+        if ( isTriangleFound )
+          break;
       }
     }
 
-    // destruction  of external triangles on boundary edges
-    for ( ; !aTrianglesList.IsEmpty(); aTrianglesList.RemoveFirst() )
-      DeleteTriangle( aTrianglesList.First(), aLoopEdges );
-
     // destrucrion of remaining hanging edges :
     BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
     for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
     {
-      if (myMeshData->ElemConnectedTo( aLoopEdgesIt.Key() ).IsEmpty() )
-        myMeshData->RemoveLink( aLoopEdgesIt.Key() );
+      Standard_Integer aLoopEdgeId = aLoopEdgesIt.Key();
+      if (myMeshData->ElemConnectedTo( aLoopEdgeId ).IsEmpty() )
+        myMeshData->RemoveLink( aLoopEdgeId );
     }
-      
 
     // destruction of triangles crossing the boundary edges and 
     // their replacement by makeshift triangles
-    if ( aPass == 1 )
-      aFrontierIt.Reset();
-    else
-      // in some cases there remain unused boundaries check
-      aFrontierIt.Initialize( Frontier() );
-
-    NCollection_Vector<Bnd_Box2d> aFrontBoxes;
-    for ( ; aFrontierIt.More(); aFrontierIt.Next() )
+    for ( aFrontierIt.Reset(); aFrontierIt.More(); aFrontierIt.Next() )
     {
-      const BRepMesh_Edge& anEdge = GetEdge( aFrontierIt.Key() );
-      const BRepMesh_Vertex& aV1  = GetVertex( anEdge.FirstNode() );
-      const BRepMesh_Vertex& aV2  = GetVertex( anEdge.LastNode()  );
-      fillBndBox( aFrontBoxes, aV1, aV2 );
+      Standard_Integer aFrontierId = aFrontierIt.Key();
+      if ( !myMeshData->ElemConnectedTo( aFrontierId ).IsEmpty() )
+        continue;
 
-      if ( myMeshData->ElemConnectedTo( aFrontierIt.Key() ).IsEmpty() )
-        MeshLeftPolygonOf( aFrontierIt.Key(), Standard_True ); 
+      Standard_Boolean isSuccess = 
+        meshLeftPolygonOf( aFrontierId, Standard_True, aIntFrontierEdges );
+
+      if ( aPass == 2 && !isSuccess )
+        aFailedFrontiers.Append( aFrontierId );
     }
+  }
 
-    if ( aPass == 1 ) 
+  cleanupMesh();
+
+  // When the mesh has been cleaned up, try to process frontier edges 
+  // once again to fill the possible gaps that might be occured in case of "saw" -
+  // situation when frontier edge has a triangle at a right side, but its free 
+  // links cross another frontieres  and meshLeftPolygonOf itself can't collect 
+  // a closed polygon.
+  NCollection_Vector<Standard_Integer>::Iterator aFailedFrontiersIt( aFailedFrontiers );
+  for ( ; aFailedFrontiersIt.More(); aFailedFrontiersIt.Next() )
+  {
+    Standard_Integer aFrontierId = aFailedFrontiersIt.Value();
+    if ( !myMeshData->ElemConnectedTo( aFrontierId ).IsEmpty() )
       continue;
 
-    CleanupMesh();
+    meshLeftPolygonOf( aFrontierId, Standard_True, aIntFrontierEdges );
   }
 }
 
 //=======================================================================
-//function : KillInternalTriangles
-//purpose  : Removes triangles within polygon
+//function : fillBndBox
+//purpose  : Add boundig box for edge defined by start & end point to
+//           the given vector of bounding boxes for triangulation edges
+//=======================================================================
+void BRepMesh_Delaun::fillBndBox( NCollection_Sequence<Bnd_B2d>& theBoxes,
+                                  const BRepMesh_Vertex&         theV1,
+                                  const BRepMesh_Vertex&         theV2 )
+{
+  Bnd_B2d aBox;      
+  aBox.Add( theV1.Coord() );
+  aBox.Add( theV2.Coord() );
+  theBoxes.Append( aBox );
+}
+
+//=======================================================================
+//function : meshLeftPolygonOf
+//purpose  : Collect the polygon at the left of the given edge (material side)
 //=======================================================================
-void BRepMesh_Delaun::KillInternalTriangles( Standard_Integer theEdgeId, 
-                                             const TColStd_MapOfInteger& theIgnoredEdges,
-                                             BRepMesh_MapOfIntegerInteger& theLoopEdges )
+Standard_Boolean BRepMesh_Delaun::meshLeftPolygonOf( 
+  const Standard_Integer                theStartEdgeId,
+  const Standard_Boolean                isForward,
+  BRepMesh_Delaun::HandleOfMapOfInteger theSkipped )
 {
-  if ( theIgnoredEdges.Contains( theEdgeId ) )
-    return;   
+  if ( !theSkipped.IsNull() && theSkipped->Contains( theStartEdgeId ) )
+    return Standard_True;
+
+  const BRepMesh_Edge& aRefEdge = GetEdge( theStartEdgeId );
 
-  const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( theEdgeId );
-  for ( Standard_Integer i = 1; i <= aPair.Extent(); ++i )
+  TColStd_SequenceOfInteger aPolygon;
+  Standard_Integer aStartNode, aPivotNode;
+  if ( isForward )
   {
-    Standard_Integer anElemId = aPair.Index(i);
-    if( anElemId < 0 )
-      continue;
+    aPolygon.Append( theStartEdgeId );
+    aStartNode = aRefEdge.FirstNode();
+    aPivotNode = aRefEdge.LastNode();
+  }
+  else
+  {
+    aPolygon.Append( -theStartEdgeId );
+    aStartNode = aRefEdge.LastNode();
+    aPivotNode = aRefEdge.FirstNode();
+  }
+
+
+  const BRepMesh_Vertex& aStartEdgeVertexS = GetVertex( aStartNode );
+  BRepMesh_Vertex        aPivotVertex      = GetVertex( aPivotNode );
+
+  gp_Vec2d aRefLinkDir( aPivotVertex.Coord() - 
+    aStartEdgeVertexS.Coord() );
+
+  if ( aRefLinkDir.SquareMagnitude() < Precision2 )
+    return Standard_True;
+
+  // Auxilary structures.
+  // Bounding boxes of polygon links to be used for preliminary
+  // analysis of intersections
+  NCollection_Sequence<Bnd_B2d> aBoxes;
+  fillBndBox( aBoxes, aStartEdgeVertexS, aPivotVertex );
+
+  // Hanging ends
+  NCollection_Map<Standard_Integer> aDeadLinks;
+
+  // Links are temporarily excluded from consideration
+  NCollection_Map<Standard_Integer> aLeprousLinks;
+  aLeprousLinks.Add( theStartEdgeId );
+
+  Standard_Boolean isSkipLeprous = Standard_True;
+  Standard_Integer aFirstNode    = aStartNode;
+  while ( aPivotNode != aFirstNode )
+  {
+    Bnd_B2d          aNextLinkBndBox;
+    gp_Vec2d         aNextLinkDir;
+    Standard_Integer aNextPivotNode = 0;
+
+    Standard_Integer aNextLinkId = findNextPolygonLink(
+      aFirstNode,
+      aPivotNode,     aPivotVertex,  aRefLinkDir, 
+      aBoxes,         aPolygon,      theSkipped,
+      isSkipLeprous,  aLeprousLinks, aDeadLinks, 
+      aNextPivotNode, aNextLinkDir,  aNextLinkBndBox );
+
+    if ( aNextLinkId != 0 )
+    {
+      aStartNode   = aPivotNode;
+      aRefLinkDir  = aNextLinkDir;
+
+      aPivotNode   = aNextPivotNode;
+      aPivotVertex = GetVertex( aNextPivotNode );
 
-    Standard_Integer e[3];
-    Standard_Boolean o[3];
-    GetTriangle( anElemId ).Edges( e[0], e[1], e[2],
-                                   o[0], o[1], o[2] );
+      aBoxes.Append  ( aNextLinkBndBox );
+      aPolygon.Append( aNextLinkId );
 
-    for ( Standard_Integer anIndex = 0; anIndex < 3; ++anIndex )
+      isSkipLeprous = Standard_True;
+    }
+    else
     {
-      if ( e[anIndex] == theEdgeId )
+      // Nothing to do
+      if ( aPolygon.Length() == 1 )
+        return Standard_False;
+
+      // Return to the previous point
+      Standard_Integer aDeadLinkId = Abs( aPolygon.Last() );
+      aDeadLinks.Add      ( aDeadLinkId );
+
+      aLeprousLinks.Remove( aDeadLinkId );
+      aPolygon.Remove     ( aPolygon.Length() );
+      aBoxes.Remove       ( aBoxes.Length() );
+
+      Standard_Integer aPrevLinkInfo = aPolygon.Last();
+      const BRepMesh_Edge& aPrevLink = GetEdge( Abs( aPrevLinkInfo ) );
+
+      if( aPrevLinkInfo > 0 )
       {
-        DeleteTriangle( anElemId, theLoopEdges );
-        for ( Standard_Integer n = 0; n < 2; ++n )
-        {
-          Standard_Integer aCurEdge = e[(anIndex + n + 1) % 3];
-          KillInternalTriangles( aCurEdge, theIgnoredEdges, theLoopEdges );
-        }
+        aStartNode = aPrevLink.FirstNode();
+        aPivotNode = aPrevLink.LastNode();
+      }
+      else
+      {
+        aStartNode = aPrevLink.LastNode();
+        aPivotNode = aPrevLink.FirstNode();
       }
+      
+      aPivotVertex = GetVertex( aPivotNode );
+      aRefLinkDir = 
+        aPivotVertex.Coord() - GetVertex( aStartNode ).Coord();
+
+      isSkipLeprous = Standard_False;
     }
   }
+
+  if ( aPolygon.Length() < 3 )
+    return Standard_False;
+
+  cleanupPolygon( aPolygon, aBoxes );
+  meshPolygon   ( aPolygon, aBoxes, theSkipped );
+
+  return Standard_True;
 }
 
 //=======================================================================
-//function : RemovePivotTriangles
-//purpose  : Removes triangles around the given pivot node
+//function : findNextPolygonLink
+//purpose  : Find next link starting from the given node and has maximum
+//           angle respect the given reference link.
+//           Each time the next link is found other neighbor links at the 
+//           pivot node are marked as leprous and will be excluded from 
+//           consideration next time until a hanging end is occured.
 //=======================================================================
-void BRepMesh_Delaun::RemovePivotTriangles( const Standard_Integer theEdgeInfo,
-                                            const Standard_Integer thePivotNode,
-                                            TColStd_MapOfInteger& theInfectedEdges,
-                                            BRepMesh_MapOfIntegerInteger& theLoopEdges,
-                                            const Standard_Boolean isFirstPass )
+Standard_Integer BRepMesh_Delaun::findNextPolygonLink(
+  const Standard_Integer&                     theFirstNode,
+  const Standard_Integer&                     thePivotNode,
+  const BRepMesh_Vertex&                      thePivotVertex,
+  const gp_Vec2d&                             theRefLinkDir,
+  const NCollection_Sequence<Bnd_B2d>&        theBoxes,
+  const TColStd_SequenceOfInteger&            thePolygon,
+  const BRepMesh_Delaun::HandleOfMapOfInteger theSkipped,
+  const Standard_Boolean&                     isSkipLeprous,
+  NCollection_Map<Standard_Integer>&          theLeprousLinks,
+  NCollection_Map<Standard_Integer>&          theDeadLinks,
+  Standard_Integer&                           theNextPivotNode,
+  gp_Vec2d&                                   theNextLinkDir,
+  Bnd_B2d&                                    theNextLinkBndBox )
 {
-  Standard_Integer e[3];
-  Standard_Boolean o[3];
-  Standard_Integer aGeneralEdgeId = -1;
+  // Find the next link having the greatest angle 
+  // respect to a direction of a reference one
+  Standard_Real aMaxAngle = myIsPositiveOrientation ?
+    RealFirst() : RealLast();
 
-  Standard_Integer anMainEdgeId = Abs( theEdgeInfo );
-  Standard_Boolean anIsForward = theEdgeInfo > 0;
-  const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( anMainEdgeId );
-  for ( Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; ++j )
+  Standard_Integer aNextLinkId = 0;
+  BRepMesh_ListOfInteger::Iterator aLinkIt( myMeshData->LinkNeighboursOf( thePivotNode ) );
+  for ( ; aLinkIt.More(); aLinkIt.Next() )
   {
-    Standard_Integer anElemId = aPair.Index(j);
-    if( anElemId < 0 )
+    const Standard_Integer& aNeighbourLinkInfo = aLinkIt.Value();
+    Standard_Integer        aNeighbourLinkId   = Abs( aNeighbourLinkInfo );
+
+    if ( theDeadLinks.Contains( aNeighbourLinkId ) ||
+       ( !theSkipped.IsNull() && theSkipped->Contains( aNeighbourLinkId ) ) )
+    {
+      continue;
+    }
+      
+    Standard_Boolean isLeprous = theLeprousLinks.Contains( aNeighbourLinkId );
+    if ( isSkipLeprous && isLeprous )
       continue;
 
-    GetTriangle( anElemId ).Edges( e[0], e[1], e[2],
-                                   o[0], o[1], o[2] );
+    const BRepMesh_Edge& aNeighbourLink = GetEdge( aNeighbourLinkId );
 
-    for ( Standard_Integer anIndex = 0; anIndex < 3; ++anIndex )
+    // Determine whether the link belongs to the mesh
+    if ( aNeighbourLink.Movability() == BRepMesh_Free &&
+         myMeshData->ElemConnectedTo( aNeighbourLinkInfo ).IsEmpty() )
     {
-      if ( e[anIndex] == anMainEdgeId && o[anIndex] == anIsForward )
-      {
-        // triangle detected
-        DeleteTriangle( anElemId, theLoopEdges );
-        aGeneralEdgeId = anIndex;
-        break;
-      }
+      theDeadLinks.Add( aNeighbourLinkId );
+      continue;
     }
-       
-    if ( aGeneralEdgeId != -1 )
-      break;
-  }
 
-  if ( aGeneralEdgeId != -1 )
-  {
-    // delete triangles around of aPivotNode starting from the bounding one
-    // define edge connected to vertes
-    Standard_Integer anEdgeId = -1;
-    for ( Standard_Integer i = 0; i < 2; ++i )
+    Standard_Integer anOtherNode = aNeighbourLink.FirstNode();
+    if ( anOtherNode == thePivotNode )
+      anOtherNode = aNeighbourLink.LastNode();
+
+    gp_Vec2d aCurLinkDir( GetVertex( anOtherNode ).Coord() - 
+      thePivotVertex.Coord() );
+
+    if ( aCurLinkDir.SquareMagnitude() < Precision2 )
     {
-      Standard_Integer aTmpEdgeId = e[(aGeneralEdgeId + i + 1) % 3];
-      const BRepMesh_Edge& anEdge = GetEdge( aTmpEdgeId );
-      if ( anEdge.FirstNode() == thePivotNode || 
-           anEdge.LastNode()  == thePivotNode )
+      theDeadLinks.Add( aNeighbourLinkId );
+      continue;
+    }
+
+    if ( !isLeprous )
+      theLeprousLinks.Add( aNeighbourLinkId ); 
+
+    Standard_Real    anAngle    = theRefLinkDir.Angle( aCurLinkDir );
+    Standard_Boolean isFrontier = 
+      ( aNeighbourLink.Movability() == BRepMesh_Frontier );
+
+    Standard_Boolean isCheckPointOnEdge = Standard_True;
+    if ( isFrontier )
+    {
+      if ( Abs( Abs(anAngle) - M_PI ) < Precision::Angular() )
       {
-        anEdgeId = aTmpEdgeId;
+        // Glued constrains - don't check intersection
+        isCheckPointOnEdge = Standard_False;
+        anAngle            = Abs( anAngle );
       }
+    }
 
-      if ( theInfectedEdges.Contains( aTmpEdgeId ) )
-        theInfectedEdges.Remove( aTmpEdgeId );
-      else
-        theInfectedEdges.Add( aTmpEdgeId );
+    if ( ( myIsPositiveOrientation && anAngle <= aMaxAngle ) ||
+         (!myIsPositiveOrientation && anAngle >= aMaxAngle ) )
+    {
+      continue;
     }
 
-    if ( isFirstPass )
-      return;
+    Standard_Boolean isCheckEndPoints = ( anOtherNode != theFirstNode );
 
-    while ( anEdgeId > 0 )
+    Bnd_B2d aBox;
+    Standard_Boolean isNotIntersect =
+      checkIntersection( aNeighbourLink, thePolygon, theBoxes, 
+      isCheckEndPoints, isCheckPointOnEdge, Standard_True, aBox );
+    
+    if( isNotIntersect )
     {
-      const BRepMesh_PairOfIndex& aFreePair = myMeshData->ElemConnectedTo( anEdgeId );
-      Standard_Integer anOldEdge = anEdgeId;
-      anEdgeId = -1;
+      aMaxAngle         = anAngle;
 
-      for ( Standard_Integer aPairIndex = 1; aPairIndex <= aFreePair.Extent(); ++aPairIndex )
-      {
-        Standard_Integer anElemId = aFreePair.Index( aPairIndex );
-        if( anElemId < 0 )
-          continue;
-          
-        Standard_Integer e1[3];
-        Standard_Boolean o1[3];
-        GetTriangle( anElemId ).Edges( e1[0], e1[1], e1[2],
-                                       o1[0], o1[1], o1[2] );
-        
-        DeleteTriangle( anElemId, theLoopEdges );
+      theNextLinkDir    = aCurLinkDir;
+      theNextPivotNode  = anOtherNode;
+      theNextLinkBndBox = aBox;
 
-        for ( Standard_Integer i = 0; i < 3; ++i )
-        {
-          if ( e1[i] == anOldEdge )
-          {
-            for ( Standard_Integer j = 0; j < 2; ++j )
-            {              
-              Standard_Integer aTmpEdgeId = e1[(j + i + 1) % 3];
-              const BRepMesh_Edge& anEdge = GetEdge( aTmpEdgeId );
-              if ( anEdge.FirstNode() == thePivotNode || 
-                   anEdge.LastNode()  == thePivotNode )
-              {
-                anEdgeId = aTmpEdgeId;
-              }
-            
-              if ( theInfectedEdges.Contains( aTmpEdgeId ) )
-                theInfectedEdges.Remove( aTmpEdgeId );
-              else
-                theInfectedEdges.Add( aTmpEdgeId );
-            }
-            break;
-          }
-        }
-      }
+      aNextLinkId       = ( aNeighbourLink.FirstNode() == thePivotNode ) ?
+        aNeighbourLinkId : -aNeighbourLinkId;
     }
   }
+
+  return aNextLinkId;
 }
 
 //=======================================================================
-//function : CleanupPolygon
-//purpose  : Remove internal triangles from the given polygon
+//function : checkIntersection
+//purpose  : Check is the given link intersects the polygon boundaries.
+//           Returns bounding box for the given link trough the 
+//           <theLinkBndBox> parameter.
 //=======================================================================
-void BRepMesh_Delaun::CleanupPolygon( const TColStd_SequenceOfInteger& thePolygon,
-                                      TColStd_MapOfInteger& theInfectedEdges,
-                                      BRepMesh_MapOfIntegerInteger& theLoopEdges )
+Standard_Boolean BRepMesh_Delaun::checkIntersection( 
+  const BRepMesh_Edge&                 theLink,
+  const TColStd_SequenceOfInteger&     thePolygon,
+  const NCollection_Sequence<Bnd_B2d>& thePolyBoxes,
+  const Standard_Boolean               isConsiderEndPointTouch,
+  const Standard_Boolean               isConsiderPointOnEdge,
+  const Standard_Boolean               isSkipLastEdge,
+  Bnd_B2d&                             theLinkBndBox ) const
 {
-  Standard_Integer aPolyLen = thePolygon.Length();
-
-  TColStd_MapOfInteger anIgnoredEdges;
-  NCollection_Map<Standard_Integer> aPolyVertices;
-  for ( Standard_Integer i = 1; i <= aPolyLen; ++i )
-  {
-    Standard_Integer aPolyEdgeId = Abs( thePolygon(i) );
-    anIgnoredEdges.Add( aPolyEdgeId );
+  theLinkBndBox.Add( GetVertex( theLink.FirstNode() ).Coord() );
+  theLinkBndBox.Add( GetVertex( theLink.LastNode()  ).Coord() );
 
-    const BRepMesh_Edge& aPolyEdge = GetEdge( aPolyEdgeId );
-    aPolyVertices.Add( aPolyEdge.FirstNode() );
-    aPolyVertices.Add( aPolyEdge.LastNode() );
+  Standard_Integer aPolyLen = thePolygon.Length();
+  // Don't check intersection with the last link
+  if ( isSkipLastEdge )
+    --aPolyLen;
 
-    if ( theInfectedEdges.Contains( aPolyEdgeId ) )
-      theInfectedEdges.Remove( aPolyEdgeId );
-  }
+  Standard_Boolean isFrontier = 
+    ( theLink.Movability() == BRepMesh_Frontier );
 
-  Standard_Real aRefTotalAngle = 2 * M_PI;
-  TColStd_MapIteratorOfMapOfInteger anInfectedIt( theInfectedEdges );
-  for ( ; anInfectedIt.More(); anInfectedIt.Next() )
+  for ( Standard_Integer aPolyIt = 1; aPolyIt <= aPolyLen; ++aPolyIt )
   {
-    Standard_Integer anEdgeId = anInfectedIt.Key();
-    const BRepMesh_Edge& anInfectedEdge = GetEdge( anEdgeId );
-    Standard_Integer anEdgeVertices[2] = { anInfectedEdge.FirstNode(),
-                                           anInfectedEdge.LastNode() };
-
-    Standard_Boolean isToExclude = Standard_False;
-    for ( Standard_Integer i = 0; i < 2; ++i )
+    if ( !theLinkBndBox.IsOut( thePolyBoxes.Value( aPolyIt ) ) )
     {
-      Standard_Integer aCurVertex = anEdgeVertices[i];
-      if ( aPolyVertices.Contains( aCurVertex ) )
-        continue;
+      // intersection is possible...
+      Standard_Integer aPolyLinkId   = Abs( thePolygon( aPolyIt ) );
+      const BRepMesh_Edge& aPolyLink = GetEdge( aPolyLinkId );
 
-      gp_XY aCenterPointXY = GetVertex( aCurVertex ).Coord();
-      Standard_Real aTotalAng = 0.0;
-
-      for ( Standard_Integer i = 1; i <= aPolyLen; ++i )
-      {
-        Standard_Integer aPolyEdgeId = thePolygon(i);
-        const BRepMesh_Edge& aPolyEdge = GetEdge( Abs( aPolyEdgeId ) );
-        Standard_Integer aStartPoint, anEndPoint;
-        if ( aPolyEdgeId >= 0 )
-        {
-          aStartPoint = aPolyEdge.FirstNode();
-          anEndPoint  = aPolyEdge.LastNode();
-        }
-        else
-        {
-          aStartPoint = aPolyEdge.LastNode();
-          anEndPoint  = aPolyEdge.FirstNode();
-        }
+      // skip intersections between frontier edges
+      if ( aPolyLink.Movability() == BRepMesh_Frontier && isFrontier )
+        continue;
 
-        gp_XY aStartV = GetVertex( aStartPoint ).Coord() - aCenterPointXY;
-        gp_XY anEndV  = GetVertex( anEndPoint ).Coord()  - aCenterPointXY;
+      gp_Pnt2d anIntPnt;
+      IntFlag aIntFlag = intSegSeg( theLink, aPolyLink, 
+        isConsiderEndPointTouch, isConsiderPointOnEdge, anIntPnt );
 
-        Standard_Real anAngle = gp_Vec2d( anEndV ).Angle( gp_Vec2d( aStartV ) );
-        aTotalAng += anAngle;
-      }
-      
-      if ( Abs( aRefTotalAngle - aTotalAng ) > Precision::Angular() )
-        isToExclude = Standard_True;
+      if ( aIntFlag != BRepMesh_Delaun::NoIntersection )
+        return Standard_False;
     }
-
-    if ( isToExclude )
-      anIgnoredEdges.Add( anEdgeId );
-  }
-
-
-  anInfectedIt.Initialize( theInfectedEdges );
-  for ( ; anInfectedIt.More(); anInfectedIt.Next() )
-  {
-    Standard_Integer anEdgeId = anInfectedIt.Key();
-    KillInternalTriangles( anEdgeId, anIgnoredEdges, theLoopEdges);
   }
 
-  BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( theLoopEdges );
-  for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
-  {
-    if ( myMeshData->ElemConnectedTo( aLoopEdgesIt.Key() ).IsEmpty() )
-      myMeshData->RemoveLink( aLoopEdgesIt.Key() );
-  }
+  // Ok, no intersection
+  return Standard_True;
 }
 
 //=======================================================================
-//function : MeshLeftPolygonOf
-//purpose  : Triangulation of the aPolygonon to the left of <theEdgeIndex>.(material side)
+//function : addTriangle
+//purpose  : Add a triangle based on the given oriented edges into mesh
 //=======================================================================
-void BRepMesh_Delaun::MeshLeftPolygonOf( const Standard_Integer theEdgeIndex,
-                                         const Standard_Boolean isForward )
+inline void BRepMesh_Delaun::addTriangle( const Standard_Integer (&theEdgesId)[3],
+                                          const Standard_Boolean (&theEdgesOri)[3],
+                                          const Standard_Integer (&theNodesId)[3] )
 {
-  const BRepMesh_Edge& anEdge = GetEdge( theEdgeIndex );
-  NCollection_Map<Standard_Integer> aDealLinks;
-  TColStd_SequenceOfInteger aPolygon;
-  BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
-  TColStd_MapOfInteger anUsedEdges;
-  TColStd_MapOfInteger anInfectedEdges;
-  
-  // Find the aPolygonon
-  anUsedEdges.Add( theEdgeIndex );
-  Standard_Integer aFirstNode, aStartNode, aPivotNode;
-  
-  if ( isForward )
+  Standard_Integer aNewTriangleId = 
+    myMeshData->AddElement( BRepMesh_Triangle( 
+      theEdgesId[0],  theEdgesId[1],  theEdgesId[2], 
+      theEdgesOri[0], theEdgesOri[1], theEdgesOri[2],
+      BRepMesh_Free ) );
+
+  Standard_Boolean isAdded = myCircles.Add( 
+    GetVertex( theNodesId[0] ).Coord(), 
+    GetVertex( theNodesId[1] ).Coord(),
+    GetVertex( theNodesId[2] ).Coord(),
+    aNewTriangleId );
+    
+  if ( !isAdded )
+    myMeshData->RemoveElement( aNewTriangleId );
+}
+
+//=======================================================================
+//function : cleanupPolygon
+//purpose  : Remove internal triangles from the given polygon
+//=======================================================================
+void BRepMesh_Delaun::cleanupPolygon( const TColStd_SequenceOfInteger&     thePolygon,
+                                      const NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
+{
+  Standard_Integer aPolyLen = thePolygon.Length();
+  if ( aPolyLen < 3 )
+    return;
+
+  BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
+  NCollection_Map<Standard_Integer>    anIgnoredEdges;
+  NCollection_Map<Standard_Integer>    aPolyVerticesFindMap;
+  NCollection_Vector<Standard_Integer> aPolyVertices;
+  // Collect boundary vertices of the polygon
+  for ( Standard_Integer aPolyIt = 1; aPolyIt <= aPolyLen; ++aPolyIt )
   {
-    aPolygon.Append( theEdgeIndex );
-    aStartNode = anEdge.FirstNode();
-    aPivotNode = anEdge.LastNode();
+    Standard_Integer aPolyEdgeInfo = thePolygon( aPolyIt );
+    Standard_Integer aPolyEdgeId   = Abs( aPolyEdgeInfo );
+    anIgnoredEdges.Add( aPolyEdgeId );
+
+    Standard_Boolean isForward = ( aPolyEdgeInfo > 0 );
+    const BRepMesh_PairOfIndex& aPair = 
+      myMeshData->ElemConnectedTo( aPolyEdgeId );
+
+    Standard_Integer anElemIt = 1;
+    for ( ; anElemIt <= aPair.Extent(); ++anElemIt )
+    {
+      Standard_Integer anElemId = aPair.Index( anElemIt );
+      if ( anElemId < 0 )
+        continue;
+
+      Standard_Integer anEdges[3];
+      Standard_Boolean anEdgesOri[3];
+      GetTriangle( anElemId ).Edges( anEdges[0], anEdges[1], anEdges[2],
+        anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
+
+      Standard_Integer isTriangleFound = Standard_False;
+      for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
+      {
+        if ( anEdges[anEdgeIt]    == aPolyEdgeId && 
+             anEdgesOri[anEdgeIt] == isForward )
+        {
+          isTriangleFound = Standard_True;
+          deleteTriangle( anElemId, aLoopEdges );
+          break;
+        }
+      }
+
+      if ( isTriangleFound )
+        break;
+    }
+
+    // Skip a neighbor link to extract unique vertices each time
+    if ( aPolyIt % 2 )
+    {
+      const BRepMesh_Edge& aPolyEdge = GetEdge( aPolyEdgeId );
+      Standard_Integer aFirstVertex  = aPolyEdge.FirstNode();
+      Standard_Integer aLastVertex   = aPolyEdge.LastNode();
+
+      aPolyVerticesFindMap.Add( aFirstVertex );
+      aPolyVerticesFindMap.Add( aLastVertex );
+
+      if ( aPolyEdgeInfo > 0 )
+      {
+        aPolyVertices.Append( aFirstVertex );
+        aPolyVertices.Append( aLastVertex );
+      }
+      else
+      {
+        aPolyVertices.Append( aLastVertex );
+        aPolyVertices.Append( aFirstVertex );
+      }
+    }
   }
-  else
+
+  // Make closed sequence
+  if ( aPolyVertices.First() != aPolyVertices.Last() )
+    aPolyVertices.Append( aPolyVertices.First() );
+
+  NCollection_Map<Standard_Integer> aSurvivedLinks( anIgnoredEdges );
+
+  Standard_Integer aPolyVertIt          = 0;
+  Standard_Integer anUniqueVerticesNum  = aPolyVertices.Length() - 1;
+  for ( ; aPolyVertIt < anUniqueVerticesNum; ++aPolyVertIt )
+  {
+    killTrianglesAroundVertex( aPolyVertices( aPolyVertIt ),
+      aPolyVertices, aPolyVerticesFindMap, thePolygon,
+      thePolyBoxes, aSurvivedLinks, aLoopEdges );
+  }
+
+  BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
+  for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
   {
-    aPolygon.Append( -theEdgeIndex );
-    aStartNode = anEdge.LastNode();
-    aPivotNode = anEdge.FirstNode();
+    const Standard_Integer& aLoopEdgeId = aLoopEdgesIt.Key();
+    if ( anIgnoredEdges.Contains( aLoopEdgeId ) )
+      continue;
+
+    if ( myMeshData->ElemConnectedTo( aLoopEdgeId ).IsEmpty() )
+      myMeshData->RemoveLink( aLoopEdgesIt.Key() );
   }
-  aFirstNode = aStartNode;
+}
 
-  const BRepMesh_Vertex& aStartVertex = GetVertex( aFirstNode );
-  const BRepMesh_Vertex& anEndVertex  = GetVertex( aPivotNode );
+//=======================================================================
+//function : killTrianglesAroundVertex
+//purpose  : Remove all triangles and edges that are placed
+//           inside the polygon or crossed it.
+//=======================================================================
+void BRepMesh_Delaun::killTrianglesAroundVertex( 
+  const Standard_Integer                        theZombieNodeId,
+  const NCollection_Vector<Standard_Integer>&   thePolyVertices,
+  const NCollection_Map<Standard_Integer>&      thePolyVerticesFindMap,
+  const TColStd_SequenceOfInteger&              thePolygon,
+  const NCollection_Sequence<Bnd_B2d>&          thePolyBoxes,
+  NCollection_Map<Standard_Integer>&            theSurvivedLinks,
+  BRepMesh_MapOfIntegerInteger&                 theLoopEdges )
+{
+  BRepMesh_ListOfInteger::Iterator aNeighborsIt = 
+      myMeshData->LinkNeighboursOf( theZombieNodeId );
 
-  Standard_Integer aRefOtherNode = 0, anOtherNode = 0;
-  // Check the presence of the previous edge <theEdgeIndex> :
-  BRepMesh_ListOfInteger::Iterator aLinkIt( myMeshData->LinkNeighboursOf( aStartNode ) );
-  for ( ; aLinkIt.More(); aLinkIt.Next() )
+  // Try to infect neighbor nodes
+  NCollection_Vector<Standard_Integer> aVictimNodes;
+  for ( ; aNeighborsIt.More(); aNeighborsIt.Next() )
   {
-    if ( aLinkIt.Value() != theEdgeIndex )
+    const Standard_Integer& aNeighborLinkId = aNeighborsIt.Value();
+    if ( theSurvivedLinks.Contains( aNeighborLinkId ) )
+      continue;
+
+    const BRepMesh_Edge& aNeighborLink = GetEdge( aNeighborLinkId );
+    if ( aNeighborLink.Movability() == BRepMesh_Frontier )
     {
-      const BRepMesh_Edge& aNextEdge = GetEdge( aLinkIt.Value() );
-      anOtherNode = aNextEdge.LastNode();
-      if ( anOtherNode == aStartNode )
-        anOtherNode = aNextEdge.FirstNode();
-      break;
+      // Though, if it lies onto the polygon boundary -
+      // take its triangles
+      Bnd_B2d aBox;
+      Standard_Boolean isNotIntersect =
+        checkIntersection( aNeighborLink, thePolygon, 
+        thePolyBoxes, Standard_False, Standard_True, 
+        Standard_False, aBox );
+
+      if ( isNotIntersect )
+      {
+        // Don't touch it
+        theSurvivedLinks.Add( aNeighborLinkId );
+        continue;
+      }
     }
+    else
+    {
+      Standard_Integer anOtherNode = aNeighborLink.FirstNode();
+      if ( anOtherNode == theZombieNodeId )
+        anOtherNode = aNeighborLink.LastNode();
+
+      // Possible sacrifice
+      if ( !thePolyVerticesFindMap.Contains( anOtherNode ) )
+      {
+        if ( isVertexInsidePolygon( anOtherNode, thePolyVertices ) )
+        {
+          // Got you!
+          aVictimNodes.Append( anOtherNode );
+        }
+        else
+        {
+          // Lucky. But it may intersect the polygon boundary.
+          // Let's check it!
+          killTrianglesOnIntersectingLinks( aNeighborLinkId, 
+            aNeighborLink, anOtherNode, thePolygon, 
+            thePolyBoxes, theSurvivedLinks, theLoopEdges );
+
+          continue;
+        }
+      }
+    }
+
+    // Add link to the survivers to avoid cycling
+    theSurvivedLinks.Add( aNeighborLinkId );
+    killLinkTriangles( aNeighborLinkId, theLoopEdges );
   }
-  if ( anOtherNode == 0 )
-    return;
 
+  // Go and do your job!
+  NCollection_Vector<Standard_Integer>::Iterator aVictimIt( aVictimNodes );
+  for ( ; aVictimIt.More(); aVictimIt.Next() )
+  {
+    killTrianglesAroundVertex( aVictimIt.Value(), thePolyVertices,
+      thePolyVerticesFindMap, thePolygon, thePolyBoxes,
+      theSurvivedLinks, theLoopEdges );
+  }
+}
+
+//=======================================================================
+//function : isVertexInsidePolygon
+//purpose  : Checks is the given vertex lies inside the polygon
+//=======================================================================
+Standard_Boolean BRepMesh_Delaun::isVertexInsidePolygon( 
+  const Standard_Integer&                     theVertexId,
+  const NCollection_Vector<Standard_Integer>& thePolygonVertices ) const
+{
+  Standard_Integer aPolyLen = thePolygonVertices.Length();
+  if ( aPolyLen < 3 )
+    return Standard_False;
+
+
+  const gp_XY aCenterPointXY = GetVertex( theVertexId ).Coord();
+
+  const BRepMesh_Vertex& aFirstVertex = GetVertex( thePolygonVertices( 0 ) );
+  gp_Vec2d aPrevVertexDir( aFirstVertex.Coord() - aCenterPointXY );
+  if ( aPrevVertexDir.SquareMagnitude() < Precision2 )
+    return Standard_True;
+
+  Standard_Real aTotalAng = 0.0;
+  for ( Standard_Integer aPolyIt = 1; aPolyIt < aPolyLen; ++aPolyIt )
+  {
+    const BRepMesh_Vertex& aPolyVertex = GetVertex( thePolygonVertices( aPolyIt ) );
+
+    gp_Vec2d aCurVertexDir( aPolyVertex.Coord() - aCenterPointXY );
+    if ( aCurVertexDir.SquareMagnitude() < Precision2 )
+      return Standard_True;
 
-  gp_XY aVEdge( anEndVertex.Coord() );
-  aVEdge.Subtract( aStartVertex.Coord() );
-  gp_XY aPrevVEdge( aVEdge );
-  gp_XY aRefCurrVEdge, aCurrVEdge;
+    aTotalAng     += aCurVertexDir.Angle( aPrevVertexDir );
+    aPrevVertexDir = aCurVertexDir;
+  }
   
-  Standard_Integer aCurrEdgeId = theEdgeIndex;
-  Standard_Integer aNextEdgeId;
+  if ( Abs( Angle2PI - aTotalAng ) > Precision::Angular() )
+    return Standard_False;
 
-  // Find the nearest to <theEdgeIndex> closed aPolygonon :
-  Standard_Boolean isInMesh, isNotInters;
-  Standard_Real anAngle, aRefAngle;
+  return Standard_True;
+}
 
-  NCollection_Vector <Bnd_Box2d> aBoxes;
-  fillBndBox( aBoxes, aStartVertex, anEndVertex );
-    
-  while ( aPivotNode != aFirstNode )
+//=======================================================================
+//function : killTrianglesOnIntersectingLinks
+//purpose  : Checks is the given link crosses the polygon boundary.
+//           If yes, kills its triangles and checks neighbor links on
+//           boundary intersection. Does nothing elsewhere.
+//=======================================================================
+void BRepMesh_Delaun::killTrianglesOnIntersectingLinks( 
+  const Standard_Integer&               theLinkToCheckId, 
+  const BRepMesh_Edge&                  theLinkToCheck,
+  const Standard_Integer&               theEndPoint,
+  const TColStd_SequenceOfInteger&      thePolygon,
+  const NCollection_Sequence<Bnd_B2d>&  thePolyBoxes,
+  NCollection_Map<Standard_Integer>&    theSurvivedLinks,
+  BRepMesh_MapOfIntegerInteger&         theLoopEdges )
+{
+  if ( theSurvivedLinks.Contains( theLinkToCheckId ) )
+    return;
+
+  Bnd_B2d aBox;
+  Standard_Boolean isNotIntersect =
+    checkIntersection( theLinkToCheck, thePolygon, 
+      thePolyBoxes, Standard_False, Standard_False, 
+      Standard_False, aBox );
+
+  theSurvivedLinks.Add( theLinkToCheckId );
+
+  if ( isNotIntersect )
+    return;
+
+  killLinkTriangles( theLinkToCheckId, theLoopEdges );
+
+  BRepMesh_ListOfInteger::Iterator aNeighborsIt = 
+    myMeshData->LinkNeighboursOf( theEndPoint );
+
+  for ( ; aNeighborsIt.More(); aNeighborsIt.Next() )
   {
-    aNextEdgeId = 0;
-    if ( myPositiveOrientation )
-      aRefAngle = RealFirst();
-    else
-      aRefAngle = RealLast();
+    const Standard_Integer& aNeighborLinkId = aNeighborsIt.Value();
+    const BRepMesh_Edge&    aNeighborLink   = GetEdge( aNeighborLinkId );
+    Standard_Integer anOtherNode = aNeighborLink.FirstNode();
+    if ( anOtherNode == theEndPoint )
+      anOtherNode = aNeighborLink.LastNode();
+
+    killTrianglesOnIntersectingLinks( aNeighborLinkId, 
+      aNeighborLink, anOtherNode, thePolygon, 
+      thePolyBoxes, theSurvivedLinks, theLoopEdges );
+  }
+}
+
+//=======================================================================
+//function : killLinkTriangles
+//purpose  : Kill triangles bound to the given link.
+//=======================================================================
+void BRepMesh_Delaun::killLinkTriangles( 
+  const Standard_Integer&       theLinkId, 
+  BRepMesh_MapOfIntegerInteger& theLoopEdges )
+{
+  const BRepMesh_PairOfIndex& aPair = 
+      myMeshData->ElemConnectedTo( theLinkId );
 
-    const BRepMesh_Vertex& aPivotVertex = GetVertex( aPivotNode );
+  Standard_Integer anElemNb = aPair.Extent();
+  for ( Standard_Integer aPairIt = 1; aPairIt <= anElemNb; ++aPairIt )
+  {
+    Standard_Integer anElemId = aPair.FirstIndex();
+    if ( anElemId < 0 )
+      continue;
 
-    // Find the next edge with the greatest angle with <theEdgeIndex>
-    // and non intersecting <theEdgeIndex> :
-    
-    aLinkIt.Init( myMeshData->LinkNeighboursOf( aPivotNode ) );
-    for ( ; aLinkIt.More(); aLinkIt.Next() )
+    deleteTriangle( anElemId, theLoopEdges );
+  }
+}
+
+//=======================================================================
+//function : getOrientedNodes
+//purpose  : Returns start and end nodes of the given edge in respect to 
+//           its orientation.
+//=======================================================================
+void BRepMesh_Delaun::getOrientedNodes(const BRepMesh_Edge&   theEdge,
+                                       const Standard_Boolean isForward,
+                                       Standard_Integer       *theNodes) const
+{
+  if ( isForward )
+  {
+    theNodes[0] = theEdge.FirstNode();
+    theNodes[1] = theEdge.LastNode();
+  }
+  else
+  {
+    theNodes[0] = theEdge.LastNode();
+    theNodes[1] = theEdge.FirstNode();
+  }
+}
+
+//=======================================================================
+//function : processLoop
+//purpose  : Processes loop within the given polygon formed by range of 
+//           its links specified by start and end link indices.
+//=======================================================================
+void BRepMesh_Delaun::processLoop(const Standard_Integer               theLinkFrom,
+                                  const Standard_Integer               theLinkTo,
+                                  const TColStd_SequenceOfInteger&     thePolygon,
+                                  const NCollection_Sequence<Bnd_B2d>& thePolyBoxes)
+{
+  Standard_Integer aNbOfLinksInLoop = theLinkTo - theLinkFrom - 1;
+  if ( aNbOfLinksInLoop < 3 )
+    return;
+
+  TColStd_SequenceOfInteger     aPolygon;
+  NCollection_Sequence<Bnd_B2d> aPolyBoxes;
+  for ( ; aNbOfLinksInLoop > 0; --aNbOfLinksInLoop )
+  {
+    Standard_Integer aLoopLinkIndex = theLinkFrom + aNbOfLinksInLoop;
+    aPolygon  .Prepend( thePolygon  ( aLoopLinkIndex ) );
+    aPolyBoxes.Prepend( thePolyBoxes( aLoopLinkIndex ) );
+  }
+  meshPolygon( aPolygon, aPolyBoxes );
+}
+
+//=======================================================================
+//function : createAndReplacePolygonLink
+//purpose  : Creates new link based on the given nodes and updates the 
+//           given polygon.
+//=======================================================================
+Standard_Integer BRepMesh_Delaun::createAndReplacePolygonLink(
+  const Standard_Integer         *theNodes,
+  const gp_Pnt2d                 *thePnts,
+  const Standard_Integer         theRootIndex,
+  const ReplaceFlag              theReplaceFlag,
+  TColStd_SequenceOfInteger&     thePolygon,
+  NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
+{
+  Standard_Integer aNewEdgeId = 
+    myMeshData->AddLink( BRepMesh_Edge(
+    theNodes[0], theNodes[1], BRepMesh_Free ) );
+
+  Bnd_B2d aNewBox;
+  aNewBox.Add( thePnts[0] );
+  aNewBox.Add( thePnts[1] );
+
+  switch ( theReplaceFlag )
+  {
+  case BRepMesh_Delaun::Replace:
+    thePolygon  .SetValue( theRootIndex, aNewEdgeId );
+    thePolyBoxes.SetValue( theRootIndex, aNewBox );
+    break;
+
+  case BRepMesh_Delaun::InsertAfter:
+    thePolygon  .InsertAfter( theRootIndex, aNewEdgeId );
+    thePolyBoxes.InsertAfter( theRootIndex, aNewBox );
+    break;
+
+  case BRepMesh_Delaun::InsertBefore:
+    thePolygon  .InsertBefore( theRootIndex, aNewEdgeId );
+    thePolyBoxes.InsertBefore( theRootIndex, aNewBox );
+    break;
+  }
+
+  return aNewEdgeId;
+}
+
+//=======================================================================
+//function : meshPolygon
+//purpose  : 
+//=======================================================================
+void BRepMesh_Delaun::meshPolygon( TColStd_SequenceOfInteger&            thePolygon,
+                                   NCollection_Sequence<Bnd_B2d>&        thePolyBoxes,
+                                   BRepMesh_Delaun::HandleOfMapOfInteger theSkipped )
+{
+  // Check is the source polygon elementary
+  if ( meshElementaryPolygon( thePolygon ) )
+    return;
+
+  // Check and correct boundary edges
+  Standard_Integer aPolyLen = thePolygon.Length();
+  const Standard_Real aPolyArea      = Abs( polyArea( thePolygon, 1, aPolyLen ) );
+  const Standard_Real aSmallLoopArea = 0.001 * aPolyArea;
+  for ( Standard_Integer aPolyIt = 1; aPolyIt < aPolyLen; ++aPolyIt )
+  {
+    Standard_Integer aCurEdgeInfo = thePolygon( aPolyIt );
+    Standard_Integer aCurEdgeId   = Abs( aCurEdgeInfo );
+    const BRepMesh_Edge* aCurEdge = &GetEdge( aCurEdgeId );
+    if ( aCurEdge->Movability() != BRepMesh_Frontier )
+      continue;
+
+    Standard_Integer aCurNodes[2];
+    getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aCurNodes );
+
+    gp_Pnt2d aCurPnts[2] = { 
+      GetVertex(aCurNodes[0]).Coord(),
+      GetVertex(aCurNodes[1]).Coord()
+    };
+
+    gp_Vec2d aCurVec( aCurPnts[0], aCurPnts[1] );
+
+    // check further links
+    Standard_Integer aNextPolyIt = aPolyIt + 1;
+    for ( ; aNextPolyIt <= aPolyLen; ++aNextPolyIt )
     {
-      Standard_Integer aNextLink = aLinkIt.Value();
-      Standard_Integer aNextLinkId = Abs( aNextLink );
-      if ( aDealLinks.Contains( aNextLinkId ) )
+      Standard_Integer aNextEdgeInfo = thePolygon( aNextPolyIt );
+      Standard_Integer aNextEdgeId   = Abs( aNextEdgeInfo );
+      const BRepMesh_Edge* aNextEdge = &GetEdge( aNextEdgeId );
+      if ( aNextEdge->Movability() != BRepMesh_Frontier )
+        continue;
+
+      Standard_Integer aNextNodes[2];
+      getOrientedNodes( *aNextEdge, aNextEdgeInfo > 0, aNextNodes );
+
+      gp_Pnt2d aNextPnts[2] = { 
+        GetVertex(aNextNodes[0]).Coord(),
+        GetVertex(aNextNodes[1]).Coord() 
+      };
+
+      gp_Pnt2d anIntPnt;
+      IntFlag aIntFlag = intSegSeg( *aCurEdge, *aNextEdge, 
+        Standard_False, Standard_True, anIntPnt );
+
+      if ( aIntFlag == BRepMesh_Delaun::NoIntersection )
         continue;
 
-      if ( aNextLinkId != aCurrEdgeId )
+      Standard_Boolean isRemoveFromFirst  = Standard_False;
+      Standard_Boolean isAddReplacingEdge = Standard_True;
+      Standard_Integer aIndexToRemoveTo   = aNextPolyIt;
+      if ( aIntFlag == BRepMesh_Delaun::Cross )
       {
-        isNotInters = Standard_True;
-        const BRepMesh_Edge& aNextEdge = GetEdge( aNextLinkId );
+        Standard_Real aLoopArea = polyArea( thePolygon, aPolyIt + 1, aNextPolyIt );
+        gp_Vec2d aVec1( anIntPnt, aCurPnts [1] );
+        gp_Vec2d aVec2( anIntPnt, aNextPnts[0] );
 
-        isInMesh = Standard_True;
-        
-        //define if link exist in mesh      
-        if ( aNextEdge.Movability() == BRepMesh_Free )
+        aLoopArea += ( aVec1 ^ aVec2 ) / 2.;
+        if ( Abs( aLoopArea ) > aSmallLoopArea )
         {
-          if ( myMeshData->ElemConnectedTo( aLinkIt.Value() ).IsEmpty() ) 
-            isInMesh = Standard_False;
+          aNextNodes[1] = aCurNodes[0];
+          aNextPnts [1] = aCurPnts [0];
+
+          aNextEdgeId = Abs( createAndReplacePolygonLink( aNextNodes, aNextPnts,
+            aNextPolyIt, BRepMesh_Delaun::Replace, thePolygon, thePolyBoxes ) );
+
+          processLoop( aPolyIt, aNextPolyIt, thePolygon, thePolyBoxes );
+          return;
         }
 
-        if ( isInMesh )
-        {
-          anOtherNode = aNextEdge.FirstNode();
-          if ( anOtherNode == aPivotNode )
-            anOtherNode = aNextEdge.LastNode();
+        Standard_Real aDist1 = anIntPnt.SquareDistance(aNextPnts[0]);
+        Standard_Real aDist2 = anIntPnt.SquareDistance(aNextPnts[1]);
 
-          aCurrVEdge = GetVertex( anOtherNode ).Coord();
-          aCurrVEdge.Subtract( aPivotVertex.Coord() );
-          
-          if ( aCurrVEdge.Modulus() >= gp::Resolution() &&
-               aPrevVEdge.Modulus() >= gp::Resolution() )
+        // Choose node with lower distance
+        const Standard_Boolean isCloseToStart = ( aDist1 < aDist2 );
+        const Standard_Integer aEndPointIndex = isCloseToStart ? 0 : 1;
+        aCurNodes[1] = aNextNodes[aEndPointIndex];
+        aCurPnts [1] = aNextPnts [aEndPointIndex];
+
+        if ( isCloseToStart )
+          --aIndexToRemoveTo;
+
+        // In this context only intersections between frontier edges
+        // are possible. If intersection between edges of different
+        // types occured - treat this case as invalid (i.e. result 
+        // might not reflect the expectations).
+        if ( !theSkipped.IsNull() )
+        {
+          Standard_Integer aSkippedLinkIt = aPolyIt;
+          for ( ; aSkippedLinkIt <= aIndexToRemoveTo; ++aSkippedLinkIt )
+            theSkipped->Add( Abs( thePolygon( aSkippedLinkIt ) ) );
+        }
+      }
+      else if ( aIntFlag == BRepMesh_Delaun::PointOnEdge )
+      {
+        // Indentify chopping link 
+        Standard_Boolean isFirstChopping = Standard_False;
+        Standard_Integer aCheckPointIt = 0;
+        for ( ; aCheckPointIt < 2; ++aCheckPointIt )
+        {
+          gp_Pnt2d& aRefPoint = aCurPnts[aCheckPointIt];
+          // Check is second link touches the first one
+          gp_Vec2d aVec1( aRefPoint, aNextPnts[0] );
+          gp_Vec2d aVec2( aRefPoint, aNextPnts[1] );
+          if ( Abs( aVec1 ^ aVec2 ) < Precision )
+          {
+            isFirstChopping = Standard_True;
+            break;
+          }
+        }
+        if ( isFirstChopping )
+        {
+          // Split second link
+          isAddReplacingEdge = Standard_False;
+          isRemoveFromFirst  = ( aCheckPointIt == 0 );
+
+          Standard_Integer aSplitLink[3] = {
+            aNextNodes[0],
+            aCurNodes [aCheckPointIt],
+            aNextNodes[1]
+          };
+
+          gp_Pnt2d aSplitPnts[3] = {
+            aNextPnts[0],
+            aCurPnts [aCheckPointIt],
+            aNextPnts[1]
+          };
+
+          Standard_Integer aSplitLinkIt = 0;
+          for ( ; aSplitLinkIt < 2; ++aSplitLinkIt )
           {
-            anAngle = gp_Vec2d( aPrevVEdge ).Angle( gp_Vec2d( aCurrVEdge ) );
-
-            if ( ( myPositiveOrientation && anAngle > aRefAngle ) ||
-                ( !myPositiveOrientation && anAngle < aRefAngle ) )
-            {
-              Bnd_Box2d aBox;
-              aBox.Add( gp_Pnt2d( GetVertex( aNextEdge.FirstNode() ).Coord() ) );
-              aBox.Add( gp_Pnt2d( GetVertex( aNextEdge.LastNode()  ).Coord() ) );
-              
-              for ( Standard_Integer k = 0, aLen = aPolygon.Length(); k < aLen; ++k )
-              {
-                if ( !aBox.IsOut( aBoxes.Value( k ) ) )
-                {
-                  // intersection is possible...
-                  if ( IntSegSeg( aNextEdge, GetEdge( Abs( aPolygon( k + 1 ) ) ) ) )
-                  {
-                    isNotInters = Standard_False;
-                    break;
-                  }
-                }
-              }
-              
-              if( isNotInters )
-              {              
-                aRefAngle = anAngle;
-                aRefCurrVEdge = aCurrVEdge;
-                
-                if ( aNextEdge.FirstNode() == aPivotNode )
-                  aNextEdgeId =  aLinkIt.Value();
-                else
-                  aNextEdgeId = -aLinkIt.Value();
-                  
-                aRefOtherNode = anOtherNode;
-              }
-            }
+            createAndReplacePolygonLink( &aSplitLink[aSplitLinkIt],
+              &aSplitPnts[aSplitLinkIt], aNextPolyIt, ( aSplitLinkIt == 0 ) ? 
+              BRepMesh_Delaun::Replace : BRepMesh_Delaun::InsertAfter,
+              thePolygon, thePolyBoxes );
           }
+
+          processLoop( aPolyIt + aCheckPointIt, aIndexToRemoveTo,
+            thePolygon, thePolyBoxes );
+        }
+        else
+        {
+          // Split first link
+          Standard_Integer aSplitLinkNodes[2] = {
+            aNextNodes[1],
+            aCurNodes [1]
+          };
+
+          gp_Pnt2d aSplitLinkPnts[2] = {
+            aNextPnts[1],
+            aCurPnts [1]
+          };
+          createAndReplacePolygonLink( aSplitLinkNodes, aSplitLinkPnts,
+            aPolyIt, BRepMesh_Delaun::InsertAfter, thePolygon, thePolyBoxes );
+
+          aCurNodes[1] = aNextNodes[1];
+          aCurPnts [1] = aNextPnts [1];
+          ++aIndexToRemoveTo;
+
+          processLoop( aPolyIt + 1, aIndexToRemoveTo, 
+            thePolygon, thePolyBoxes );
         }
       }
-    }
-    if ( aNextEdgeId != 0 )
-    {
-      if ( Abs( aNextEdgeId ) != theEdgeIndex && Abs( aNextEdgeId ) != aCurrEdgeId )
+      else if ( aIntFlag == BRepMesh_Delaun::Glued )
       {
-        aCurrEdgeId = Abs( aNextEdgeId );
-
-        if ( !anUsedEdges.Add( aCurrEdgeId ) )
+        if ( aCurNodes[1] == aNextNodes[0] )
         {
-          //if this edge has already been added to the aPolygon, 
-          //there is a risk of looping (attention to open contours)
-          //remove last edge and continue
-          
-          aDealLinks.Add( aCurrEdgeId );
-
-          //roll back
-          continue;
-        }        
-
-        aPolygon.Append( aNextEdgeId );
-
-        const BRepMesh_Edge& aCurrEdge = GetEdge( aCurrEdgeId );
-        Standard_Integer aVert1 = aCurrEdge.FirstNode();
-        Standard_Integer aVert2 = aCurrEdge.LastNode();
-        fillBndBox( aBoxes, GetVertex( aVert1 ), GetVertex( aVert2 ) );
+          aCurNodes[1] = aNextNodes[1];
+          aCurPnts [1] = aNextPnts [1];
+        }
+        // TODO: Non-adjacent glued links within the polygon
+      }
+      else if ( aIntFlag == BRepMesh_Delaun::Same )
+      {
+        processLoop( aPolyIt, aNextPolyIt, thePolygon, thePolyBoxes );
 
-        RemovePivotTriangles( aNextEdgeId, aPivotNode,
-          anInfectedEdges, aLoopEdges, (aPolygon.Length() == 2) );
+        isRemoveFromFirst  = Standard_True;
+        isAddReplacingEdge = Standard_False;
       }
-    }
-    else
-    {
-      //hanging end
-      if ( aPolygon.Length() == 1 )
-        return;
+      else
+        continue; // Not supported type
 
-      Standard_Integer aDeadEdgeId = Abs( aPolygon.Last() );
+      if ( isAddReplacingEdge )
+      {
+        aCurEdgeId = Abs( createAndReplacePolygonLink( aCurNodes, aCurPnts,
+          aPolyIt, BRepMesh_Delaun::Replace, thePolygon, thePolyBoxes ) );
 
-      aDealLinks.Add( aDeadEdgeId );
+        aCurEdge = &GetEdge( aCurEdgeId );
+        aCurVec  = gp_Vec2d( aCurPnts[0], aCurPnts[1] );
+      }
 
-      anUsedEdges.Remove( aDeadEdgeId );
-      aPolygon.Remove( aPolygon.Length() );
+      Standard_Integer aIndexToRemoveFrom = 
+        isRemoveFromFirst ? aPolyIt : aPolyIt + 1;
 
-      // return to previous point
-      Standard_Integer aLastValidEdge = aPolygon.Last();
-      const BRepMesh_Edge& aLastEdge = GetEdge( Abs( aLastValidEdge ) );
+      thePolygon  .Remove( aIndexToRemoveFrom, aIndexToRemoveTo );
+      thePolyBoxes.Remove( aIndexToRemoveFrom, aIndexToRemoveTo );
 
-      if( aLastValidEdge > 0 )
-      {
-        aStartNode = aLastEdge.FirstNode();
-        aPivotNode = aLastEdge.LastNode();
-      }
-      else
+      aPolyLen = thePolygon.Length();
+      if ( isRemoveFromFirst )
       {
-        aStartNode = aLastEdge.LastNode();
-        aPivotNode = aLastEdge.FirstNode();
+        --aPolyIt;
+        break;
       }
-      
-      const BRepMesh_Vertex& dV = GetVertex( aStartNode );
-      const BRepMesh_Vertex& fV = GetVertex( aPivotNode );
-      aPrevVEdge = fV.Coord() - dV.Coord();
-      continue;
+
+      aNextPolyIt = aPolyIt;
     }
-    
-    aStartNode = aPivotNode;
-    aPivotNode = aRefOtherNode;
-    aPrevVEdge = aRefCurrVEdge;
   }
 
-  CleanupPolygon( aPolygon, anInfectedEdges, aLoopEdges );
-  MeshPolygon( aPolygon );
+  meshSimplePolygon( thePolygon, thePolyBoxes );
 }
 
 //=======================================================================
-//function : MeshPolygon
-//purpose  : Triangulatiion of a closed aPolygonon described by the list of indexes of
-//           its edges in the structure. (negative index == reversed)
+//function : meshElementaryPolygon
+//purpose  : Triangulation of closed polygon containing only three edges.
 //=======================================================================
-void BRepMesh_Delaun::MeshPolygon( TColStd_SequenceOfInteger& thePoly )
+inline Standard_Boolean BRepMesh_Delaun::meshElementaryPolygon( 
+  const TColStd_SequenceOfInteger& thePolygon )
 {
-  Standard_Integer aPivotNode, aVertex3 = 0;
-  Standard_Integer aFirstNode, aLastNode;
-  Standard_Integer aTriId;
+  Standard_Integer aPolyLen = thePolygon.Length();
+  if ( aPolyLen < 3 )
+    return Standard_True;
+  else if ( aPolyLen > 3 )
+    return Standard_False;
 
-  if ( thePoly.Length() == 3 )
-  {
-    aTriId = myMeshData->AddElement( BRepMesh_Triangle( 
-      Abs( thePoly(1) ), Abs( thePoly(2) ), Abs( thePoly(3) ), 
-      thePoly(1) > 0,    thePoly(2) > 0,    thePoly(3) > 0,
-      BRepMesh_Free ) );
-      
-    myCircles.MocAdd( aTriId );
-    const BRepMesh_Edge& anEdge1 = GetEdge( Abs( thePoly(1) ) );
-    const BRepMesh_Edge& anEdge2 = GetEdge( Abs( thePoly(2) ) );
-    
-    if ( thePoly(1) > 0)
-    {
-      aFirstNode = anEdge1.FirstNode();
-      aLastNode  = anEdge1.LastNode();
-    }
-    else
-    {
-      aFirstNode = anEdge1.LastNode();
-      aLastNode  = anEdge1.FirstNode();
-    }
-    
-    if ( thePoly(2) > 0 ) 
-      aVertex3 = anEdge2.LastNode();
-    else
-      aVertex3 = anEdge2.FirstNode();
+  // Just create a triangle
+  Standard_Integer anEdges[3];
+  Standard_Boolean anEdgesOri[3];
 
-    Standard_Boolean isAdded = myCircles.Add( GetVertex( aFirstNode ).Coord(), 
-      GetVertex( aLastNode ).Coord(), GetVertex( aVertex3 ).Coord(), aTriId );
-      
-    if ( !isAdded )
-      myMeshData->RemoveElement( aTriId );
+  for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
+  {
+    Standard_Integer anEdgeInfo = thePolygon( anEdgeIt + 1 );
+    anEdges[anEdgeIt]           = Abs( anEdgeInfo );
+    anEdgesOri[anEdgeIt]        = ( anEdgeInfo > 0 );
   }
-  else if ( thePoly.Length() > 3 )
+    
+  const BRepMesh_Edge& anEdge1 = GetEdge( anEdges[0] );
+  const BRepMesh_Edge& anEdge2 = GetEdge( anEdges[1] );
+  
+  Standard_Integer aNodes[3] = { anEdge1.FirstNode(),
+                                 anEdge1.LastNode(),
+                                 anEdge2.FirstNode() };
+  if ( aNodes[2] == aNodes[0] || 
+       aNodes[2] == aNodes[1] )
   {
-    NCollection_Vector <Bnd_Box2d> aBoxes;
-    Standard_Integer aPolyIdx, aUsedIdx = 0;
-    Standard_Integer aPolyLen = thePoly.Length();
+    aNodes[2] = anEdge2.LastNode();
+  }
 
-    for ( aPolyIdx = 1; aPolyIdx <= aPolyLen; ++aPolyIdx )
-    {
-      const BRepMesh_Edge& anEdge = GetEdge( Abs( thePoly( aPolyIdx ) ) );
-      aFirstNode = anEdge.FirstNode();
-      aLastNode  = anEdge.LastNode();
-      fillBndBox( aBoxes, GetVertex( aFirstNode ), GetVertex( aLastNode ) );
-    }
-    
-    const BRepMesh_Edge& anEdge = GetEdge( Abs( thePoly(1) ) );
-    Standard_Real aMinDist = RealLast();
+  addTriangle( anEdges, anEdgesOri, aNodes );
+  return Standard_True;
+}
 
-    if ( thePoly(1) > 0 )
+//=======================================================================
+//function : meshSimplePolygon
+//purpose  : Triangulatiion of a closed simple polygon (polygon without 
+//           glued edges and loops) described by the list of indexes of 
+//           its edges in the structure.
+//           (negative index means reversed edge)
+//=======================================================================
+void BRepMesh_Delaun::meshSimplePolygon( TColStd_SequenceOfInteger&     thePolygon,
+                                         NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
+{
+  // Check is the given polygon elementary
+  if ( meshElementaryPolygon( thePolygon ) )
+    return;
+
+
+  // Polygon contains more than 3 links
+  Standard_Integer aFirstEdgeInfo = thePolygon(1);
+  const BRepMesh_Edge& aFirstEdge = GetEdge( Abs( aFirstEdgeInfo ) );
+
+  Standard_Integer aNodes[3];
+  getOrientedNodes( aFirstEdge, aFirstEdgeInfo > 0, aNodes );
+  
+  gp_Pnt2d aRefVertices[3];
+  aRefVertices[0] = GetVertex( aNodes[0] ).Coord();
+  aRefVertices[1] = GetVertex( aNodes[1] ).Coord();
+
+  gp_Vec2d aRefEdgeDir( aRefVertices[0], aRefVertices[1] );
+
+  Standard_Real aRefEdgeLen = aRefEdgeDir.Magnitude();
+  if ( aRefEdgeLen < Precision )
+    return;
+
+  aRefEdgeDir /= aRefEdgeLen;
+
+  // Find a point with minimum distance respect
+  // the end of reference link
+  Standard_Integer aUsedLinkId = 0;
+  Standard_Real    aOptAngle   = 0.0;
+  Standard_Real    aMinDist    = RealLast();
+  Standard_Integer aPivotNode  = aNodes[1];
+  Standard_Integer aPolyLen    = thePolygon.Length();
+  for ( Standard_Integer aLinkIt = 3; aLinkIt <= aPolyLen; ++aLinkIt )
+  {
+    Standard_Integer aLinkInfo = thePolygon( aLinkIt );
+    const BRepMesh_Edge& aNextEdge = GetEdge( Abs( aLinkInfo ) );
+
+    aPivotNode = aLinkInfo > 0 ? 
+      aNextEdge.FirstNode() : 
+      aNextEdge.LastNode();
+
+    gp_Pnt2d aPivotVertex = GetVertex( aPivotNode ).Coord();
+    gp_Vec2d aDistanceDir( aRefVertices[1], aPivotVertex );
+
+    Standard_Real aDist     = aRefEdgeDir ^ aDistanceDir;
+    Standard_Real aAngle    = Abs( aRefEdgeDir.Angle(aDistanceDir) );
+    Standard_Real anAbsDist = Abs( aDist );
+    if ( ( anAbsDist < Precision                 ) ||
+         ( myIsPositiveOrientation && aDist < 0. ) ||
+         (!myIsPositiveOrientation && aDist > 0. ) )
     {
-      aFirstNode = anEdge.FirstNode();
-      aLastNode  = anEdge.LastNode();
+      continue;
     }
-    else
+
+    if ( ( anAbsDist >= aMinDist                             ) &&
+         ( aAngle <= aOptAngle || aAngle > AngDeviation90Deg ) )
     {
-      aFirstNode = anEdge.LastNode();
-      aLastNode  = anEdge.FirstNode();
+      continue;
     }
-    
-    gp_XY aVEdge( GetVertex( aLastNode  ).Coord() -
-                  GetVertex( aFirstNode ).Coord() );
-                  
-    Standard_Real aVEdgeLen = aVEdge.Modulus();
-    if ( aVEdgeLen > 0.)
+
+    // Check is the test link crosses the polygon boudaries
+    Standard_Boolean isIntersect = Standard_False;
+    for ( Standard_Integer aRefLinkNodeIt = 0; aRefLinkNodeIt < 2; ++aRefLinkNodeIt )
     {
-      aVEdge.SetCoord( aVEdge.X() / aVEdgeLen,
-                       aVEdge.Y() / aVEdgeLen );
+      const Standard_Integer& aLinkFirstNode   = aNodes[aRefLinkNodeIt];
+      const gp_Pnt2d&         aLinkFirstVertex = aRefVertices[aRefLinkNodeIt];
 
-      for ( aPolyIdx = 3; aPolyIdx <= aPolyLen; ++aPolyIdx )
-      {
-        const BRepMesh_Edge& aNextEdge = GetEdge( Abs( thePoly( aPolyIdx ) ) );
-        if ( thePoly( aPolyIdx ) > 0)
-          aPivotNode = aNextEdge.FirstNode();
-        else
-          aPivotNode = aNextEdge.LastNode();
-          
-        gp_XY aVEdgePivot( GetVertex( aPivotNode ).Coord() -
-                           GetVertex( aLastNode  ).Coord() );
+      Bnd_B2d aBox;
+      aBox.Add( aLinkFirstVertex );
+      aBox.Add( aPivotVertex );
+
+      BRepMesh_Edge aCheckLink( aLinkFirstNode, aPivotNode, BRepMesh_Free );
 
-        Standard_Real aDist = aVEdge ^ aVEdgePivot;
-        if ( Abs( aDist ) > Precision::PConfusion() )
+      Standard_Integer aCheckLinkIt = 2;
+      for ( ; aCheckLinkIt <= aPolyLen; ++aCheckLinkIt )
+      {
+        if( aCheckLinkIt == aLinkIt )
+          continue;
+        
+        if ( !aBox.IsOut( thePolyBoxes.Value( aCheckLinkIt ) ) )
         {
-          if ( ( aDist > 0. &&  myPositiveOrientation ) || 
-               ( aDist < 0. && !myPositiveOrientation ) )
-          { 
-            if ( Abs( aDist ) < aMinDist )
-            {
-              Bnd_Box2d aBoxFirst, aBoxLast;
-              aBoxFirst.Add( gp_Pnt2d( GetVertex( aPivotNode ).Coord() ) );
-              aBoxFirst.Add( gp_Pnt2d( GetVertex( aLastNode  ).Coord() ) );
-              
-              aBoxLast.Add( gp_Pnt2d( GetVertex( aPivotNode ).Coord() ) );
-              aBoxLast.Add( gp_Pnt2d( GetVertex( aFirstNode ).Coord() ) );
-              
-              BRepMesh_Edge aCheckE1( aLastNode,  aPivotNode, BRepMesh_Free );
-              BRepMesh_Edge aCheckE2( aFirstNode, aPivotNode, BRepMesh_Free );
-              
-              Standard_Boolean isIntersect = Standard_False;
-              for ( Standard_Integer aPolyIdx1 = 2; aPolyIdx1 <= aPolyLen; ++aPolyIdx1 )
-              {
-                if( aPolyIdx1 == aPolyIdx )
-                  continue;
-                
-                const BRepMesh_Edge& aNextEdge1 = GetEdge( Abs( thePoly( aPolyIdx1 ) ) );
-                if ( !aBoxFirst.IsOut( aBoxes.Value( aPolyIdx1 - 1 ) ) )
-                {                    
-                  // intersection is possible...                  
-                  if( IntSegSeg( aCheckE1, aNextEdge1 ) )
-                  {
-                    isIntersect = Standard_True;
-                    break;
-                  }
-                }
-                
-                if ( !aBoxLast.IsOut( aBoxes.Value( aPolyIdx1 - 1 ) ) &&
-                     !aCheckE2.IsEqual( aNextEdge1 ) )
-                {
-                  // intersection is possible...                  
-                  if( IntSegSeg( aCheckE2, aNextEdge1 ) )
-                  {
-                    isIntersect = Standard_True;
-                    break;
-                  }
-                }
-              }
-              
-              if( !isIntersect )
-              {
-                aMinDist = aDist;
-                aVertex3 = aPivotNode;
-                aUsedIdx = aPolyIdx;
-              }
-            }
+          const BRepMesh_Edge& aPolyLink = 
+            GetEdge( Abs( thePolygon( aCheckLinkIt ) ) );
+
+          if ( aCheckLink.IsEqual( aPolyLink ) )
+            continue;
+
+          // intersection is possible...                  
+          gp_Pnt2d anIntPnt;
+          IntFlag aIntFlag = intSegSeg( aCheckLink, aPolyLink, 
+            Standard_False, Standard_False, anIntPnt );
+
+          if( aIntFlag != BRepMesh_Delaun::NoIntersection )
+          {
+            isIntersect = Standard_True;
+            break;
           }
         }
       }
+
+      if ( isIntersect )
+        break;
     }
 
-    Standard_Integer aNewEdge2, aNewEdge3;
-    if ( aMinDist < RealLast() )
-    {
-      aNewEdge2 = myMeshData->AddLink( BRepMesh_Edge( aLastNode, aVertex3,   BRepMesh_Free ) );
-      aNewEdge3 = myMeshData->AddLink( BRepMesh_Edge( aVertex3,  aFirstNode, BRepMesh_Free ) );
-      aTriId    = myMeshData->AddElement( BRepMesh_Triangle( 
-        Abs( thePoly(1) ), Abs( aNewEdge2 ), Abs( aNewEdge3 ), 
-        thePoly(1) > 0,    aNewEdge2 > 0,    aNewEdge3 > 0,
-        BRepMesh_Free ) );
-
-      Standard_Boolean isAdded = myCircles.Add( GetVertex( aFirstNode ).Coord(), 
-        GetVertex( aLastNode ).Coord(), GetVertex( aVertex3 ).Coord(), aTriId );
-        
-      if ( !isAdded )
-        myMeshData->RemoveElement( aTriId );
+    if( isIntersect )
+      continue;
 
-      if ( aUsedIdx < aPolyLen )
-      {
-        TColStd_SequenceOfInteger aSuitePoly;
-        thePoly.Split( aUsedIdx, aSuitePoly );
-        aSuitePoly.Prepend( -aNewEdge3 );
-        MeshPolygon( aSuitePoly );
-      }
-      else 
-        thePoly.Remove( aPolyLen );
 
-      if ( aUsedIdx > 3 )
-      {
-        thePoly.SetValue( 1, -aNewEdge2 );
-        MeshPolygon( thePoly );
-      }
-    }
+    aOptAngle       = aAngle;
+    aMinDist        = anAbsDist;
+    aNodes[2]       = aPivotNode;
+    aRefVertices[2] = aPivotVertex;
+    aUsedLinkId     = aLinkIt;
+  }
+
+  if ( aUsedLinkId == 0 )
+    return;
+
+
+  BRepMesh_Edge aNewEdges[2] = {
+    BRepMesh_Edge( aNodes[1], aNodes[2], BRepMesh_Free ),
+    BRepMesh_Edge( aNodes[2], aNodes[0], BRepMesh_Free ) };
+
+  Standard_Integer aNewEdgesInfo[3] = {
+    aFirstEdgeInfo,
+    myMeshData->AddLink( aNewEdges[0] ),
+    myMeshData->AddLink( aNewEdges[1] ) };
+
+
+  Standard_Integer anEdges[3];
+  Standard_Boolean anEdgesOri[3];
+  for ( Standard_Integer aTriEdgeIt = 0; aTriEdgeIt < 3; ++aTriEdgeIt )
+  {
+    const Standard_Integer& anEdgeInfo = aNewEdgesInfo[aTriEdgeIt];
+    anEdges[aTriEdgeIt]    = Abs( anEdgeInfo );
+    anEdgesOri[aTriEdgeIt] = anEdgeInfo > 0;
+  }
+  addTriangle( anEdges, anEdgesOri, aNodes );
+
+  // Create triangle and split the source polygon on two 
+  // parts (if possible) and mesh each part as independent
+  // polygon.
+  if ( aUsedLinkId < aPolyLen )
+  {
+    TColStd_SequenceOfInteger aRightPolygon;
+    thePolygon.Split( aUsedLinkId, aRightPolygon );
+    aRightPolygon.Prepend( -aNewEdgesInfo[2] );
+
+    NCollection_Sequence<Bnd_B2d> aRightPolyBoxes;
+    thePolyBoxes.Split( aUsedLinkId, aRightPolyBoxes );
+
+    Bnd_B2d aBox;
+    aBox.Add( aRefVertices[0] );
+    aBox.Add( aRefVertices[2] );
+    aRightPolyBoxes.Prepend( aBox );
+
+    meshSimplePolygon( aRightPolygon, aRightPolyBoxes );
+  }
+  else
+  {
+    thePolygon.Remove  ( aPolyLen );
+    thePolyBoxes.Remove( aPolyLen );
+  }
+
+  if ( aUsedLinkId > 3 )
+  {
+    thePolygon.SetValue( 1, -aNewEdgesInfo[1] );
+
+    Bnd_B2d aBox;
+    aBox.Add( aRefVertices[1] );
+    aBox.Add( aRefVertices[2] );
+
+    thePolyBoxes.SetValue( 1, aBox );
+
+    meshSimplePolygon( thePolygon, thePolyBoxes );
   }
 }
 
@@ -1363,9 +1951,10 @@ void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
   // Loop on triangles to be destroyed :
   BRepMesh_MapOfInteger::Iterator aTriangleIt( aSelector.Elements() );
   for ( ; aTriangleIt.More(); aTriangleIt.Next() )
-    DeleteTriangle( aTriangleIt.Key(), aLoopEdges );
+    deleteTriangle( aTriangleIt.Key(), aLoopEdges );
 
-  TColStd_SequenceOfInteger aPolygon;
+  NCollection_Sequence<Bnd_B2d> aBoxes;
+  TColStd_SequenceOfInteger     aPolygon;
   Standard_Integer aLoopEdgesCount = aLoopEdges.Extent();
   BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
 
@@ -1390,6 +1979,8 @@ void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
     else
       aPolygon.Append( anEdgeId );
 
+    fillBndBox( aBoxes, GetVertex( aFirstNode ), GetVertex( aPivotNode ) );
+
     aLoopEdges.UnBind( anEdgeId );
     
     aLastNode = aFirstNode;
@@ -1413,6 +2004,8 @@ void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
           }
           else
             aPolygon.Append( anEdgeId );
+
+          fillBndBox( aBoxes, GetVertex( aCurrentNode ), GetVertex( aPivotNode ) );
             
           aPivotNode = aCurrentNode;
           aLoopEdges.UnBind( anEdgeId );
@@ -1425,7 +2018,7 @@ void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
       --aLoopEdgesCount;
     }
     
-    MeshPolygon( aPolygon );
+    meshPolygon( aPolygon, aBoxes );
   }
 }
 
@@ -1437,7 +2030,7 @@ void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
 void  BRepMesh_Delaun::AddVertices( BRepMesh_Array1OfVertexOfDelaun& theVertices )
 {
   BRepMesh_HeapSortVertexOfDelaun::Sort( theVertices, 
-    BRepMesh_ComparatorOfVertexOfDelaun( SortingDirection, Precision::PConfusion() ) );
+    BRepMesh_ComparatorOfVertexOfDelaun( SortingDirection, Precision ) );
 
   Standard_Integer aLower  = theVertices.Lower();
   Standard_Integer anUpper = theVertices.Upper();
@@ -1446,14 +2039,14 @@ void  BRepMesh_Delaun::AddVertices( BRepMesh_Array1OfVertexOfDelaun& theVertices
   for ( Standard_Integer i = aLower; i <= anUpper; ++i )     
     aVertexIndexes(i) = myMeshData->AddNode( theVertices(i) );
 
-  CreateTrianglesOnNewVertices( aVertexIndexes );
+  createTrianglesOnNewVertices( aVertexIndexes );
 }
 
 //=======================================================================
 //function : UseEdge
 //purpose  : Modify mesh to use the edge. Return True if done
 //=======================================================================
-Standard_Boolean BRepMesh_Delaun::UseEdge( const Standard_Integer /*theIndex*/ )
+Standard_Boolean BRepMesh_Delaun::UseEdge( const Standard_Integer theIndex )
 {
   /*
   const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( theIndex );
@@ -1538,69 +2131,27 @@ Standard_Boolean BRepMesh_Delaun::UseEdge( const Standard_Integer /*theIndex*/ )
 }
 
 //=======================================================================
-//function : Result
-//purpose  : Gives the Mesh data structure
+//function : getEdgesByType
+//purpose  : Gives the list of edges with type defined by input parameter
 //=======================================================================
-const Handle( BRepMesh_DataStructureOfDelaun )& BRepMesh_Delaun::Result() const
-{
-  return myMeshData;
-}
-
-//=======================================================================
-//function : Frontier
-//purpose  : Gives the list of frontier edges
-//=======================================================================
-const BRepMesh_MapOfInteger& BRepMesh_Delaun::Frontier()
+Handle(BRepMesh_MapOfInteger) BRepMesh_Delaun::getEdgesByType(
+  const BRepMesh_DegreeOfFreedom theEdgeType ) const
 {
+  Handle(BRepMesh_MapOfInteger) aResult = new BRepMesh_MapOfInteger;
   BRepMesh_MapOfInteger::Iterator anEdgeIt( myMeshData->LinkOfDomain() );
 
-  myMapEdges.Clear();
   for ( ; anEdgeIt.More(); anEdgeIt.Next() )
   {
     Standard_Integer anEdge = anEdgeIt.Key();
-    if ( GetEdge( anEdge ).Movability() == BRepMesh_Frontier )
-      myMapEdges.Add( anEdge );
-  }
-  
-  return myMapEdges;
-}
-
-//=======================================================================
-//function : InternalEdges
-//purpose  : Gives the list of internal edges
-//=======================================================================
-const BRepMesh_MapOfInteger& BRepMesh_Delaun::InternalEdges()
-{
-  BRepMesh_MapOfInteger::Iterator anEdgeIt( myMeshData->LinkOfDomain() );
+    Standard_Boolean isToAdd = (theEdgeType == BRepMesh_Free) ?
+      (myMeshData->ElemConnectedTo( anEdge ).Extent() <= 1) :
+      (GetEdge( anEdge ).Movability() == theEdgeType);
 
-  myMapEdges.Clear();
-  for ( ; anEdgeIt.More(); anEdgeIt.Next() )
-  {
-    Standard_Integer anEdge = anEdgeIt.Key();
-    if ( GetEdge( anEdge ).Movability() == BRepMesh_Fixed )
-      myMapEdges.Add( anEdge );
+    if (isToAdd)
+      aResult->Add( anEdge );
   }
   
-  return myMapEdges;
-}
-
-//=======================================================================
-//function : FreeEdges
-//purpose  : Gives the list of free edges used only one time
-//=======================================================================
-const BRepMesh_MapOfInteger& BRepMesh_Delaun::FreeEdges()
-{
-  BRepMesh_MapOfInteger::Iterator anEdgeIt( myMeshData->LinkOfDomain() );
-
-  myMapEdges.Clear();
-  for ( ; anEdgeIt.More(); anEdgeIt.Next() )
-  {
-    Standard_Integer anEdge = anEdgeIt.Key();
-    if ( myMeshData->ElemConnectedTo( anEdge ).Extent() <= 1 )
-      myMapEdges.Add( anEdge );
-  }
-  
-  return myMapEdges;
+  return aResult;
 }
 
 //=======================================================================
@@ -1608,13 +2159,13 @@ const BRepMesh_MapOfInteger& BRepMesh_Delaun::FreeEdges()
 //purpose  : Calculates distances between the given point and edges of
 //           triangle
 //=======================================================================
-static Standard_Real calculateDist( const gp_XY theVEdges[3],
-                                    const gp_XY thePoints[3],
-                                    const Standard_Integer theEdgesId[3],
-                                    const BRepMesh_Vertex& theVertex,
-                                    Standard_Real theDistance[3],
-                                    Standard_Real theSqModulus[3],
-                                    Standard_Integer& theEdgeOn )
+Standard_Real BRepMesh_Delaun::calculateDist( const gp_XY            theVEdges[3],
+                                              const gp_XY            thePoints[3],
+                                              const Standard_Integer theEdgesId[3],
+                                              const BRepMesh_Vertex& theVertex,
+                                              Standard_Real          theDistance[3],
+                                              Standard_Real          theSqModulus[3],
+                                              Standard_Integer&      theEdgeOn ) const
 {
   Standard_Real aMinDist = -1;
   if ( !theVEdges || !thePoints || !theEdgesId || 
@@ -1624,7 +2175,7 @@ static Standard_Real calculateDist( const gp_XY theVEdges[3],
   for( Standard_Integer i = 0; i < 3; ++i )
   {
     theSqModulus[i] = theVEdges[i].SquareModulus();
-    if ( theSqModulus[i] <= EPSEPS )
+    if ( theSqModulus[i] <= Precision2 )
       return -1;
       
     theDistance[i] = theVEdges[i] ^ ( theVertex.Coord() - thePoints[i] );
@@ -1650,7 +2201,7 @@ static Standard_Real calculateDist( const gp_XY theVEdges[3],
 //=======================================================================
 Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId,
                                             const BRepMesh_Vertex& theVertex,
-                                            Standard_Integer& theEdgeOn ) const
+                                            Standard_Integer&      theEdgeOn ) const
 {
   theEdgeOn = 0;
   
@@ -1703,7 +2254,7 @@ Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId
   if ( aMinDist < 0 )
     return Standard_False;
       
-  if ( aMinDist > EPSEPS )
+  if ( aMinDist > Precision2 )
   {
     Standard_Integer anEdgeId = theEdgeOn;
     theEdgeOn = 0;
@@ -1730,28 +2281,26 @@ Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId
 
 
 //=============================================================================
-// Function: classifyPoint
-// This function is used for point classifying in case of coincidence of two 
-// vectors.
-// Returns zero value if point is out of segment and non zero value if point is
-// between the first and the second point of segment.
-//
-// thePoint1       - the start point of a segment (base point)
-// thePoint2       - the end point of a segment
-// thePointToCheck - the point to classify
+//function : classifyPoint
+//purpose  : Classifies the point in case of coincidence of two vectors.
+//           Returns zero value if point is out of segment and non zero 
+//           value if point is between the first and the second point of segment.
+//           thePoint1       - the start point of a segment (base point)
+//           thePoint2       - the end point of a segment
+//           thePointToCheck - the point to classify
 //=============================================================================
-static Standard_Integer classifyPoint( const gp_XY& thePoint1,
-                                       const gp_XY& thePoint2,
-                                       const gp_XY& thePointToCheck )
+Standard_Integer BRepMesh_Delaun::classifyPoint( const gp_XY& thePoint1,
+                                                 const gp_XY& thePoint2,
+                                                 const gp_XY& thePointToCheck ) const
 {
   gp_XY aP1 = thePoint2       - thePoint1;
   gp_XY aP2 = thePointToCheck - thePoint1;
   
   Standard_Real aDist = Abs( aP1 ^ aP2 );
-  if ( aDist >= Precision::PConfusion() )
+  if ( aDist >= Precision )
   {
     aDist = ( aDist * aDist ) / aP1.SquareModulus();
-    if ( aDist >= EPSEPS )
+    if ( aDist >= Precision2 )
       return 0; //out
   }
     
@@ -1762,18 +2311,22 @@ static Standard_Integer classifyPoint( const gp_XY& thePoint1,
   if ( aP1.SquareModulus() < aP2.SquareModulus() )
     return 0; //out
     
-  if ( thePointToCheck.IsEqual( thePoint1, Precision::PConfusion() ) || 
-       thePointToCheck.IsEqual( thePoint2, Precision::PConfusion() ) )
-    return -1; //end point
+  if ( thePointToCheck.IsEqual( thePoint1, Precision ) || 
+       thePointToCheck.IsEqual( thePoint2, Precision ) )
+    return -1; //coinsides with an end point
     
   return 1;
 }
 
 //=============================================================================
-// Function: IntSegSeg
+//function : intSegSeg
+//purpose  : Checks intersection between the two segments.
 //=============================================================================
-Standard_Boolean BRepMesh_Delaun::IntSegSeg( const BRepMesh_Edge& theEdg1,
-                                             const BRepMesh_Edge& theEdg2 )
+BRepMesh_Delaun::IntFlag BRepMesh_Delaun::intSegSeg( const BRepMesh_Edge&   theEdg1,
+                                                     const BRepMesh_Edge&   theEdg2,
+                                                     const Standard_Boolean isConsiderEndPointTouch,
+                                                     const Standard_Boolean isConsiderPointOnEdge,
+                                                     gp_Pnt2d&              theIntPnt) const
 {
   gp_XY p1, p2, p3, p4;
   p1 = GetVertex( theEdg1.FirstNode() ).Coord();
@@ -1785,42 +2338,152 @@ Standard_Boolean BRepMesh_Delaun::IntSegSeg( const BRepMesh_Edge& theEdg1,
   Standard_Integer aPoint2 = classifyPoint( p1, p2, p4 );
   Standard_Integer aPoint3 = classifyPoint( p3, p4, p1 );
   Standard_Integer aPoint4 = classifyPoint( p3, p4, p2 );
+
+  // Consider case when edges have shared vertex
+  if ( isConsiderEndPointTouch )
+  {
+    if ( aPoint1 < 0 || aPoint2 < 0 )
+      return BRepMesh_Delaun::EndPointTouch;
+  }
+
+  Standard_Integer aPosHash = 
+    aPoint1 + aPoint2 + aPoint3 + aPoint4;
+
+  /*=========================================*/
+  /*  1) hash code == 1:
+
+                    0+
+                    /
+           0      1/         0
+           +======+==========+
   
-  if ( aPoint1 > 0 || aPoint2 > 0 ||
-       aPoint3 > 0 || aPoint4 > 0 )
-    return Standard_True;
-                    
-  gp_XY aPl1 = p2 - p1;
-  gp_XY aPl2 = p4 - p3;
-  gp_XY aPl3 = p1 - p3;
+      2) hash code == 2:
+
+           0    1        1   0
+        a) +----+========+---+
+
+           0       1   1     0
+        b) +-------+===+=====+
+
+                                             */
+  /*=========================================*/
+  if ( aPosHash == 1 )
+  {
+    return isConsiderPointOnEdge ? 
+      BRepMesh_Delaun::PointOnEdge :
+      BRepMesh_Delaun::NoIntersection;
+  }
+  else if ( aPosHash == 2 )
+    return BRepMesh_Delaun::Glued;
+
+  gp_XY aVec1           = p2 - p1;
+  gp_XY aVec2           = p4 - p3;
+  gp_XY aVecStartPoints = p3 - p1;
     
-  Standard_Real aCrossD1D2 = aPl1 ^ aPl2;
-  Standard_Real aCrossD1D3 = aPl1 ^ aPl3;
-  if ( Abs( aCrossD1D2 ) < Precision::PConfusion() )
+  Standard_Real aCrossD1D2 = aVec1           ^ aVec2;
+  Standard_Real aCrossD1D3 = aVecStartPoints ^ aVec2;
+
+  // is edgegs codirectional
+  if ( Abs( aCrossD1D2 ) < Precision )
   {
-    if( Abs( aCrossD1D3 ) < Precision::PConfusion() )
+    // just a parallel case?
+    if( Abs( aCrossD1D3 ) < Precision )
     {
-      Standard_Integer aPosHash = aPoint1 + aPoint2;
-      if ( ( !aPosHash && aPoint3 ) ) //|| aPosHash < -1 )
-        return Standard_True;
-        
-      return Standard_False;
+      /*=========================================*/
+      /*  Here the following cases are possible:
+          1) hash code == -4:
+
+               -1                -1
+                +=================+
+               -1                -1
+
+          2) hash code == -2:
+
+                0       -1        0
+                +--------+========+
+                        -1
+
+          3) hash code == -1:
+
+                0        1        -1
+                +--------+========+
+                                  -1
+
+          4) hash code == 0:
+
+                0      0  0       0
+                +------+  +=======+
+                0      0  0       0
+                                                 */
+      /*=========================================*/
+
+      if ( aPosHash < -2 )
+        return BRepMesh_Delaun::Same;
+      else if ( aPosHash == -1 )
+        return BRepMesh_Delaun::Glued;
+
+      return BRepMesh_Delaun::NoIntersection;
     }
     else
-      //parallel case
-      return Standard_False;
+      return BRepMesh_Delaun::NoIntersection;
   }
 
   Standard_Real aPar = aCrossD1D3 / aCrossD1D2;
   // inrersects out of first segment range
-  if( aPar < Precision::Angular() || aPar > 1 - Precision::Angular() )
-    return Standard_False;
+  if( aPar < Precision || aPar > EndPrecision )
+    return BRepMesh_Delaun::NoIntersection;
  
-  Standard_Real aCrossD2D3 = aPl2 ^ aPl3;
-  aPar = aCrossD2D3 / aCrossD1D2;
+  Standard_Real aCrossD2D3 = aVecStartPoints.Reversed() ^ aVec1;
+  aPar = aCrossD2D3 / -aCrossD1D2;
   // inrersects out of second segment range
-  if( aPar < Precision::Angular() || aPar > 1 - Precision::Angular() )
-    return Standard_False;
+  if( aPar < Precision || aPar > EndPrecision )
+    return BRepMesh_Delaun::NoIntersection;
  
-  return Standard_True;
+  theIntPnt = p3 + aPar * aVec2;
+  return BRepMesh_Delaun::Cross;
+}
+
+//=============================================================================
+//function : polyArea
+//purpose  : Returns area of the loop of the given polygon defined by indices 
+//           of its start and end links.
+//=============================================================================
+Standard_Real BRepMesh_Delaun::polyArea( const TColStd_SequenceOfInteger& thePolygon,
+                                         const Standard_Integer           theStartIndex,
+                                         const Standard_Integer           theEndIndex ) const
+{
+  Standard_Real aArea = 0.0;
+  Standard_Integer aPolyLen = thePolygon.Length();
+  if ( theStartIndex >= theEndIndex ||
+       theStartIndex > aPolyLen )
+  {
+    return aArea;
+  }
+
+  Standard_Integer aEndIndex = (theEndIndex > aPolyLen) ? 
+    aPolyLen : theEndIndex;
+
+  Standard_Integer aCurEdgeInfo = thePolygon( theStartIndex );
+  Standard_Integer aCurEdgeId   = Abs( aCurEdgeInfo );
+  const BRepMesh_Edge* aCurEdge = &GetEdge( aCurEdgeId );
+
+  Standard_Integer aNodes[2];
+  getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aNodes );
+
+  gp_Pnt2d aRefPnt = GetVertex( aNodes[0] ).Coord();
+  Standard_Integer aPolyIt = theStartIndex + 1;
+  for ( ; aPolyIt <= theEndIndex; ++aPolyIt )
+  {
+    aCurEdgeInfo = thePolygon( aPolyIt );
+    aCurEdgeId   = Abs( aCurEdgeInfo );
+    aCurEdge     = &GetEdge( aCurEdgeId );
+
+    getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aNodes );
+    gp_Vec2d aVec1( aRefPnt, GetVertex( aNodes[0] ).Coord() );
+    gp_Vec2d aVec2( aRefPnt, GetVertex( aNodes[1] ).Coord() );
+
+    aArea += aVec1 ^ aVec2;
+  }
+
+  return aArea / 2.;
 }