0023105: Exception during Meshing / Missing triangles
[occt.git] / src / BRepMesh / BRepMesh_Delaun.cxx
1 // Created on: 1993-05-12
2 // Created by: Didier PIFFAULT
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #include <BRepMesh_Delaun.hxx>
23
24 #include <gp.hxx>
25 #include <gp_XY.hxx>
26 #include <gp_Pnt2d.hxx>
27 #include <gp_Vec2d.hxx>
28
29 #include <Precision.hxx>
30 #include <Bnd_Box2d.hxx>
31 #include <Bnd_B2d.hxx>
32
33 #include <TColStd_MapOfInteger.hxx>
34 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
35 #include <TColStd_Array1OfBoolean.hxx>
36 #include <TColStd_ListIteratorOfListOfInteger.hxx>
37 #include <TColStd_ListOfInteger.hxx>
38 #include <TColStd_Array1OfReal.hxx>
39 #include <TColStd_Array1OfInteger.hxx>
40 #include <TColStd_SequenceOfInteger.hxx>
41
42 #include <BRepMesh_MapOfIntegerInteger.hxx>
43 #include <BRepMesh_HeapSortIndexedVertexOfDelaun.hxx>
44 #include <BRepMesh_ComparatorOfIndexedVertexOfDelaun.hxx>
45 #include <BRepMesh_HeapSortIndexedVertexOfDelaun.hxx>
46 #include <BRepMesh_SelectorOfDataStructureOfDelaun.hxx>
47 #include <BRepMesh_HeapSortVertexOfDelaun.hxx>
48 #include <BRepMesh_ComparatorOfVertexOfDelaun.hxx>
49 #include <BRepMesh_Array1OfVertexOfDelaun.hxx>
50
51 #include <BRepMesh_Edge.hxx>
52 #include <BRepMesh_Vertex.hxx>
53 #include <BRepMesh_Triangle.hxx>
54
55 #include <NCollection_Vector.hxx>
56
57 typedef TColStd_ListIteratorOfListOfInteger  IteratorOnListOfInteger;
58 typedef TColStd_ListOfInteger                ListOfInteger;
59
60 const Standard_Real AngDeviation1Deg  = M_PI/180.;
61 const Standard_Real AngDeviation90Deg = 90 * AngDeviation1Deg;
62 const Standard_Real Angle2PI          = 2 * M_PI;
63
64 const Standard_Real Precision    = Precision::PConfusion();
65 const Standard_Real EndPrecision = 1 - Precision;
66 const Standard_Real Precision2   = Precision * Precision;
67 const gp_XY SortingDirection(M_SQRT1_2, M_SQRT1_2);
68
69 //=======================================================================
70 //function : BRepMesh_Delaun
71 //purpose  : Creates the triangulation with an empty Mesh data structure
72 //=======================================================================
73 BRepMesh_Delaun::BRepMesh_Delaun( BRepMesh_Array1OfVertexOfDelaun& theVertices,
74                                   const Standard_Boolean           isPositive )
75 : myIsPositiveOrientation( isPositive ),
76   myCircles( theVertices.Length(), new NCollection_IncAllocator() )
77 {
78   if ( theVertices.Length() > 2 )
79   {
80     myMeshData = new BRepMesh_DataStructureOfDelaun( new NCollection_IncAllocator(),
81                                                      theVertices.Length() );
82     Init( theVertices );
83   }
84 }
85
86 //=======================================================================
87 //function : BRepMesh_Delaun
88 //purpose  : Creates the triangulation with and existent Mesh data structure
89 //=======================================================================
90 BRepMesh_Delaun::BRepMesh_Delaun( const Handle( BRepMesh_DataStructureOfDelaun )& theOldMesh,
91                                   BRepMesh_Array1OfVertexOfDelaun&                theVertices,
92                                   const Standard_Boolean                          isPositive )
93  : myIsPositiveOrientation( isPositive ),
94    myCircles( theVertices.Length(), theOldMesh->Allocator() )
95 {
96   myMeshData = theOldMesh;
97   if ( theVertices.Length() > 2 )
98     Init( theVertices );
99 }
100
101 //=======================================================================
102 //function : BRepMesh_Delaun
103 //purpose  : Creates the triangulation with and existent Mesh data structure
104 //=======================================================================
105 BRepMesh_Delaun::BRepMesh_Delaun( const Handle( BRepMesh_DataStructureOfDelaun )& theOldMesh, 
106                                   TColStd_Array1OfInteger&                        theVertexIndices,
107                                   const Standard_Boolean                          isPositive )
108  : myIsPositiveOrientation( isPositive ),
109    myCircles( theVertexIndices.Length(), theOldMesh->Allocator() )
110 {
111   myMeshData = theOldMesh;
112   if ( theVertexIndices.Length() > 2 )
113   {
114     Bnd_Box2d aBox;
115     Standard_Integer anIndex = theVertexIndices.Lower();
116     Standard_Integer anUpper = theVertexIndices.Upper();
117     for ( ; anIndex <= anUpper; ++anIndex )
118       aBox.Add( gp_Pnt2d( GetVertex( theVertexIndices( anIndex) ).Coord() ) );
119
120     perform( aBox, theVertexIndices );
121   }
122 }
123
124 //=======================================================================
125 //function : Init
126 //purpose  : Initializes the triangulation with an Array of Vertex
127 //=======================================================================
128 void BRepMesh_Delaun::Init( BRepMesh_Array1OfVertexOfDelaun& theVertices )
129 {
130   Bnd_Box2d aBox;
131   Standard_Integer aLowerIdx  = theVertices.Lower();
132   Standard_Integer anUpperIdx = theVertices.Upper();
133   TColStd_Array1OfInteger aVertexIndexes( aLowerIdx, anUpperIdx );
134   
135   Standard_Integer anIndex = aLowerIdx;
136   for ( ; anIndex <= anUpperIdx; ++anIndex )
137   {
138     aBox.Add( gp_Pnt2d( theVertices( anIndex ).Coord() ) );
139     aVertexIndexes( anIndex ) = myMeshData->AddNode( theVertices( anIndex ) );
140   }
141
142   perform( aBox, aVertexIndexes );
143 }
144
145 //=======================================================================
146 //function : perform
147 //purpose  : Create super mesh and run triangulation procedure
148 //=======================================================================
149 void BRepMesh_Delaun::perform( Bnd_Box2d&               theBndBox,
150                                TColStd_Array1OfInteger& theVertexIndexes )
151 {
152   theBndBox.Enlarge( Precision );
153   superMesh( theBndBox );
154
155   BRepMesh_HeapSortIndexedVertexOfDelaun::Sort( theVertexIndexes, 
156       BRepMesh_ComparatorOfIndexedVertexOfDelaun( SortingDirection,
157         Precision, myMeshData ) );
158
159   compute( theVertexIndexes );
160 }
161
162 //=======================================================================
163 //function : superMesh
164 //purpose  : Build the super mesh
165 //=======================================================================
166 void BRepMesh_Delaun::superMesh( const Bnd_Box2d& theBox )
167 {
168   Standard_Real aMinX, aMinY, aMaxX, aMaxY;
169   theBox.Get( aMinX, aMinY, aMaxX, aMaxY );
170   Standard_Real aDeltaX = aMaxX - aMinX;
171   Standard_Real aDeltaY = aMaxY - aMinY;
172
173   Standard_Real aDeltaMin = Min( aDeltaX, aDeltaY );
174   Standard_Real aDeltaMax = Max( aDeltaX, aDeltaY );
175   Standard_Real aDelta    = aDeltaX + aDeltaY;
176   
177   myCircles.SetMinMaxSize( gp_XY( aMinX, aMinY ), gp_XY( aMaxX, aMaxY ) );
178
179   Standard_Integer aScaler = 2;
180   if ( myMeshData->NbNodes() > 100 )
181     aScaler = 5;
182   else if( myMeshData->NbNodes() > 1000 )
183     aScaler = 7;
184
185   myCircles.SetCellSize( aDeltaX / aScaler,
186                          aDeltaY / aScaler );
187
188   mySupVert[0] = myMeshData->AddNode(
189     BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) );
190     
191   mySupVert[1] = myMeshData->AddNode(
192     BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
193     
194   mySupVert[2] = myMeshData->AddNode(
195     BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
196
197   if ( !myIsPositiveOrientation )
198   {
199     Standard_Integer aTmp;
200     aTmp         = mySupVert[1];
201     mySupVert[1] = mySupVert[2];
202     mySupVert[2] = aTmp;
203   }
204
205   Standard_Integer anEdgeId[3];
206   
207   for (Standard_Integer aNodeId = 0; aNodeId < 3; ++aNodeId)
208   {
209     Standard_Integer aFirstNode = aNodeId;
210     Standard_Integer aLastNode  = (aNodeId + 1) % 3;
211     anEdgeId[aNodeId] = myMeshData->AddLink( BRepMesh_Edge( 
212       mySupVert[aFirstNode], mySupVert[aLastNode], BRepMesh_Free ) );
213   }
214   
215   mySupTrian = BRepMesh_Triangle( 
216     Abs( anEdgeId[0] ),  Abs( anEdgeId[1] ),  Abs( anEdgeId[2] ), 
217     ( anEdgeId[0] > 0 ), ( anEdgeId[1] > 0 ), ( anEdgeId[2] > 0),
218     BRepMesh_Free);
219 }
220
221 //=======================================================================
222 //function : deleteTriangle
223 //purpose  : Deletes the triangle with the given index and adds the free
224 //           edges into the map.
225 //           When an edge is suppressed more than one time it is destroyed.
226 //=======================================================================
227 void  BRepMesh_Delaun::deleteTriangle( const Standard_Integer        theIndex, 
228                                        BRepMesh_MapOfIntegerInteger& theLoopEdges )
229 {
230   myCircles.Delete( theIndex );
231
232   Standard_Integer e[3];
233   Standard_Boolean o[3];
234   GetTriangle( theIndex ).Edges( e[0], e[1], e[2],
235                                  o[0], o[1], o[2] );
236   
237   myMeshData->RemoveElement( theIndex );
238
239   for ( Standard_Integer i = 0; i < 3; ++i )
240   {
241     if ( !theLoopEdges.Bind( e[i], o[i] ) ) 
242     {
243       theLoopEdges.UnBind( e[i] );
244       myMeshData->RemoveLink( e[i] );
245     }
246   }
247 }
248
249 //=======================================================================
250 //function : compute
251 //purpose  : Computes the triangulation and add the vertices edges and 
252 //           triangles to the Mesh data structure
253 //=======================================================================
254 void BRepMesh_Delaun::compute( TColStd_Array1OfInteger& theVertexIndexes )
255 {
256   // Insertion of edges of super triangles in the list of free edges: 
257   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
258   Standard_Integer e[3];
259   Standard_Boolean o[3];
260   mySupTrian.Edges( e[0], e[1], e[2],
261                     o[0], o[1], o[2] );
262                     
263   aLoopEdges.Bind( e[0], Standard_True );
264   aLoopEdges.Bind( e[1], Standard_True );
265   aLoopEdges.Bind( e[2], Standard_True );
266
267   if ( theVertexIndexes.Length() > 0 )
268   {
269     // Creation of 3 trianglers with the first node and the edges of the super triangle:
270     Standard_Integer anVertexIdx = theVertexIndexes.Lower();
271     createTriangles( theVertexIndexes( anVertexIdx ), aLoopEdges );
272
273     // Add other nodes to the mesh
274     createTrianglesOnNewVertices( theVertexIndexes );
275   }
276
277   // Destruction of triangles containing a top of the super triangle
278   BRepMesh_SelectorOfDataStructureOfDelaun aSelector( myMeshData );
279   for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
280     aSelector.NeighboursOfNode( mySupVert[aSupVertId] );
281   
282   aLoopEdges.Clear();
283   BRepMesh_MapOfInteger::Iterator aFreeTriangles( aSelector.Elements() );
284   for ( ; aFreeTriangles.More(); aFreeTriangles.Next() )
285     deleteTriangle( aFreeTriangles.Key(), aLoopEdges );
286
287   // All edges that remain free are removed from aLoopEdges;
288   // only the boundary edges of the triangulation remain there
289   BRepMesh_MapOfIntegerInteger::Iterator aFreeEdges( aLoopEdges );
290   for ( ; aFreeEdges.More(); aFreeEdges.Next() )
291   {
292     if ( myMeshData->ElemConnectedTo( aFreeEdges.Key() ).IsEmpty() )
293       myMeshData->RemoveLink( aFreeEdges.Key() );
294   }
295
296   // The tops of the super triangle are destroyed
297   for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
298     myMeshData->RemoveNode( mySupVert[aSupVertId] );
299 }
300
301 //=======================================================================
302 //function : createTriangles
303 //purpose  : Creates the triangles beetween the node and the polyline.
304 //=======================================================================
305 void BRepMesh_Delaun::createTriangles ( const Standard_Integer        theVertexIndex,  
306                                         BRepMesh_MapOfIntegerInteger& thePoly )
307 {
308   ListOfInteger aLoopEdges, anExternalEdges;
309   const gp_XY& aVertexCoord = myMeshData->GetNode( theVertexIndex ).Coord();
310   
311   BRepMesh_MapOfIntegerInteger::Iterator anEdges( thePoly );
312   for ( ; anEdges.More(); anEdges.Next() )
313   {
314     Standard_Integer     anEdgeId = anEdges.Key();
315     const BRepMesh_Edge& anEdge   = GetEdge( anEdgeId );
316
317     Standard_Boolean isPositive = (Standard_Boolean)thePoly( anEdgeId );
318
319     Standard_Integer aNodes[3];
320     if ( isPositive )
321     {
322       aNodes[0] = anEdge.FirstNode();
323       aNodes[2] = anEdge.LastNode();
324     }
325     else
326     {
327       aNodes[0] = anEdge.LastNode();
328       aNodes[2] = anEdge.FirstNode();
329     }
330     aNodes[1] = theVertexIndex;
331
332     const BRepMesh_Vertex& aFirstVertex = GetVertex( aNodes[0] );
333     const BRepMesh_Vertex& aLastVertex  = GetVertex( aNodes[2] );
334
335     gp_XY anEdgeDir( aLastVertex.Coord() - aFirstVertex.Coord() );
336     Standard_Real anEdgeLen = anEdgeDir.Modulus();
337     if ( anEdgeLen < Precision )
338       continue;
339
340     anEdgeDir.SetCoord( anEdgeDir.X() / anEdgeLen,
341                         anEdgeDir.Y() / anEdgeLen );
342
343     gp_XY aFirstLinkDir( aFirstVertex.Coord() - aVertexCoord );
344     gp_XY aLastLinkDir ( aVertexCoord         - aLastVertex.Coord() );
345                       
346     Standard_Real aDist12 = aFirstLinkDir ^ anEdgeDir;
347     Standard_Real aDist23 = anEdgeDir     ^ aLastLinkDir;
348     if ( Abs( aDist12 ) < Precision || 
349          Abs( aDist23 ) < Precision )
350     {
351       continue;
352     }
353
354     Standard_Boolean isSensOK;
355     if ( myIsPositiveOrientation )
356       isSensOK = ( aDist12 > 0. && aDist23 > 0.);
357     else
358       isSensOK = ( aDist12 < 0. && aDist23 < 0.);
359     
360
361     BRepMesh_Edge aFirstLink( aNodes[1], aNodes[0], BRepMesh_Free );
362     BRepMesh_Edge aLastLink ( aNodes[2], aNodes[1], BRepMesh_Free );
363
364     Standard_Integer anEdgesInfo[3] = {
365       myMeshData->AddLink( aFirstLink ),
366       isPositive ? anEdgeId : -anEdgeId,
367       myMeshData->AddLink( aLastLink ) };
368
369     if ( isSensOK )
370     {
371       Standard_Integer anEdges[3];
372       Standard_Boolean anEdgesOri[3];
373       for ( Standard_Integer aTriLinkIt = 0; aTriLinkIt < 3; ++aTriLinkIt )
374       {
375         const Standard_Integer& anEdgeInfo = anEdgesInfo[aTriLinkIt];
376         anEdges[aTriLinkIt]    = Abs( anEdgeInfo );
377         anEdgesOri[aTriLinkIt] = anEdgeInfo > 0;
378       }
379
380       addTriangle( anEdges, anEdgesOri, aNodes );
381     }
382     else
383     {
384       if ( isPositive )
385         aLoopEdges.Append(  anEdges.Key() );
386       else
387         aLoopEdges.Append( -anEdges.Key() );
388         
389       if ( aFirstLinkDir.SquareModulus() > aLastLinkDir.SquareModulus() )
390         anExternalEdges.Append( Abs( anEdgesInfo[0] ) );
391       else
392         anExternalEdges.Append( Abs( anEdgesInfo[2] ) );
393     }
394   }
395   
396   thePoly.Clear();
397   while ( !anExternalEdges.IsEmpty() )
398   {
399     const BRepMesh_PairOfIndex& aPair = 
400       myMeshData->ElemConnectedTo( Abs( anExternalEdges.First() ) );
401     
402     
403     if ( !aPair.IsEmpty() )
404       deleteTriangle( aPair.FirstIndex(), thePoly );
405       
406     anExternalEdges.RemoveFirst();
407   }
408
409   for ( anEdges.Initialize( thePoly ); anEdges.More(); anEdges.Next() )
410   {
411     if ( myMeshData->ElemConnectedTo( anEdges.Key() ).IsEmpty() )
412       myMeshData->RemoveLink( anEdges.Key() );
413   }
414
415   while ( !aLoopEdges.IsEmpty() )
416   {
417     const BRepMesh_Edge& anEdge = GetEdge( Abs( aLoopEdges.First() ) );
418     if ( anEdge.Movability() != BRepMesh_Deleted )
419     {
420       Standard_Integer anEdgeIdx = aLoopEdges.First();
421       meshLeftPolygonOf( Abs( anEdgeIdx ), ( anEdgeIdx > 0 ) );
422     }
423       
424     aLoopEdges.RemoveFirst();
425   }
426 }
427
428 //=======================================================================
429 //function : createTrianglesOnNewVertices
430 //purpose  : Creation of triangles from the new nodes
431 //=======================================================================
432 void BRepMesh_Delaun::createTrianglesOnNewVertices( TColStd_Array1OfInteger& theVertexIndexes )
433 {
434   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
435
436   // Insertion of nodes :
437   Standard_Boolean isModify = Standard_True;
438   
439   Standard_Integer anIndex = theVertexIndexes.Lower();
440   Standard_Integer anUpper = theVertexIndexes.Upper();
441   for( ; anIndex <= anUpper; ++anIndex ) 
442   {
443     aLoopEdges.Clear();
444     
445     Standard_Integer aVertexIdx = theVertexIndexes( anIndex );    
446     const BRepMesh_Vertex& aVertex = GetVertex( aVertexIdx );
447
448     // Iterator in the list of indexes of circles containing the node
449     BRepMesh_ListOfInteger& aCirclesList = myCircles.Select( aVertex.Coord() );
450     
451     Standard_Integer onEgdeId = 0, aTriangleId = 0;
452     BRepMesh_ListOfInteger::Iterator aCircleIt( aCirclesList );
453     for ( ; aCircleIt.More(); aCircleIt.Next() )
454     {
455       // To add a node in the mesh it is necessary to check conditions: 
456       // - the node should be within the boundaries of the mesh and so in an existing triangle
457       // - all adjacent triangles should belong to a component connected with this triangle
458       if ( Contains( aCircleIt.Value(), aVertex, onEgdeId ) )
459       {
460         if ( onEgdeId == 0 )
461         {
462           aTriangleId = aCircleIt.Value();
463           aCirclesList.Remove( aCircleIt );
464           break;
465         }
466         else if ( GetEdge( onEgdeId ).Movability() == BRepMesh_Free )
467         {
468           aTriangleId = aCircleIt.Value();
469           aCirclesList.Remove( aCircleIt );
470           break;
471         }
472       }
473     }
474
475     if ( aTriangleId > 0 )
476     {
477       deleteTriangle( aTriangleId, aLoopEdges );
478
479       isModify = Standard_True;    
480       while ( isModify && !aCirclesList.IsEmpty() )
481       {
482         isModify = Standard_False;
483         BRepMesh_ListOfInteger::Iterator aCircleIt1( aCirclesList );
484         for ( ; aCircleIt1.More(); aCircleIt1.Next() )
485         {
486           Standard_Integer e[3];
487           Standard_Boolean o[3];
488           GetTriangle( aCircleIt1.Value() ).Edges( e[0], e[1], e[2],
489                                                    o[0], o[1], o[2] );
490                                                    
491           if ( aLoopEdges.IsBound( e[0] ) || 
492                aLoopEdges.IsBound( e[1] ) || 
493                aLoopEdges.IsBound( e[2] ) )
494           {
495             isModify = Standard_True;
496             deleteTriangle( aCircleIt1.Value(), aLoopEdges );
497             aCirclesList.Remove( aCircleIt1 );
498             break;
499           }
500         }
501       }
502
503       // Creation of triangles with the current node and free edges
504       // and removal of these edges from the list of free edges
505       createTriangles( aVertexIdx, aLoopEdges );
506     }
507   }
508   // Check that internal edges are not crossed by triangles
509   Handle(BRepMesh_MapOfInteger) anInternalEdges = InternalEdges();
510
511   // Destruction of triancles intersecting internal edges 
512   // and their replacement by makeshift triangles
513   BRepMesh_MapOfInteger::Iterator anInernalEdgesIt( *anInternalEdges );
514   for ( ; anInernalEdgesIt.More(); anInernalEdgesIt.Next() )
515   {
516     Standard_Integer aNbC;
517     aNbC = myMeshData->ElemConnectedTo( anInernalEdgesIt.Key() ).Extent();
518     if ( aNbC == 0 )
519     {
520       meshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_True  ); 
521       meshLeftPolygonOf( anInernalEdgesIt.Key(), Standard_False ); 
522     }
523   }
524
525   // Adjustment of meshes to boundary edges
526   frontierAdjust();
527 }
528
529 //=======================================================================
530 //function : isBoundToFrontier
531 //purpose  : Goes through the neighbour triangles around the given node 
532 //           started from the given link, returns TRUE if some triangle 
533 //           has a bounding frontier edge or FALSE elsewhere.
534 //           Stop link is used to prevent cycles.
535 //           Previous element Id is used to identify next neighbor element.
536 //=======================================================================
537 Standard_Boolean BRepMesh_Delaun::isBoundToFrontier(
538   const Standard_Integer theRefNodeId,
539   const Standard_Integer theRefLinkId,
540   const Standard_Integer theStopLinkId,
541   const Standard_Integer thePrevElementId)
542 {
543   const BRepMesh_PairOfIndex& aPair = 
544     myMeshData->ElemConnectedTo( theRefLinkId );
545   if ( aPair.IsEmpty() )
546     return Standard_False;
547
548   Standard_Integer aNbElements = aPair.Extent();
549   for ( Standard_Integer anElemIt = 1; anElemIt <= aNbElements; ++anElemIt )
550   {
551     const Standard_Integer aTriId = aPair.Index(anElemIt);
552     if ( aTriId < 0 || aTriId == thePrevElementId )
553       continue;
554
555     Standard_Integer anEdges[3];
556     Standard_Boolean anEdgesOri[3];
557     GetTriangle( aTriId ).Edges( anEdges[0], anEdges[1], anEdges[2],
558       anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
559
560     for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
561     {
562       const Standard_Integer anEdgeId = anEdges[anEdgeIt];
563       if ( anEdgeId == theRefLinkId )
564         continue;
565
566       if ( anEdgeId == theStopLinkId )
567         return Standard_False;
568
569       const BRepMesh_Edge& anEdge = GetEdge( anEdgeId );
570       if ( anEdge.FirstNode() != theRefNodeId &&
571            anEdge.LastNode()  != theRefNodeId )
572       {
573         continue;
574       }
575
576       if ( anEdge.Movability() != BRepMesh_Free )
577         return Standard_True;
578
579       return isBoundToFrontier( theRefNodeId, anEdgeId, theStopLinkId, aTriId );
580     }
581   }
582
583   return Standard_False;
584 }
585
586 //=======================================================================
587 //function : cleanupMesh
588 //purpose  : Cleanup mesh from the free triangles
589 //=======================================================================
590 void BRepMesh_Delaun::cleanupMesh()
591 {
592   while ( Standard_True )
593   {
594     BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
595     NCollection_Map<Standard_Integer> aDelTriangles;
596
597     Handle(BRepMesh_MapOfInteger) aFreeEdges = FreeEdges();
598     BRepMesh_MapOfInteger::Iterator aFreeEdgesIt( *aFreeEdges );
599     for ( ; aFreeEdgesIt.More(); aFreeEdgesIt.Next() )
600     {
601       const Standard_Integer& aFreeEdgeId = aFreeEdgesIt.Key();
602       const BRepMesh_Edge&    anEdge      = GetEdge( aFreeEdgeId );
603       if ( anEdge.Movability() == BRepMesh_Frontier )
604         continue;
605
606       const BRepMesh_PairOfIndex& aPair = 
607         myMeshData->ElemConnectedTo( aFreeEdgeId );
608       if ( aPair.IsEmpty() )
609       {
610         aLoopEdges.Bind( aFreeEdgeId, Standard_True );
611         continue;
612       }
613
614       Standard_Integer aTriId = aPair.FirstIndex();
615
616       // Check that the connected triangle is not surrounded by another triangles
617       Standard_Integer anEdges[3];
618       Standard_Boolean anEdgesOri[3];
619       GetTriangle( aTriId ).Edges( anEdges[0], anEdges[1], anEdges[2],
620         anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
621
622       Standard_Boolean isCanNotBeRemoved = Standard_True;
623       for ( Standard_Integer aCurEdgeIdx = 0; aCurEdgeIdx < 3; ++aCurEdgeIdx )
624       {
625         if ( anEdges[aCurEdgeIdx] != aFreeEdgeId )
626           continue;
627
628         for ( Standard_Integer anOtherEdgeIt = 1; anOtherEdgeIt <= 2; ++anOtherEdgeIt )
629         {
630           Standard_Integer anOtherEdgeId = ( aCurEdgeIdx + anOtherEdgeIt ) % 3;
631           const BRepMesh_PairOfIndex& anOtherEdgePair = 
632             myMeshData->ElemConnectedTo( anEdges[anOtherEdgeId] );
633
634           if ( anOtherEdgePair.Extent() < 2 )
635           {
636             isCanNotBeRemoved = Standard_False;
637             break;
638           }
639         }
640
641         break;
642       }
643
644       if ( isCanNotBeRemoved )
645         continue;
646
647       Standard_Boolean isConnected[2] = { Standard_False, Standard_False };
648       for ( Standard_Integer aLinkNodeIt = 0; aLinkNodeIt < 2; ++aLinkNodeIt )
649       {
650         isConnected[aLinkNodeIt] = isBoundToFrontier( ( aLinkNodeIt == 0 ) ? 
651           anEdge.FirstNode() : anEdge.LastNode(), 
652           aFreeEdgeId, aFreeEdgeId, -1);
653       }
654
655       if ( !isConnected[0] || !isConnected[1] )
656         aDelTriangles.Add( aTriId );
657     }
658
659     // Destruction of triangles :
660     Standard_Integer aDeletedTrianglesNb = 0;
661     NCollection_Map<Standard_Integer>::Iterator aDelTrianglesIt( aDelTriangles );
662     for ( ; aDelTrianglesIt.More(); aDelTrianglesIt.Next() )
663     {
664       deleteTriangle( aDelTrianglesIt.Key(), aLoopEdges );
665       aDeletedTrianglesNb++;
666     }
667
668     // Destruction of remaining hanging edges
669     BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
670     for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
671     {
672       if ( myMeshData->ElemConnectedTo( aLoopEdgesIt.Key() ).IsEmpty() )
673         myMeshData->RemoveLink( aLoopEdgesIt.Key() );
674     }
675
676     if ( aDeletedTrianglesNb == 0 )
677       break;
678   }
679 }
680
681 //=======================================================================
682 //function : frontierAdjust
683 //purpose  : Adjust the mesh on the frontier
684 //=======================================================================
685 void BRepMesh_Delaun::frontierAdjust()
686 {
687   Handle(BRepMesh_MapOfInteger)         aFrontier = Frontier();
688   NCollection_Vector<Standard_Integer>  aFailedFrontiers;
689   BRepMesh_MapOfIntegerInteger          aLoopEdges( 10, myMeshData->Allocator() );
690   BRepMesh_Delaun::HandleOfMapOfInteger aIntFrontierEdges = new NCollection_Map<Standard_Integer>();
691   for ( Standard_Integer aPass = 1; aPass <= 2; ++aPass )
692   {      
693     // 1 pass): find external triangles on boundary edges;
694     // 2 pass): find external triangles on boundary edges appeared 
695     //          during triangles replacement.
696     
697     BRepMesh_MapOfInteger::Iterator aFrontierIt( *aFrontier );
698     for ( ; aFrontierIt.More(); aFrontierIt.Next() )
699     {
700       Standard_Integer aFrontierId = aFrontierIt.Key();
701       const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( aFrontierId );
702       Standard_Integer aNbElem = aPair.Extent();
703       for( Standard_Integer aElemIt = 1; aElemIt <= aNbElem; ++aElemIt )
704       {
705         const Standard_Integer aPriorElemId = aPair.Index( aElemIt );
706         if( aPriorElemId < 0 )
707           continue;
708             
709         Standard_Integer e[3];
710         Standard_Boolean o[3];
711         GetTriangle( aPriorElemId ).Edges( e[0], e[1], e[2],
712                                            o[0], o[1], o[2] );
713
714         Standard_Boolean isTriangleFound = Standard_False;
715         for ( Standard_Integer n = 0; n < 3; ++n )
716         {
717           if ( aFrontierId == e[n] && !o[n] )
718           {
719             // Destruction  of external triangles on boundary edges
720             isTriangleFound = Standard_True;
721             deleteTriangle( aPriorElemId, aLoopEdges );
722             break;
723           }
724         }
725
726         if ( isTriangleFound )
727           break;
728       }
729     }
730
731     // destrucrion of remaining hanging edges :
732     BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
733     for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
734     {
735       Standard_Integer aLoopEdgeId = aLoopEdgesIt.Key();
736       if (myMeshData->ElemConnectedTo( aLoopEdgeId ).IsEmpty() )
737         myMeshData->RemoveLink( aLoopEdgeId );
738     }
739
740     // destruction of triangles crossing the boundary edges and 
741     // their replacement by makeshift triangles
742     for ( aFrontierIt.Reset(); aFrontierIt.More(); aFrontierIt.Next() )
743     {
744       Standard_Integer aFrontierId = aFrontierIt.Key();
745       if ( !myMeshData->ElemConnectedTo( aFrontierId ).IsEmpty() )
746         continue;
747
748       Standard_Boolean isSuccess = 
749         meshLeftPolygonOf( aFrontierId, Standard_True, aIntFrontierEdges );
750
751       if ( aPass == 2 && !isSuccess )
752         aFailedFrontiers.Append( aFrontierId );
753     }
754   }
755
756   cleanupMesh();
757
758   // When the mesh has been cleaned up, try to process frontier edges 
759   // once again to fill the possible gaps that might be occured in case of "saw" -
760   // situation when frontier edge has a triangle at a right side, but its free 
761   // links cross another frontieres  and meshLeftPolygonOf itself can't collect 
762   // a closed polygon.
763   NCollection_Vector<Standard_Integer>::Iterator aFailedFrontiersIt( aFailedFrontiers );
764   for ( ; aFailedFrontiersIt.More(); aFailedFrontiersIt.Next() )
765   {
766     Standard_Integer aFrontierId = aFailedFrontiersIt.Value();
767     if ( !myMeshData->ElemConnectedTo( aFrontierId ).IsEmpty() )
768       continue;
769
770     meshLeftPolygonOf( aFrontierId, Standard_True, aIntFrontierEdges );
771   }
772 }
773
774 //=======================================================================
775 //function : fillBndBox
776 //purpose  : Add boundig box for edge defined by start & end point to
777 //           the given vector of bounding boxes for triangulation edges
778 //=======================================================================
779 void BRepMesh_Delaun::fillBndBox( NCollection_Sequence<Bnd_B2d>& theBoxes,
780                                   const BRepMesh_Vertex&         theV1,
781                                   const BRepMesh_Vertex&         theV2 )
782 {
783   Bnd_B2d aBox;      
784   aBox.Add( theV1.Coord() );
785   aBox.Add( theV2.Coord() );
786   theBoxes.Append( aBox );
787 }
788
789 //=======================================================================
790 //function : meshLeftPolygonOf
791 //purpose  : Collect the polygon at the left of the given edge (material side)
792 //=======================================================================
793 Standard_Boolean BRepMesh_Delaun::meshLeftPolygonOf( 
794   const Standard_Integer                theStartEdgeId,
795   const Standard_Boolean                isForward,
796   BRepMesh_Delaun::HandleOfMapOfInteger theSkipped )
797 {
798   if ( !theSkipped.IsNull() && theSkipped->Contains( theStartEdgeId ) )
799     return Standard_True;
800
801   const BRepMesh_Edge& aRefEdge = GetEdge( theStartEdgeId );
802
803   TColStd_SequenceOfInteger aPolygon;
804   Standard_Integer aStartNode, aPivotNode;
805   if ( isForward )
806   {
807     aPolygon.Append( theStartEdgeId );
808     aStartNode = aRefEdge.FirstNode();
809     aPivotNode = aRefEdge.LastNode();
810   }
811   else
812   {
813     aPolygon.Append( -theStartEdgeId );
814     aStartNode = aRefEdge.LastNode();
815     aPivotNode = aRefEdge.FirstNode();
816   }
817
818
819   const BRepMesh_Vertex& aStartEdgeVertexS = GetVertex( aStartNode );
820   BRepMesh_Vertex        aPivotVertex      = GetVertex( aPivotNode );
821
822   gp_Vec2d aRefLinkDir( aPivotVertex.Coord() - 
823     aStartEdgeVertexS.Coord() );
824
825   if ( aRefLinkDir.SquareMagnitude() < Precision2 )
826     return Standard_True;
827
828   // Auxilary structures.
829   // Bounding boxes of polygon links to be used for preliminary
830   // analysis of intersections
831   NCollection_Sequence<Bnd_B2d> aBoxes;
832   fillBndBox( aBoxes, aStartEdgeVertexS, aPivotVertex );
833
834   // Hanging ends
835   NCollection_Map<Standard_Integer> aDeadLinks;
836
837   // Links are temporarily excluded from consideration
838   NCollection_Map<Standard_Integer> aLeprousLinks;
839   aLeprousLinks.Add( theStartEdgeId );
840
841   Standard_Boolean isSkipLeprous = Standard_True;
842   Standard_Integer aFirstNode    = aStartNode;
843   while ( aPivotNode != aFirstNode )
844   {
845     Bnd_B2d          aNextLinkBndBox;
846     gp_Vec2d         aNextLinkDir;
847     Standard_Integer aNextPivotNode = 0;
848
849     Standard_Integer aNextLinkId = findNextPolygonLink(
850       aFirstNode,
851       aPivotNode,     aPivotVertex,  aRefLinkDir, 
852       aBoxes,         aPolygon,      theSkipped,
853       isSkipLeprous,  aLeprousLinks, aDeadLinks, 
854       aNextPivotNode, aNextLinkDir,  aNextLinkBndBox );
855
856     if ( aNextLinkId != 0 )
857     {
858       aStartNode   = aPivotNode;
859       aRefLinkDir  = aNextLinkDir;
860
861       aPivotNode   = aNextPivotNode;
862       aPivotVertex = GetVertex( aNextPivotNode );
863
864       aBoxes.Append  ( aNextLinkBndBox );
865       aPolygon.Append( aNextLinkId );
866
867       isSkipLeprous = Standard_True;
868     }
869     else
870     {
871       // Nothing to do
872       if ( aPolygon.Length() == 1 )
873         return Standard_False;
874
875       // Return to the previous point
876       Standard_Integer aDeadLinkId = Abs( aPolygon.Last() );
877       aDeadLinks.Add      ( aDeadLinkId );
878
879       aLeprousLinks.Remove( aDeadLinkId );
880       aPolygon.Remove     ( aPolygon.Length() );
881       aBoxes.Remove       ( aBoxes.Length() );
882
883       Standard_Integer aPrevLinkInfo = aPolygon.Last();
884       const BRepMesh_Edge& aPrevLink = GetEdge( Abs( aPrevLinkInfo ) );
885
886       if( aPrevLinkInfo > 0 )
887       {
888         aStartNode = aPrevLink.FirstNode();
889         aPivotNode = aPrevLink.LastNode();
890       }
891       else
892       {
893         aStartNode = aPrevLink.LastNode();
894         aPivotNode = aPrevLink.FirstNode();
895       }
896       
897       aPivotVertex = GetVertex( aPivotNode );
898       aRefLinkDir = 
899         aPivotVertex.Coord() - GetVertex( aStartNode ).Coord();
900
901       isSkipLeprous = Standard_False;
902     }
903   }
904
905   if ( aPolygon.Length() < 3 )
906     return Standard_False;
907
908   cleanupPolygon( aPolygon, aBoxes );
909   meshPolygon   ( aPolygon, aBoxes, theSkipped );
910
911   return Standard_True;
912 }
913
914 //=======================================================================
915 //function : findNextPolygonLink
916 //purpose  : Find next link starting from the given node and has maximum
917 //           angle respect the given reference link.
918 //           Each time the next link is found other neighbor links at the 
919 //           pivot node are marked as leprous and will be excluded from 
920 //           consideration next time until a hanging end is occured.
921 //=======================================================================
922 Standard_Integer BRepMesh_Delaun::findNextPolygonLink(
923   const Standard_Integer&                     theFirstNode,
924   const Standard_Integer&                     thePivotNode,
925   const BRepMesh_Vertex&                      thePivotVertex,
926   const gp_Vec2d&                             theRefLinkDir,
927   const NCollection_Sequence<Bnd_B2d>&        theBoxes,
928   const TColStd_SequenceOfInteger&            thePolygon,
929   const BRepMesh_Delaun::HandleOfMapOfInteger theSkipped,
930   const Standard_Boolean&                     isSkipLeprous,
931   NCollection_Map<Standard_Integer>&          theLeprousLinks,
932   NCollection_Map<Standard_Integer>&          theDeadLinks,
933   Standard_Integer&                           theNextPivotNode,
934   gp_Vec2d&                                   theNextLinkDir,
935   Bnd_B2d&                                    theNextLinkBndBox )
936 {
937   // Find the next link having the greatest angle 
938   // respect to a direction of a reference one
939   Standard_Real aMaxAngle = myIsPositiveOrientation ?
940     RealFirst() : RealLast();
941
942   Standard_Integer aNextLinkId = 0;
943   BRepMesh_ListOfInteger::Iterator aLinkIt( myMeshData->LinkNeighboursOf( thePivotNode ) );
944   for ( ; aLinkIt.More(); aLinkIt.Next() )
945   {
946     const Standard_Integer& aNeighbourLinkInfo = aLinkIt.Value();
947     Standard_Integer        aNeighbourLinkId   = Abs( aNeighbourLinkInfo );
948
949     if ( theDeadLinks.Contains( aNeighbourLinkId ) ||
950        ( !theSkipped.IsNull() && theSkipped->Contains( aNeighbourLinkId ) ) )
951     {
952       continue;
953     }
954       
955     Standard_Boolean isLeprous = theLeprousLinks.Contains( aNeighbourLinkId );
956     if ( isSkipLeprous && isLeprous )
957       continue;
958
959     const BRepMesh_Edge& aNeighbourLink = GetEdge( aNeighbourLinkId );
960
961     // Determine whether the link belongs to the mesh
962     if ( aNeighbourLink.Movability() == BRepMesh_Free &&
963          myMeshData->ElemConnectedTo( aNeighbourLinkInfo ).IsEmpty() )
964     {
965       theDeadLinks.Add( aNeighbourLinkId );
966       continue;
967     }
968
969     Standard_Integer anOtherNode = aNeighbourLink.FirstNode();
970     if ( anOtherNode == thePivotNode )
971       anOtherNode = aNeighbourLink.LastNode();
972
973     gp_Vec2d aCurLinkDir( GetVertex( anOtherNode ).Coord() - 
974       thePivotVertex.Coord() );
975
976     if ( aCurLinkDir.SquareMagnitude() < Precision2 )
977     {
978       theDeadLinks.Add( aNeighbourLinkId );
979       continue;
980     }
981
982     if ( !isLeprous )
983       theLeprousLinks.Add( aNeighbourLinkId ); 
984
985     Standard_Real    anAngle    = theRefLinkDir.Angle( aCurLinkDir );
986     Standard_Boolean isFrontier = 
987       ( aNeighbourLink.Movability() == BRepMesh_Frontier );
988
989     Standard_Boolean isCheckPointOnEdge = Standard_True;
990     if ( isFrontier )
991     {
992       if ( Abs( Abs(anAngle) - M_PI ) < Precision::Angular() )
993       {
994         // Glued constrains - don't check intersection
995         isCheckPointOnEdge = Standard_False;
996         anAngle            = Abs( anAngle );
997       }
998     }
999
1000     if ( ( myIsPositiveOrientation && anAngle <= aMaxAngle ) ||
1001          (!myIsPositiveOrientation && anAngle >= aMaxAngle ) )
1002     {
1003       continue;
1004     }
1005
1006     Standard_Boolean isCheckEndPoints = ( anOtherNode != theFirstNode );
1007
1008     Bnd_B2d aBox;
1009     Standard_Boolean isNotIntersect =
1010       checkIntersection( aNeighbourLink, thePolygon, theBoxes, 
1011       isCheckEndPoints, isCheckPointOnEdge, Standard_True, aBox );
1012     
1013     if( isNotIntersect )
1014     {
1015       aMaxAngle         = anAngle;
1016
1017       theNextLinkDir    = aCurLinkDir;
1018       theNextPivotNode  = anOtherNode;
1019       theNextLinkBndBox = aBox;
1020
1021       aNextLinkId       = ( aNeighbourLink.FirstNode() == thePivotNode ) ?
1022         aNeighbourLinkId : -aNeighbourLinkId;
1023     }
1024   }
1025
1026   return aNextLinkId;
1027 }
1028
1029 //=======================================================================
1030 //function : checkIntersection
1031 //purpose  : Check is the given link intersects the polygon boundaries.
1032 //           Returns bounding box for the given link trough the 
1033 //           <theLinkBndBox> parameter.
1034 //=======================================================================
1035 Standard_Boolean BRepMesh_Delaun::checkIntersection( 
1036   const BRepMesh_Edge&                 theLink,
1037   const TColStd_SequenceOfInteger&     thePolygon,
1038   const NCollection_Sequence<Bnd_B2d>& thePolyBoxes,
1039   const Standard_Boolean               isConsiderEndPointTouch,
1040   const Standard_Boolean               isConsiderPointOnEdge,
1041   const Standard_Boolean               isSkipLastEdge,
1042   Bnd_B2d&                             theLinkBndBox ) const
1043 {
1044   theLinkBndBox.Add( GetVertex( theLink.FirstNode() ).Coord() );
1045   theLinkBndBox.Add( GetVertex( theLink.LastNode()  ).Coord() );
1046
1047   Standard_Integer aPolyLen = thePolygon.Length();
1048   // Don't check intersection with the last link
1049   if ( isSkipLastEdge )
1050     --aPolyLen;
1051
1052   Standard_Boolean isFrontier = 
1053     ( theLink.Movability() == BRepMesh_Frontier );
1054
1055   for ( Standard_Integer aPolyIt = 1; aPolyIt <= aPolyLen; ++aPolyIt )
1056   {
1057     if ( !theLinkBndBox.IsOut( thePolyBoxes.Value( aPolyIt ) ) )
1058     {
1059       // intersection is possible...
1060       Standard_Integer aPolyLinkId   = Abs( thePolygon( aPolyIt ) );
1061       const BRepMesh_Edge& aPolyLink = GetEdge( aPolyLinkId );
1062
1063       // skip intersections between frontier edges
1064       if ( aPolyLink.Movability() == BRepMesh_Frontier && isFrontier )
1065         continue;
1066
1067       gp_Pnt2d anIntPnt;
1068       IntFlag aIntFlag = intSegSeg( theLink, aPolyLink, 
1069         isConsiderEndPointTouch, isConsiderPointOnEdge, anIntPnt );
1070
1071       if ( aIntFlag != BRepMesh_Delaun::NoIntersection )
1072         return Standard_False;
1073     }
1074   }
1075
1076   // Ok, no intersection
1077   return Standard_True;
1078 }
1079
1080 //=======================================================================
1081 //function : addTriangle
1082 //purpose  : Add a triangle based on the given oriented edges into mesh
1083 //=======================================================================
1084 inline void BRepMesh_Delaun::addTriangle( const Standard_Integer (&theEdgesId)[3],
1085                                           const Standard_Boolean (&theEdgesOri)[3],
1086                                           const Standard_Integer (&theNodesId)[3] )
1087 {
1088   Standard_Integer aNewTriangleId = 
1089     myMeshData->AddElement( BRepMesh_Triangle( 
1090       theEdgesId[0],  theEdgesId[1],  theEdgesId[2], 
1091       theEdgesOri[0], theEdgesOri[1], theEdgesOri[2],
1092       BRepMesh_Free ) );
1093
1094   Standard_Boolean isAdded = myCircles.Add( 
1095     GetVertex( theNodesId[0] ).Coord(), 
1096     GetVertex( theNodesId[1] ).Coord(),
1097     GetVertex( theNodesId[2] ).Coord(),
1098     aNewTriangleId );
1099     
1100   if ( !isAdded )
1101     myMeshData->RemoveElement( aNewTriangleId );
1102 }
1103
1104 //=======================================================================
1105 //function : cleanupPolygon
1106 //purpose  : Remove internal triangles from the given polygon
1107 //=======================================================================
1108 void BRepMesh_Delaun::cleanupPolygon( const TColStd_SequenceOfInteger&     thePolygon,
1109                                       const NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
1110 {
1111   Standard_Integer aPolyLen = thePolygon.Length();
1112   if ( aPolyLen < 3 )
1113     return;
1114
1115   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
1116   NCollection_Map<Standard_Integer>    anIgnoredEdges;
1117   NCollection_Map<Standard_Integer>    aPolyVerticesFindMap;
1118   NCollection_Vector<Standard_Integer> aPolyVertices;
1119   // Collect boundary vertices of the polygon
1120   for ( Standard_Integer aPolyIt = 1; aPolyIt <= aPolyLen; ++aPolyIt )
1121   {
1122     Standard_Integer aPolyEdgeInfo = thePolygon( aPolyIt );
1123     Standard_Integer aPolyEdgeId   = Abs( aPolyEdgeInfo );
1124     anIgnoredEdges.Add( aPolyEdgeId );
1125
1126     Standard_Boolean isForward = ( aPolyEdgeInfo > 0 );
1127     const BRepMesh_PairOfIndex& aPair = 
1128       myMeshData->ElemConnectedTo( aPolyEdgeId );
1129
1130     Standard_Integer anElemIt = 1;
1131     for ( ; anElemIt <= aPair.Extent(); ++anElemIt )
1132     {
1133       Standard_Integer anElemId = aPair.Index( anElemIt );
1134       if ( anElemId < 0 )
1135         continue;
1136
1137       Standard_Integer anEdges[3];
1138       Standard_Boolean anEdgesOri[3];
1139       GetTriangle( anElemId ).Edges( anEdges[0], anEdges[1], anEdges[2],
1140         anEdgesOri[0], anEdgesOri[1], anEdgesOri[2] );
1141
1142       Standard_Integer isTriangleFound = Standard_False;
1143       for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
1144       {
1145         if ( anEdges[anEdgeIt]    == aPolyEdgeId && 
1146              anEdgesOri[anEdgeIt] == isForward )
1147         {
1148           isTriangleFound = Standard_True;
1149           deleteTriangle( anElemId, aLoopEdges );
1150           break;
1151         }
1152       }
1153
1154       if ( isTriangleFound )
1155         break;
1156     }
1157
1158     // Skip a neighbor link to extract unique vertices each time
1159     if ( aPolyIt % 2 )
1160     {
1161       const BRepMesh_Edge& aPolyEdge = GetEdge( aPolyEdgeId );
1162       Standard_Integer aFirstVertex  = aPolyEdge.FirstNode();
1163       Standard_Integer aLastVertex   = aPolyEdge.LastNode();
1164
1165       aPolyVerticesFindMap.Add( aFirstVertex );
1166       aPolyVerticesFindMap.Add( aLastVertex );
1167
1168       if ( aPolyEdgeInfo > 0 )
1169       {
1170         aPolyVertices.Append( aFirstVertex );
1171         aPolyVertices.Append( aLastVertex );
1172       }
1173       else
1174       {
1175         aPolyVertices.Append( aLastVertex );
1176         aPolyVertices.Append( aFirstVertex );
1177       }
1178     }
1179   }
1180
1181   // Make closed sequence
1182   if ( aPolyVertices.First() != aPolyVertices.Last() )
1183     aPolyVertices.Append( aPolyVertices.First() );
1184
1185   NCollection_Map<Standard_Integer> aSurvivedLinks( anIgnoredEdges );
1186
1187   Standard_Integer aPolyVertIt          = 0;
1188   Standard_Integer anUniqueVerticesNum  = aPolyVertices.Length() - 1;
1189   for ( ; aPolyVertIt < anUniqueVerticesNum; ++aPolyVertIt )
1190   {
1191     killTrianglesAroundVertex( aPolyVertices( aPolyVertIt ),
1192       aPolyVertices, aPolyVerticesFindMap, thePolygon,
1193       thePolyBoxes, aSurvivedLinks, aLoopEdges );
1194   }
1195
1196   BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
1197   for ( ; aLoopEdgesIt.More(); aLoopEdgesIt.Next() )
1198   {
1199     const Standard_Integer& aLoopEdgeId = aLoopEdgesIt.Key();
1200     if ( anIgnoredEdges.Contains( aLoopEdgeId ) )
1201       continue;
1202
1203     if ( myMeshData->ElemConnectedTo( aLoopEdgeId ).IsEmpty() )
1204       myMeshData->RemoveLink( aLoopEdgesIt.Key() );
1205   }
1206 }
1207
1208 //=======================================================================
1209 //function : killTrianglesAroundVertex
1210 //purpose  : Remove all triangles and edges that are placed
1211 //           inside the polygon or crossed it.
1212 //=======================================================================
1213 void BRepMesh_Delaun::killTrianglesAroundVertex( 
1214   const Standard_Integer                        theZombieNodeId,
1215   const NCollection_Vector<Standard_Integer>&   thePolyVertices,
1216   const NCollection_Map<Standard_Integer>&      thePolyVerticesFindMap,
1217   const TColStd_SequenceOfInteger&              thePolygon,
1218   const NCollection_Sequence<Bnd_B2d>&          thePolyBoxes,
1219   NCollection_Map<Standard_Integer>&            theSurvivedLinks,
1220   BRepMesh_MapOfIntegerInteger&                 theLoopEdges )
1221 {
1222   BRepMesh_ListOfInteger::Iterator aNeighborsIt = 
1223       myMeshData->LinkNeighboursOf( theZombieNodeId );
1224
1225   // Try to infect neighbor nodes
1226   NCollection_Vector<Standard_Integer> aVictimNodes;
1227   for ( ; aNeighborsIt.More(); aNeighborsIt.Next() )
1228   {
1229     const Standard_Integer& aNeighborLinkId = aNeighborsIt.Value();
1230     if ( theSurvivedLinks.Contains( aNeighborLinkId ) )
1231       continue;
1232
1233     const BRepMesh_Edge& aNeighborLink = GetEdge( aNeighborLinkId );
1234     if ( aNeighborLink.Movability() == BRepMesh_Frontier )
1235     {
1236       // Though, if it lies onto the polygon boundary -
1237       // take its triangles
1238       Bnd_B2d aBox;
1239       Standard_Boolean isNotIntersect =
1240         checkIntersection( aNeighborLink, thePolygon, 
1241         thePolyBoxes, Standard_False, Standard_True, 
1242         Standard_False, aBox );
1243
1244       if ( isNotIntersect )
1245       {
1246         // Don't touch it
1247         theSurvivedLinks.Add( aNeighborLinkId );
1248         continue;
1249       }
1250     }
1251     else
1252     {
1253       Standard_Integer anOtherNode = aNeighborLink.FirstNode();
1254       if ( anOtherNode == theZombieNodeId )
1255         anOtherNode = aNeighborLink.LastNode();
1256
1257       // Possible sacrifice
1258       if ( !thePolyVerticesFindMap.Contains( anOtherNode ) )
1259       {
1260         if ( isVertexInsidePolygon( anOtherNode, thePolyVertices ) )
1261         {
1262           // Got you!
1263           aVictimNodes.Append( anOtherNode );
1264         }
1265         else
1266         {
1267           // Lucky. But it may intersect the polygon boundary.
1268           // Let's check it!
1269           killTrianglesOnIntersectingLinks( aNeighborLinkId, 
1270             aNeighborLink, anOtherNode, thePolygon, 
1271             thePolyBoxes, theSurvivedLinks, theLoopEdges );
1272
1273           continue;
1274         }
1275       }
1276     }
1277
1278     // Add link to the survivers to avoid cycling
1279     theSurvivedLinks.Add( aNeighborLinkId );
1280     killLinkTriangles( aNeighborLinkId, theLoopEdges );
1281   }
1282
1283   // Go and do your job!
1284   NCollection_Vector<Standard_Integer>::Iterator aVictimIt( aVictimNodes );
1285   for ( ; aVictimIt.More(); aVictimIt.Next() )
1286   {
1287     killTrianglesAroundVertex( aVictimIt.Value(), thePolyVertices,
1288       thePolyVerticesFindMap, thePolygon, thePolyBoxes,
1289       theSurvivedLinks, theLoopEdges );
1290   }
1291 }
1292
1293 //=======================================================================
1294 //function : isVertexInsidePolygon
1295 //purpose  : Checks is the given vertex lies inside the polygon
1296 //=======================================================================
1297 Standard_Boolean BRepMesh_Delaun::isVertexInsidePolygon( 
1298   const Standard_Integer&                     theVertexId,
1299   const NCollection_Vector<Standard_Integer>& thePolygonVertices ) const
1300 {
1301   Standard_Integer aPolyLen = thePolygonVertices.Length();
1302   if ( aPolyLen < 3 )
1303     return Standard_False;
1304
1305
1306   const gp_XY aCenterPointXY = GetVertex( theVertexId ).Coord();
1307
1308   const BRepMesh_Vertex& aFirstVertex = GetVertex( thePolygonVertices( 0 ) );
1309   gp_Vec2d aPrevVertexDir( aFirstVertex.Coord() - aCenterPointXY );
1310   if ( aPrevVertexDir.SquareMagnitude() < Precision2 )
1311     return Standard_True;
1312
1313   Standard_Real aTotalAng = 0.0;
1314   for ( Standard_Integer aPolyIt = 1; aPolyIt < aPolyLen; ++aPolyIt )
1315   {
1316     const BRepMesh_Vertex& aPolyVertex = GetVertex( thePolygonVertices( aPolyIt ) );
1317
1318     gp_Vec2d aCurVertexDir( aPolyVertex.Coord() - aCenterPointXY );
1319     if ( aCurVertexDir.SquareMagnitude() < Precision2 )
1320       return Standard_True;
1321
1322     aTotalAng     += aCurVertexDir.Angle( aPrevVertexDir );
1323     aPrevVertexDir = aCurVertexDir;
1324   }
1325   
1326   if ( Abs( Angle2PI - aTotalAng ) > Precision::Angular() )
1327     return Standard_False;
1328
1329   return Standard_True;
1330 }
1331
1332 //=======================================================================
1333 //function : killTrianglesOnIntersectingLinks
1334 //purpose  : Checks is the given link crosses the polygon boundary.
1335 //           If yes, kills its triangles and checks neighbor links on
1336 //           boundary intersection. Does nothing elsewhere.
1337 //=======================================================================
1338 void BRepMesh_Delaun::killTrianglesOnIntersectingLinks( 
1339   const Standard_Integer&               theLinkToCheckId, 
1340   const BRepMesh_Edge&                  theLinkToCheck,
1341   const Standard_Integer&               theEndPoint,
1342   const TColStd_SequenceOfInteger&      thePolygon,
1343   const NCollection_Sequence<Bnd_B2d>&  thePolyBoxes,
1344   NCollection_Map<Standard_Integer>&    theSurvivedLinks,
1345   BRepMesh_MapOfIntegerInteger&         theLoopEdges )
1346 {
1347   if ( theSurvivedLinks.Contains( theLinkToCheckId ) )
1348     return;
1349
1350   Bnd_B2d aBox;
1351   Standard_Boolean isNotIntersect =
1352     checkIntersection( theLinkToCheck, thePolygon, 
1353       thePolyBoxes, Standard_False, Standard_False, 
1354       Standard_False, aBox );
1355
1356   theSurvivedLinks.Add( theLinkToCheckId );
1357
1358   if ( isNotIntersect )
1359     return;
1360
1361   killLinkTriangles( theLinkToCheckId, theLoopEdges );
1362
1363   BRepMesh_ListOfInteger::Iterator aNeighborsIt = 
1364     myMeshData->LinkNeighboursOf( theEndPoint );
1365
1366   for ( ; aNeighborsIt.More(); aNeighborsIt.Next() )
1367   {
1368     const Standard_Integer& aNeighborLinkId = aNeighborsIt.Value();
1369     const BRepMesh_Edge&    aNeighborLink   = GetEdge( aNeighborLinkId );
1370     Standard_Integer anOtherNode = aNeighborLink.FirstNode();
1371     if ( anOtherNode == theEndPoint )
1372       anOtherNode = aNeighborLink.LastNode();
1373
1374     killTrianglesOnIntersectingLinks( aNeighborLinkId, 
1375       aNeighborLink, anOtherNode, thePolygon, 
1376       thePolyBoxes, theSurvivedLinks, theLoopEdges );
1377   }
1378 }
1379
1380 //=======================================================================
1381 //function : killLinkTriangles
1382 //purpose  : Kill triangles bound to the given link.
1383 //=======================================================================
1384 void BRepMesh_Delaun::killLinkTriangles( 
1385   const Standard_Integer&       theLinkId, 
1386   BRepMesh_MapOfIntegerInteger& theLoopEdges )
1387 {
1388   const BRepMesh_PairOfIndex& aPair = 
1389       myMeshData->ElemConnectedTo( theLinkId );
1390
1391   Standard_Integer anElemNb = aPair.Extent();
1392   for ( Standard_Integer aPairIt = 1; aPairIt <= anElemNb; ++aPairIt )
1393   {
1394     Standard_Integer anElemId = aPair.FirstIndex();
1395     if ( anElemId < 0 )
1396       continue;
1397
1398     deleteTriangle( anElemId, theLoopEdges );
1399   }
1400 }
1401
1402 //=======================================================================
1403 //function : getOrientedNodes
1404 //purpose  : Returns start and end nodes of the given edge in respect to 
1405 //           its orientation.
1406 //=======================================================================
1407 void BRepMesh_Delaun::getOrientedNodes(const BRepMesh_Edge&   theEdge,
1408                                        const Standard_Boolean isForward,
1409                                        Standard_Integer       *theNodes) const
1410 {
1411   if ( isForward )
1412   {
1413     theNodes[0] = theEdge.FirstNode();
1414     theNodes[1] = theEdge.LastNode();
1415   }
1416   else
1417   {
1418     theNodes[0] = theEdge.LastNode();
1419     theNodes[1] = theEdge.FirstNode();
1420   }
1421 }
1422
1423 //=======================================================================
1424 //function : processLoop
1425 //purpose  : Processes loop within the given polygon formed by range of 
1426 //           its links specified by start and end link indices.
1427 //=======================================================================
1428 void BRepMesh_Delaun::processLoop(const Standard_Integer               theLinkFrom,
1429                                   const Standard_Integer               theLinkTo,
1430                                   const TColStd_SequenceOfInteger&     thePolygon,
1431                                   const NCollection_Sequence<Bnd_B2d>& thePolyBoxes)
1432 {
1433   Standard_Integer aNbOfLinksInLoop = theLinkTo - theLinkFrom - 1;
1434   if ( aNbOfLinksInLoop < 3 )
1435     return;
1436
1437   TColStd_SequenceOfInteger     aPolygon;
1438   NCollection_Sequence<Bnd_B2d> aPolyBoxes;
1439   for ( ; aNbOfLinksInLoop > 0; --aNbOfLinksInLoop )
1440   {
1441     Standard_Integer aLoopLinkIndex = theLinkFrom + aNbOfLinksInLoop;
1442     aPolygon  .Prepend( thePolygon  ( aLoopLinkIndex ) );
1443     aPolyBoxes.Prepend( thePolyBoxes( aLoopLinkIndex ) );
1444   }
1445   meshPolygon( aPolygon, aPolyBoxes );
1446 }
1447
1448 //=======================================================================
1449 //function : createAndReplacePolygonLink
1450 //purpose  : Creates new link based on the given nodes and updates the 
1451 //           given polygon.
1452 //=======================================================================
1453 Standard_Integer BRepMesh_Delaun::createAndReplacePolygonLink(
1454   const Standard_Integer         *theNodes,
1455   const gp_Pnt2d                 *thePnts,
1456   const Standard_Integer         theRootIndex,
1457   const ReplaceFlag              theReplaceFlag,
1458   TColStd_SequenceOfInteger&     thePolygon,
1459   NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
1460 {
1461   Standard_Integer aNewEdgeId = 
1462     myMeshData->AddLink( BRepMesh_Edge(
1463     theNodes[0], theNodes[1], BRepMesh_Free ) );
1464
1465   Bnd_B2d aNewBox;
1466   aNewBox.Add( thePnts[0] );
1467   aNewBox.Add( thePnts[1] );
1468
1469   switch ( theReplaceFlag )
1470   {
1471   case BRepMesh_Delaun::Replace:
1472     thePolygon  .SetValue( theRootIndex, aNewEdgeId );
1473     thePolyBoxes.SetValue( theRootIndex, aNewBox );
1474     break;
1475
1476   case BRepMesh_Delaun::InsertAfter:
1477     thePolygon  .InsertAfter( theRootIndex, aNewEdgeId );
1478     thePolyBoxes.InsertAfter( theRootIndex, aNewBox );
1479     break;
1480
1481   case BRepMesh_Delaun::InsertBefore:
1482     thePolygon  .InsertBefore( theRootIndex, aNewEdgeId );
1483     thePolyBoxes.InsertBefore( theRootIndex, aNewBox );
1484     break;
1485   }
1486
1487   return aNewEdgeId;
1488 }
1489
1490 //=======================================================================
1491 //function : meshPolygon
1492 //purpose  : 
1493 //=======================================================================
1494 void BRepMesh_Delaun::meshPolygon( TColStd_SequenceOfInteger&            thePolygon,
1495                                    NCollection_Sequence<Bnd_B2d>&        thePolyBoxes,
1496                                    BRepMesh_Delaun::HandleOfMapOfInteger theSkipped )
1497 {
1498   // Check is the source polygon elementary
1499   if ( meshElementaryPolygon( thePolygon ) )
1500     return;
1501
1502   // Check and correct boundary edges
1503   Standard_Integer aPolyLen = thePolygon.Length();
1504   const Standard_Real aPolyArea      = Abs( polyArea( thePolygon, 1, aPolyLen ) );
1505   const Standard_Real aSmallLoopArea = 0.001 * aPolyArea;
1506   for ( Standard_Integer aPolyIt = 1; aPolyIt < aPolyLen; ++aPolyIt )
1507   {
1508     Standard_Integer aCurEdgeInfo = thePolygon( aPolyIt );
1509     Standard_Integer aCurEdgeId   = Abs( aCurEdgeInfo );
1510     const BRepMesh_Edge* aCurEdge = &GetEdge( aCurEdgeId );
1511     if ( aCurEdge->Movability() != BRepMesh_Frontier )
1512       continue;
1513
1514     Standard_Integer aCurNodes[2];
1515     getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aCurNodes );
1516
1517     gp_Pnt2d aCurPnts[2] = { 
1518       GetVertex(aCurNodes[0]).Coord(),
1519       GetVertex(aCurNodes[1]).Coord()
1520     };
1521
1522     gp_Vec2d aCurVec( aCurPnts[0], aCurPnts[1] );
1523
1524     // check further links
1525     Standard_Integer aNextPolyIt = aPolyIt + 1;
1526     for ( ; aNextPolyIt <= aPolyLen; ++aNextPolyIt )
1527     {
1528       Standard_Integer aNextEdgeInfo = thePolygon( aNextPolyIt );
1529       Standard_Integer aNextEdgeId   = Abs( aNextEdgeInfo );
1530       const BRepMesh_Edge* aNextEdge = &GetEdge( aNextEdgeId );
1531       if ( aNextEdge->Movability() != BRepMesh_Frontier )
1532         continue;
1533
1534       Standard_Integer aNextNodes[2];
1535       getOrientedNodes( *aNextEdge, aNextEdgeInfo > 0, aNextNodes );
1536
1537       gp_Pnt2d aNextPnts[2] = { 
1538         GetVertex(aNextNodes[0]).Coord(),
1539         GetVertex(aNextNodes[1]).Coord() 
1540       };
1541
1542       gp_Pnt2d anIntPnt;
1543       IntFlag aIntFlag = intSegSeg( *aCurEdge, *aNextEdge, 
1544         Standard_False, Standard_True, anIntPnt );
1545
1546       if ( aIntFlag == BRepMesh_Delaun::NoIntersection )
1547         continue;
1548
1549       Standard_Boolean isRemoveFromFirst  = Standard_False;
1550       Standard_Boolean isAddReplacingEdge = Standard_True;
1551       Standard_Integer aIndexToRemoveTo   = aNextPolyIt;
1552       if ( aIntFlag == BRepMesh_Delaun::Cross )
1553       {
1554         Standard_Real aLoopArea = polyArea( thePolygon, aPolyIt + 1, aNextPolyIt );
1555         gp_Vec2d aVec1( anIntPnt, aCurPnts [1] );
1556         gp_Vec2d aVec2( anIntPnt, aNextPnts[0] );
1557
1558         aLoopArea += ( aVec1 ^ aVec2 ) / 2.;
1559         if ( Abs( aLoopArea ) > aSmallLoopArea )
1560         {
1561           aNextNodes[1] = aCurNodes[0];
1562           aNextPnts [1] = aCurPnts [0];
1563
1564           aNextEdgeId = Abs( createAndReplacePolygonLink( aNextNodes, aNextPnts,
1565             aNextPolyIt, BRepMesh_Delaun::Replace, thePolygon, thePolyBoxes ) );
1566
1567           processLoop( aPolyIt, aNextPolyIt, thePolygon, thePolyBoxes );
1568           return;
1569         }
1570
1571         Standard_Real aDist1 = anIntPnt.SquareDistance(aNextPnts[0]);
1572         Standard_Real aDist2 = anIntPnt.SquareDistance(aNextPnts[1]);
1573
1574         // Choose node with lower distance
1575         const Standard_Boolean isCloseToStart = ( aDist1 < aDist2 );
1576         const Standard_Integer aEndPointIndex = isCloseToStart ? 0 : 1;
1577         aCurNodes[1] = aNextNodes[aEndPointIndex];
1578         aCurPnts [1] = aNextPnts [aEndPointIndex];
1579
1580         if ( isCloseToStart )
1581           --aIndexToRemoveTo;
1582
1583         // In this context only intersections between frontier edges
1584         // are possible. If intersection between edges of different
1585         // types occured - treat this case as invalid (i.e. result 
1586         // might not reflect the expectations).
1587         if ( !theSkipped.IsNull() )
1588         {
1589           Standard_Integer aSkippedLinkIt = aPolyIt;
1590           for ( ; aSkippedLinkIt <= aIndexToRemoveTo; ++aSkippedLinkIt )
1591             theSkipped->Add( Abs( thePolygon( aSkippedLinkIt ) ) );
1592         }
1593       }
1594       else if ( aIntFlag == BRepMesh_Delaun::PointOnEdge )
1595       {
1596         // Indentify chopping link 
1597         Standard_Boolean isFirstChopping = Standard_False;
1598         Standard_Integer aCheckPointIt = 0;
1599         for ( ; aCheckPointIt < 2; ++aCheckPointIt )
1600         {
1601           gp_Pnt2d& aRefPoint = aCurPnts[aCheckPointIt];
1602           // Check is second link touches the first one
1603           gp_Vec2d aVec1( aRefPoint, aNextPnts[0] );
1604           gp_Vec2d aVec2( aRefPoint, aNextPnts[1] );
1605           if ( Abs( aVec1 ^ aVec2 ) < Precision )
1606           {
1607             isFirstChopping = Standard_True;
1608             break;
1609           }
1610         }
1611  
1612         if ( isFirstChopping )
1613         {
1614           // Split second link
1615           isAddReplacingEdge = Standard_False;
1616           isRemoveFromFirst  = ( aCheckPointIt == 0 );
1617
1618           Standard_Integer aSplitLink[3] = {
1619             aNextNodes[0],
1620             aCurNodes [aCheckPointIt],
1621             aNextNodes[1]
1622           };
1623
1624           gp_Pnt2d aSplitPnts[3] = {
1625             aNextPnts[0],
1626             aCurPnts [aCheckPointIt],
1627             aNextPnts[1]
1628           };
1629
1630           Standard_Integer aSplitLinkIt = 0;
1631           for ( ; aSplitLinkIt < 2; ++aSplitLinkIt )
1632           {
1633             createAndReplacePolygonLink( &aSplitLink[aSplitLinkIt],
1634               &aSplitPnts[aSplitLinkIt], aNextPolyIt, ( aSplitLinkIt == 0 ) ? 
1635               BRepMesh_Delaun::Replace : BRepMesh_Delaun::InsertAfter,
1636               thePolygon, thePolyBoxes );
1637           }
1638
1639           processLoop( aPolyIt + aCheckPointIt, aIndexToRemoveTo,
1640             thePolygon, thePolyBoxes );
1641         }
1642         else
1643         {
1644           // Split first link
1645           Standard_Integer aSplitLinkNodes[2] = {
1646             aNextNodes[1],
1647             aCurNodes [1]
1648           };
1649
1650           gp_Pnt2d aSplitLinkPnts[2] = {
1651             aNextPnts[1],
1652             aCurPnts [1]
1653           };
1654           createAndReplacePolygonLink( aSplitLinkNodes, aSplitLinkPnts,
1655             aPolyIt, BRepMesh_Delaun::InsertAfter, thePolygon, thePolyBoxes );
1656
1657           aCurNodes[1] = aNextNodes[1];
1658           aCurPnts [1] = aNextPnts [1];
1659           ++aIndexToRemoveTo;
1660
1661           processLoop( aPolyIt + 1, aIndexToRemoveTo, 
1662             thePolygon, thePolyBoxes );
1663         }
1664       }
1665       else if ( aIntFlag == BRepMesh_Delaun::Glued )
1666       {
1667         if ( aCurNodes[1] == aNextNodes[0] )
1668         {
1669           aCurNodes[1] = aNextNodes[1];
1670           aCurPnts [1] = aNextPnts [1];
1671         }
1672         // TODO: Non-adjacent glued links within the polygon
1673       }
1674       else if ( aIntFlag == BRepMesh_Delaun::Same )
1675       {
1676         processLoop( aPolyIt, aNextPolyIt, thePolygon, thePolyBoxes );
1677
1678         isRemoveFromFirst  = Standard_True;
1679         isAddReplacingEdge = Standard_False;
1680       }
1681       else
1682         continue; // Not supported type
1683
1684       if ( isAddReplacingEdge )
1685       {
1686         aCurEdgeId = Abs( createAndReplacePolygonLink( aCurNodes, aCurPnts,
1687           aPolyIt, BRepMesh_Delaun::Replace, thePolygon, thePolyBoxes ) );
1688
1689         aCurEdge = &GetEdge( aCurEdgeId );
1690         aCurVec  = gp_Vec2d( aCurPnts[0], aCurPnts[1] );
1691       }
1692
1693       Standard_Integer aIndexToRemoveFrom = 
1694         isRemoveFromFirst ? aPolyIt : aPolyIt + 1;
1695
1696       thePolygon  .Remove( aIndexToRemoveFrom, aIndexToRemoveTo );
1697       thePolyBoxes.Remove( aIndexToRemoveFrom, aIndexToRemoveTo );
1698
1699       aPolyLen = thePolygon.Length();
1700       if ( isRemoveFromFirst )
1701       {
1702         --aPolyIt;
1703         break;
1704       }
1705
1706       aNextPolyIt = aPolyIt;
1707     }
1708   }
1709
1710   meshSimplePolygon( thePolygon, thePolyBoxes );
1711 }
1712
1713 //=======================================================================
1714 //function : meshElementaryPolygon
1715 //purpose  : Triangulation of closed polygon containing only three edges.
1716 //=======================================================================
1717 inline Standard_Boolean BRepMesh_Delaun::meshElementaryPolygon( 
1718   const TColStd_SequenceOfInteger& thePolygon )
1719 {
1720   Standard_Integer aPolyLen = thePolygon.Length();
1721   if ( aPolyLen < 3 )
1722     return Standard_True;
1723   else if ( aPolyLen > 3 )
1724     return Standard_False;
1725
1726   // Just create a triangle
1727   Standard_Integer anEdges[3];
1728   Standard_Boolean anEdgesOri[3];
1729
1730   for ( Standard_Integer anEdgeIt = 0; anEdgeIt < 3; ++anEdgeIt )
1731   {
1732     Standard_Integer anEdgeInfo = thePolygon( anEdgeIt + 1 );
1733     anEdges[anEdgeIt]           = Abs( anEdgeInfo );
1734     anEdgesOri[anEdgeIt]        = ( anEdgeInfo > 0 );
1735   }
1736     
1737   const BRepMesh_Edge& anEdge1 = GetEdge( anEdges[0] );
1738   const BRepMesh_Edge& anEdge2 = GetEdge( anEdges[1] );
1739   
1740   Standard_Integer aNodes[3] = { anEdge1.FirstNode(),
1741                                  anEdge1.LastNode(),
1742                                  anEdge2.FirstNode() };
1743   if ( aNodes[2] == aNodes[0] || 
1744        aNodes[2] == aNodes[1] )
1745   {
1746     aNodes[2] = anEdge2.LastNode();
1747   }
1748
1749   addTriangle( anEdges, anEdgesOri, aNodes );
1750   return Standard_True;
1751 }
1752
1753 //=======================================================================
1754 //function : meshSimplePolygon
1755 //purpose  : Triangulatiion of a closed simple polygon (polygon without 
1756 //           glued edges and loops) described by the list of indexes of 
1757 //           its edges in the structure.
1758 //           (negative index means reversed edge)
1759 //=======================================================================
1760 void BRepMesh_Delaun::meshSimplePolygon( TColStd_SequenceOfInteger&     thePolygon,
1761                                          NCollection_Sequence<Bnd_B2d>& thePolyBoxes )
1762 {
1763   // Check is the given polygon elementary
1764   if ( meshElementaryPolygon( thePolygon ) )
1765     return;
1766
1767
1768   // Polygon contains more than 3 links
1769   Standard_Integer aFirstEdgeInfo = thePolygon(1);
1770   const BRepMesh_Edge& aFirstEdge = GetEdge( Abs( aFirstEdgeInfo ) );
1771
1772   Standard_Integer aNodes[3];
1773   getOrientedNodes( aFirstEdge, aFirstEdgeInfo > 0, aNodes );
1774   
1775   gp_Pnt2d aRefVertices[3];
1776   aRefVertices[0] = GetVertex( aNodes[0] ).Coord();
1777   aRefVertices[1] = GetVertex( aNodes[1] ).Coord();
1778
1779   gp_Vec2d aRefEdgeDir( aRefVertices[0], aRefVertices[1] );
1780
1781   Standard_Real aRefEdgeLen = aRefEdgeDir.Magnitude();
1782   if ( aRefEdgeLen < Precision )
1783     return;
1784
1785   aRefEdgeDir /= aRefEdgeLen;
1786
1787   // Find a point with minimum distance respect
1788   // the end of reference link
1789   Standard_Integer aUsedLinkId = 0;
1790   Standard_Real    aOptAngle   = 0.0;
1791   Standard_Real    aMinDist    = RealLast();
1792   Standard_Integer aPivotNode  = aNodes[1];
1793   Standard_Integer aPolyLen    = thePolygon.Length();
1794   for ( Standard_Integer aLinkIt = 3; aLinkIt <= aPolyLen; ++aLinkIt )
1795   {
1796     Standard_Integer aLinkInfo = thePolygon( aLinkIt );
1797     const BRepMesh_Edge& aNextEdge = GetEdge( Abs( aLinkInfo ) );
1798
1799     aPivotNode = aLinkInfo > 0 ? 
1800       aNextEdge.FirstNode() : 
1801       aNextEdge.LastNode();
1802
1803     gp_Pnt2d aPivotVertex = GetVertex( aPivotNode ).Coord();
1804     gp_Vec2d aDistanceDir( aRefVertices[1], aPivotVertex );
1805
1806     Standard_Real aDist     = aRefEdgeDir ^ aDistanceDir;
1807     Standard_Real aAngle    = Abs( aRefEdgeDir.Angle(aDistanceDir) );
1808     Standard_Real anAbsDist = Abs( aDist );
1809     if ( ( anAbsDist < Precision                 ) ||
1810          ( myIsPositiveOrientation && aDist < 0. ) ||
1811          (!myIsPositiveOrientation && aDist > 0. ) )
1812     {
1813       continue;
1814     }
1815
1816     if ( ( anAbsDist >= aMinDist                             ) &&
1817          ( aAngle <= aOptAngle || aAngle > AngDeviation90Deg ) )
1818     {
1819       continue;
1820     }
1821
1822     // Check is the test link crosses the polygon boudaries
1823     Standard_Boolean isIntersect = Standard_False;
1824     for ( Standard_Integer aRefLinkNodeIt = 0; aRefLinkNodeIt < 2; ++aRefLinkNodeIt )
1825     {
1826       const Standard_Integer& aLinkFirstNode   = aNodes[aRefLinkNodeIt];
1827       const gp_Pnt2d&         aLinkFirstVertex = aRefVertices[aRefLinkNodeIt];
1828
1829       Bnd_B2d aBox;
1830       aBox.Add( aLinkFirstVertex );
1831       aBox.Add( aPivotVertex );
1832
1833       BRepMesh_Edge aCheckLink( aLinkFirstNode, aPivotNode, BRepMesh_Free );
1834
1835       Standard_Integer aCheckLinkIt = 2;
1836       for ( ; aCheckLinkIt <= aPolyLen; ++aCheckLinkIt )
1837       {
1838         if( aCheckLinkIt == aLinkIt )
1839           continue;
1840         
1841         if ( !aBox.IsOut( thePolyBoxes.Value( aCheckLinkIt ) ) )
1842         {
1843           const BRepMesh_Edge& aPolyLink = 
1844             GetEdge( Abs( thePolygon( aCheckLinkIt ) ) );
1845
1846           if ( aCheckLink.IsEqual( aPolyLink ) )
1847             continue;
1848
1849           // intersection is possible...                  
1850           gp_Pnt2d anIntPnt;
1851           IntFlag aIntFlag = intSegSeg( aCheckLink, aPolyLink, 
1852             Standard_False, Standard_False, anIntPnt );
1853
1854           if( aIntFlag != BRepMesh_Delaun::NoIntersection )
1855           {
1856             isIntersect = Standard_True;
1857             break;
1858           }
1859         }
1860       }
1861
1862       if ( isIntersect )
1863         break;
1864     }
1865
1866     if( isIntersect )
1867       continue;
1868
1869
1870     aOptAngle       = aAngle;
1871     aMinDist        = anAbsDist;
1872     aNodes[2]       = aPivotNode;
1873     aRefVertices[2] = aPivotVertex;
1874     aUsedLinkId     = aLinkIt;
1875   }
1876
1877   if ( aUsedLinkId == 0 )
1878     return;
1879
1880
1881   BRepMesh_Edge aNewEdges[2] = {
1882     BRepMesh_Edge( aNodes[1], aNodes[2], BRepMesh_Free ),
1883     BRepMesh_Edge( aNodes[2], aNodes[0], BRepMesh_Free ) };
1884
1885   Standard_Integer aNewEdgesInfo[3] = {
1886     aFirstEdgeInfo,
1887     myMeshData->AddLink( aNewEdges[0] ),
1888     myMeshData->AddLink( aNewEdges[1] ) };
1889
1890
1891   Standard_Integer anEdges[3];
1892   Standard_Boolean anEdgesOri[3];
1893   for ( Standard_Integer aTriEdgeIt = 0; aTriEdgeIt < 3; ++aTriEdgeIt )
1894   {
1895     const Standard_Integer& anEdgeInfo = aNewEdgesInfo[aTriEdgeIt];
1896     anEdges[aTriEdgeIt]    = Abs( anEdgeInfo );
1897     anEdgesOri[aTriEdgeIt] = anEdgeInfo > 0;
1898   }
1899   addTriangle( anEdges, anEdgesOri, aNodes );
1900
1901   // Create triangle and split the source polygon on two 
1902   // parts (if possible) and mesh each part as independent
1903   // polygon.
1904   if ( aUsedLinkId < aPolyLen )
1905   {
1906     TColStd_SequenceOfInteger aRightPolygon;
1907     thePolygon.Split( aUsedLinkId, aRightPolygon );
1908     aRightPolygon.Prepend( -aNewEdgesInfo[2] );
1909
1910     NCollection_Sequence<Bnd_B2d> aRightPolyBoxes;
1911     thePolyBoxes.Split( aUsedLinkId, aRightPolyBoxes );
1912
1913     Bnd_B2d aBox;
1914     aBox.Add( aRefVertices[0] );
1915     aBox.Add( aRefVertices[2] );
1916     aRightPolyBoxes.Prepend( aBox );
1917
1918     meshSimplePolygon( aRightPolygon, aRightPolyBoxes );
1919   }
1920   else
1921   {
1922     thePolygon.Remove  ( aPolyLen );
1923     thePolyBoxes.Remove( aPolyLen );
1924   }
1925
1926   if ( aUsedLinkId > 3 )
1927   {
1928     thePolygon.SetValue( 1, -aNewEdgesInfo[1] );
1929
1930     Bnd_B2d aBox;
1931     aBox.Add( aRefVertices[1] );
1932     aBox.Add( aRefVertices[2] );
1933
1934     thePolyBoxes.SetValue( 1, aBox );
1935
1936     meshSimplePolygon( thePolygon, thePolyBoxes );
1937   }
1938 }
1939
1940 //=======================================================================
1941 //function : RemoveVertex
1942 //purpose  : Removes a vertex from the triangulation
1943 //=======================================================================
1944 void  BRepMesh_Delaun::RemoveVertex( const BRepMesh_Vertex& theVertex )
1945 {
1946   BRepMesh_SelectorOfDataStructureOfDelaun aSelector( myMeshData );
1947   aSelector.NeighboursOf( theVertex );
1948
1949   BRepMesh_MapOfIntegerInteger aLoopEdges( 10, myMeshData->Allocator() );
1950
1951   // Loop on triangles to be destroyed :
1952   BRepMesh_MapOfInteger::Iterator aTriangleIt( aSelector.Elements() );
1953   for ( ; aTriangleIt.More(); aTriangleIt.Next() )
1954     deleteTriangle( aTriangleIt.Key(), aLoopEdges );
1955
1956   NCollection_Sequence<Bnd_B2d> aBoxes;
1957   TColStd_SequenceOfInteger     aPolygon;
1958   Standard_Integer aLoopEdgesCount = aLoopEdges.Extent();
1959   BRepMesh_MapOfIntegerInteger::Iterator aLoopEdgesIt( aLoopEdges );
1960
1961   if ( aLoopEdgesIt.More() )
1962   {
1963     const BRepMesh_Edge& anEdge = GetEdge( aLoopEdgesIt.Key() );
1964     Standard_Integer aFirstNode = anEdge.FirstNode();
1965     Standard_Integer aLastNode;
1966     Standard_Integer aPivotNode = anEdge.LastNode();
1967     Standard_Integer anEdgeId   = aLoopEdgesIt.Key();
1968     
1969     Standard_Boolean isPositive = (Standard_Boolean)aLoopEdges( anEdgeId );
1970     if ( !isPositive )
1971     {
1972       Standard_Integer aTmp;
1973       aTmp       = aFirstNode;
1974       aFirstNode = aPivotNode;
1975       aPivotNode = aTmp;
1976       
1977       aPolygon.Append( -anEdgeId );
1978     }
1979     else
1980       aPolygon.Append( anEdgeId );
1981
1982     fillBndBox( aBoxes, GetVertex( aFirstNode ), GetVertex( aPivotNode ) );
1983
1984     aLoopEdges.UnBind( anEdgeId );
1985     
1986     aLastNode = aFirstNode;
1987     while ( aPivotNode != aLastNode )
1988     {
1989       BRepMesh_ListOfInteger::Iterator aLinkIt( myMeshData->LinkNeighboursOf( aPivotNode ) );
1990       for ( ; aLinkIt.More(); aLinkIt.Next() )
1991       {
1992         if ( aLinkIt.Value() != anEdgeId &&
1993              aLoopEdges.IsBound( aLinkIt.Value() ) )
1994         {
1995           Standard_Integer aCurrentNode;
1996           anEdgeId = aLinkIt.Value();
1997           const BRepMesh_Edge& anEdge1 = GetEdge( anEdgeId );
1998           
1999           aCurrentNode = anEdge1.LastNode();
2000           if ( aCurrentNode != aPivotNode )
2001           {
2002             aCurrentNode = anEdge1.FirstNode();
2003             aPolygon.Append( -anEdgeId );
2004           }
2005           else
2006             aPolygon.Append( anEdgeId );
2007
2008           fillBndBox( aBoxes, GetVertex( aCurrentNode ), GetVertex( aPivotNode ) );
2009             
2010           aPivotNode = aCurrentNode;
2011           aLoopEdges.UnBind( anEdgeId );
2012           break;
2013         }
2014       }
2015       
2016       if ( aLoopEdgesCount <= 0 )
2017         break;
2018       --aLoopEdgesCount;
2019     }
2020     
2021     meshPolygon( aPolygon, aBoxes );
2022   }
2023 }
2024
2025
2026 //=======================================================================
2027 //function : AddVertices
2028 //purpose  : Adds some vertices in the triangulation.
2029 //=======================================================================
2030 void  BRepMesh_Delaun::AddVertices( BRepMesh_Array1OfVertexOfDelaun& theVertices )
2031 {
2032   BRepMesh_HeapSortVertexOfDelaun::Sort( theVertices, 
2033     BRepMesh_ComparatorOfVertexOfDelaun( SortingDirection, Precision ) );
2034
2035   Standard_Integer aLower  = theVertices.Lower();
2036   Standard_Integer anUpper = theVertices.Upper();
2037     
2038   TColStd_Array1OfInteger aVertexIndexes( aLower, anUpper );
2039   for ( Standard_Integer i = aLower; i <= anUpper; ++i )     
2040     aVertexIndexes(i) = myMeshData->AddNode( theVertices(i) );
2041
2042   createTrianglesOnNewVertices( aVertexIndexes );
2043 }
2044
2045 //=======================================================================
2046 //function : UseEdge
2047 //purpose  : Modify mesh to use the edge. Return True if done
2048 //=======================================================================
2049 Standard_Boolean BRepMesh_Delaun::UseEdge( const Standard_Integer theIndex )
2050 {
2051   /*
2052   const BRepMesh_PairOfIndex& aPair = myMeshData->ElemConnectedTo( theIndex );
2053   if ( aPair.Extent() == 0 )
2054   {
2055     const BRepMesh_Edge& anEdge = GetEdge( theIndex );
2056     
2057     Standard_Integer aStartNode, aPivotNode, anOtherNode;
2058     aStartNode = anEdge.FirstNode();
2059     aPivotNode = anEdge.LastNode();
2060     
2061     const BRepMesh_ListOfInteger& aStartNodeNeighbors = myMeshData->LinkNeighboursOf( aStartNode );
2062     const BRepMesh_ListOfInteger& aPivotNodeNeighbors = myMeshData->LinkNeighboursOf( aPivotNode );
2063     
2064     if ( aStartNodeNeighbors.Extent() > 0 &&
2065          aPivotNodeNeighbors.Extent() > 0 )
2066     {
2067       const BRepMesh_Vertex& aStartVertex = GetVertex( aStartNode );
2068       const BRepMesh_Vertex& aPivotVertex = GetVertex( aPivotNode );
2069
2070       gp_XY aVEdge   ( aPivotVertex.Coord() );
2071       aVEdge.Subtract( aStartVertex.Coord() );
2072
2073       Standard_Real    anAngle    = 0.;
2074       Standard_Real    anAngleMin = RealLast();
2075       Standard_Real    anAngleMax = RealFirst();
2076       Standard_Integer aLeftEdge  = 0, aRightEdge = 0;
2077
2078       BRepMesh_ListOfInteger::Iterator aNeighborIt( aPivotNodeNeighbors );
2079       for ( ; aNeighborIt.More(); aNeighborIt.Next() )
2080       {
2081         Standard_Integer anEdgeId = aNeighborIt.Value();
2082         if ( anEdgeId != theIndex )
2083         {
2084           const BRepMesh_Edge& aNextEdge = GetEdge( anEdgeId );
2085
2086           Standard_Boolean isInMesh = Standard_True;
2087           if ( aNextEdge.Movability() == BRepMesh_Free )
2088           {
2089             if ( myMeshData->ElemConnectedTo( anEdgeId ).IsEmpty() ) 
2090               isInMesh = Standard_False;
2091           }
2092
2093           if ( isInMesh )
2094           {
2095             anOtherNode = aNextEdge.FirstNode();
2096             if ( anOtherNode == aPivotNode )
2097               anOtherNode = aNextEdge.LastNode();
2098
2099             gp_XY aVEdgeCur = GetVertex( anOtherNode ).Coord();
2100             aVEdgeCur.Subtract( aPivotVertex.Coord() );
2101
2102             anAngle = gp_Vec2d( aVEdge ).Angle( gp_Vec2d( aVEdgeCur ) );
2103           }
2104           
2105           if ( anAngle > anAngleMax )
2106           {
2107             anAngleMax = anAngle;
2108             aLeftEdge  = anEdgeId;
2109           }
2110           if ( anAngle < anAngleMin )
2111           {
2112             anAngleMin = anAngle;
2113             aRightEdge = anEdgeId;
2114           }
2115         }
2116       }
2117       
2118       if ( aLeftEdge > 0 )
2119       {
2120         if (aLeftEdge==aRightEdge)
2121         {
2122         }
2123         else
2124         {
2125         }
2126       }
2127     }
2128   }
2129   */
2130   return Standard_False;
2131 }
2132
2133 //=======================================================================
2134 //function : getEdgesByType
2135 //purpose  : Gives the list of edges with type defined by input parameter
2136 //=======================================================================
2137 Handle(BRepMesh_MapOfInteger) BRepMesh_Delaun::getEdgesByType(
2138   const BRepMesh_DegreeOfFreedom theEdgeType ) const
2139 {
2140   Handle(BRepMesh_MapOfInteger) aResult = new BRepMesh_MapOfInteger;
2141   BRepMesh_MapOfInteger::Iterator anEdgeIt( myMeshData->LinkOfDomain() );
2142
2143   for ( ; anEdgeIt.More(); anEdgeIt.Next() )
2144   {
2145     Standard_Integer anEdge = anEdgeIt.Key();
2146     Standard_Boolean isToAdd = (theEdgeType == BRepMesh_Free) ?
2147       (myMeshData->ElemConnectedTo( anEdge ).Extent() <= 1) :
2148       (GetEdge( anEdge ).Movability() == theEdgeType);
2149
2150     if (isToAdd)
2151       aResult->Add( anEdge );
2152   }
2153   
2154   return aResult;
2155 }
2156
2157 //=======================================================================
2158 //function : calculateDist
2159 //purpose  : Calculates distances between the given point and edges of
2160 //           triangle
2161 //=======================================================================
2162 Standard_Real BRepMesh_Delaun::calculateDist( const gp_XY            theVEdges[3],
2163                                               const gp_XY            thePoints[3],
2164                                               const Standard_Integer theEdgesId[3],
2165                                               const BRepMesh_Vertex& theVertex,
2166                                               Standard_Real          theDistance[3],
2167                                               Standard_Real          theSqModulus[3],
2168                                               Standard_Integer&      theEdgeOn ) const
2169 {
2170   Standard_Real aMinDist = -1;
2171   if ( !theVEdges || !thePoints || !theEdgesId || 
2172        !theDistance || !theSqModulus )
2173     return aMinDist;
2174     
2175   for( Standard_Integer i = 0; i < 3; ++i )
2176   {
2177     theSqModulus[i] = theVEdges[i].SquareModulus();
2178     if ( theSqModulus[i] <= Precision2 )
2179       return -1;
2180       
2181     theDistance[i] = theVEdges[i] ^ ( theVertex.Coord() - thePoints[i] );
2182     
2183     Standard_Real aDist = theDistance[i] * theDistance[i];
2184     aDist /= theSqModulus[i];
2185     
2186     if ( aMinDist < 0 || aDist < aMinDist )
2187     {
2188       theEdgeOn = theEdgesId[i];
2189       aMinDist  = aDist;
2190     }
2191   }
2192   
2193   return aMinDist;
2194 }
2195
2196 //=======================================================================
2197 //function : Contains
2198 //purpose  : Test if triangle of index <TrianIndex> contains geometricaly
2199 //           <theVertex>. If <theEdgeOn> is != 0 then theVertex is on Edge 
2200 //           of  index <theEdgeOn>
2201 //=======================================================================
2202 Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId,
2203                                             const BRepMesh_Vertex& theVertex,
2204                                             Standard_Integer&      theEdgeOn ) const
2205 {
2206   theEdgeOn = 0;
2207   
2208   Standard_Integer e[3];
2209   Standard_Boolean o[3];
2210   Standard_Integer p[3];
2211
2212   GetTriangle( theTriangleId ).Edges( e[0], e[1], e[2],
2213                                       o[0], o[1], o[2] );
2214                                       
2215   const BRepMesh_Edge* anEdges[3] = { &GetEdge( e[0] ),
2216                                       &GetEdge( e[1] ),
2217                                       &GetEdge( e[2] ) };
2218   if ( o[0] )
2219   {
2220     p[0] = anEdges[0]->FirstNode();
2221     p[1] = anEdges[0]->LastNode();
2222   }
2223   else
2224   {
2225     p[1] = anEdges[0]->FirstNode();
2226     p[0] = anEdges[0]->LastNode();
2227   }
2228   
2229   if ( o[2] )
2230     p[2] = anEdges[2]->FirstNode();
2231   else
2232     p[2] = anEdges[2]->LastNode();
2233
2234   gp_XY aPoints[3];
2235   aPoints[0] = GetVertex( p[0] ).Coord();
2236   aPoints[1] = GetVertex( p[1] ).Coord();
2237   aPoints[2] = GetVertex( p[2] ).Coord();
2238   
2239   gp_XY aVEdges[3];
2240   aVEdges[0] = aPoints[1]; 
2241   aVEdges[0].Subtract( aPoints[0] );
2242   
2243   aVEdges[1] = aPoints[2]; 
2244   aVEdges[1].Subtract( aPoints[1] );
2245   
2246   aVEdges[2] = aPoints[0];
2247   aVEdges[2].Subtract( aPoints[2] );
2248
2249   Standard_Real aDistance[3];
2250   Standard_Real aSqModulus[3];
2251
2252   Standard_Real aMinDist;  
2253   aMinDist = calculateDist( aVEdges, aPoints, e, theVertex, aDistance, aSqModulus, theEdgeOn );
2254   if ( aMinDist < 0 )
2255     return Standard_False;
2256       
2257   if ( aMinDist > Precision2 )
2258   {
2259     Standard_Integer anEdgeId = theEdgeOn;
2260     theEdgeOn = 0;
2261     
2262     if ( anEdgeId != 0 ) 
2263     {
2264       Standard_Integer i = 0;
2265       for ( ; i < 3; ++i )
2266       {
2267         if( e[i] == anEdgeId )
2268           break;
2269       }
2270       
2271       if( anEdges[i]->Movability() != BRepMesh_Free )
2272         if ( aDistance[i] < ( aSqModulus[i] / 5. ) )
2273           theEdgeOn = e[i];
2274     }
2275   }
2276
2277   return ( aDistance[0] + aDistance[1] + aDistance[2] != 0. &&
2278             ( ( aDistance[0] >= 0. && aDistance[1] >= 0. && aDistance[2] >= 0. ) ||
2279               ( aDistance[0] <= 0. && aDistance[1] <= 0. && aDistance[2] <= 0. ) ) );
2280 }
2281
2282
2283 //=============================================================================
2284 //function : classifyPoint
2285 //purpose  : Classifies the point in case of coincidence of two vectors.
2286 //           Returns zero value if point is out of segment and non zero 
2287 //           value if point is between the first and the second point of segment.
2288 //           thePoint1       - the start point of a segment (base point)
2289 //           thePoint2       - the end point of a segment
2290 //           thePointToCheck - the point to classify
2291 //=============================================================================
2292 Standard_Integer BRepMesh_Delaun::classifyPoint( const gp_XY& thePoint1,
2293                                                  const gp_XY& thePoint2,
2294                                                  const gp_XY& thePointToCheck ) const
2295 {
2296   gp_XY aP1 = thePoint2       - thePoint1;
2297   gp_XY aP2 = thePointToCheck - thePoint1;
2298   
2299   Standard_Real aDist = Abs( aP1 ^ aP2 );
2300   if ( aDist >= Precision )
2301   {
2302     aDist = ( aDist * aDist ) / aP1.SquareModulus();
2303     if ( aDist >= Precision2 )
2304       return 0; //out
2305   }
2306     
2307   gp_XY aMult = aP1.Multiplied( aP2 );
2308   if ( aMult.X() < 0.0 || aMult.Y() < 0.0 )
2309     return 0; //out
2310     
2311   if ( aP1.SquareModulus() < aP2.SquareModulus() )
2312     return 0; //out
2313     
2314   if ( thePointToCheck.IsEqual( thePoint1, Precision ) || 
2315        thePointToCheck.IsEqual( thePoint2, Precision ) )
2316     return -1; //coinsides with an end point
2317     
2318   return 1;
2319 }
2320
2321 //=============================================================================
2322 //function : intSegSeg
2323 //purpose  : Checks intersection between the two segments.
2324 //=============================================================================
2325 BRepMesh_Delaun::IntFlag BRepMesh_Delaun::intSegSeg( const BRepMesh_Edge&   theEdg1,
2326                                                      const BRepMesh_Edge&   theEdg2,
2327                                                      const Standard_Boolean isConsiderEndPointTouch,
2328                                                      const Standard_Boolean isConsiderPointOnEdge,
2329                                                      gp_Pnt2d&              theIntPnt) const
2330 {
2331   gp_XY p1, p2, p3, p4;
2332   p1 = GetVertex( theEdg1.FirstNode() ).Coord();
2333   p2 = GetVertex( theEdg1.LastNode()  ).Coord();
2334   p3 = GetVertex( theEdg2.FirstNode() ).Coord();
2335   p4 = GetVertex( theEdg2.LastNode()  ).Coord();
2336   
2337   Standard_Integer aPoint1 = classifyPoint( p1, p2, p3 );
2338   Standard_Integer aPoint2 = classifyPoint( p1, p2, p4 );
2339   Standard_Integer aPoint3 = classifyPoint( p3, p4, p1 );
2340   Standard_Integer aPoint4 = classifyPoint( p3, p4, p2 );
2341
2342   // Consider case when edges have shared vertex
2343   if ( isConsiderEndPointTouch )
2344   {
2345     if ( aPoint1 < 0 || aPoint2 < 0 )
2346       return BRepMesh_Delaun::EndPointTouch;
2347   }
2348
2349   Standard_Integer aPosHash = 
2350     aPoint1 + aPoint2 + aPoint3 + aPoint4;
2351
2352   /*=========================================*/
2353   /*  1) hash code == 1:
2354
2355                     0+
2356                     /
2357            0      1/         0
2358            +======+==========+
2359   
2360       2) hash code == 2:
2361
2362            0    1        1   0
2363         a) +----+========+---+
2364
2365            0       1   1     0
2366         b) +-------+===+=====+
2367
2368                                              */
2369   /*=========================================*/
2370   if ( aPosHash == 1 )
2371   {
2372     return isConsiderPointOnEdge ? 
2373       BRepMesh_Delaun::PointOnEdge :
2374       BRepMesh_Delaun::NoIntersection;
2375   }
2376   else if ( aPosHash == 2 )
2377     return BRepMesh_Delaun::Glued;
2378
2379   gp_XY aVec1           = p2 - p1;
2380   gp_XY aVec2           = p4 - p3;
2381   gp_XY aVecStartPoints = p3 - p1;
2382     
2383   Standard_Real aCrossD1D2 = aVec1           ^ aVec2;
2384   Standard_Real aCrossD1D3 = aVecStartPoints ^ aVec2;
2385
2386   // is edgegs codirectional
2387   if ( Abs( aCrossD1D2 ) < Precision )
2388   {
2389     // just a parallel case?
2390     if( Abs( aCrossD1D3 ) < Precision )
2391     {
2392       /*=========================================*/
2393       /*  Here the following cases are possible:
2394           1) hash code == -4:
2395
2396                -1                -1
2397                 +=================+
2398                -1                -1
2399
2400           2) hash code == -2:
2401
2402                 0       -1        0
2403                 +--------+========+
2404                         -1
2405
2406           3) hash code == -1:
2407
2408                 0        1        -1
2409                 +--------+========+
2410                                   -1
2411
2412           4) hash code == 0:
2413
2414                 0      0  0       0
2415                 +------+  +=======+
2416                 0      0  0       0
2417                                                  */
2418       /*=========================================*/
2419
2420       if ( aPosHash < -2 )
2421         return BRepMesh_Delaun::Same;
2422       else if ( aPosHash == -1 )
2423         return BRepMesh_Delaun::Glued;
2424
2425       return BRepMesh_Delaun::NoIntersection;
2426     }
2427     else
2428       return BRepMesh_Delaun::NoIntersection;
2429   }
2430
2431   Standard_Real aPar = aCrossD1D3 / aCrossD1D2;
2432   // inrersects out of first segment range
2433   if( aPar < Precision || aPar > EndPrecision )
2434     return BRepMesh_Delaun::NoIntersection;
2435  
2436   Standard_Real aCrossD2D3 = aVecStartPoints.Reversed() ^ aVec1;
2437   aPar = aCrossD2D3 / -aCrossD1D2;
2438   // inrersects out of second segment range
2439   if( aPar < Precision || aPar > EndPrecision )
2440     return BRepMesh_Delaun::NoIntersection;
2441  
2442   theIntPnt = p3 + aPar * aVec2;
2443   return BRepMesh_Delaun::Cross;
2444 }
2445
2446 //=============================================================================
2447 //function : polyArea
2448 //purpose  : Returns area of the loop of the given polygon defined by indices 
2449 //           of its start and end links.
2450 //=============================================================================
2451 Standard_Real BRepMesh_Delaun::polyArea( const TColStd_SequenceOfInteger& thePolygon,
2452                                          const Standard_Integer           theStartIndex,
2453                                          const Standard_Integer           theEndIndex ) const
2454 {
2455   Standard_Real aArea = 0.0;
2456   Standard_Integer aPolyLen = thePolygon.Length();
2457   if ( theStartIndex >= theEndIndex ||
2458        theStartIndex > aPolyLen )
2459   {
2460     return aArea;
2461   }
2462
2463   Standard_Integer aEndIndex = (theEndIndex > aPolyLen) ? 
2464     aPolyLen : theEndIndex;
2465
2466   Standard_Integer aCurEdgeInfo = thePolygon( theStartIndex );
2467   Standard_Integer aCurEdgeId   = Abs( aCurEdgeInfo );
2468   const BRepMesh_Edge* aCurEdge = &GetEdge( aCurEdgeId );
2469
2470   Standard_Integer aNodes[2];
2471   getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aNodes );
2472
2473   gp_Pnt2d aRefPnt = GetVertex( aNodes[0] ).Coord();
2474   Standard_Integer aPolyIt = theStartIndex + 1;
2475   for ( ; aPolyIt <= theEndIndex; ++aPolyIt )
2476   {
2477     aCurEdgeInfo = thePolygon( aPolyIt );
2478     aCurEdgeId   = Abs( aCurEdgeInfo );
2479     aCurEdge     = &GetEdge( aCurEdgeId );
2480
2481     getOrientedNodes( *aCurEdge, aCurEdgeInfo > 0, aNodes );
2482     gp_Vec2d aVec1( aRefPnt, GetVertex( aNodes[0] ).Coord() );
2483     gp_Vec2d aVec2( aRefPnt, GetVertex( aNodes[1] ).Coord() );
2484
2485     aArea += aVec1 ^ aVec2;
2486   }
2487
2488   return aArea / 2.;
2489 }