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