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