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