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