Salome HOME
Update copyrights
[modules/smesh.git] / src / SMESHUtils / SMESH_MAT2d.cxx
1 // Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22 // File      : SMESH_MAT2d.cxx
23 // Created   : Thu May 28 17:49:53 2015
24 // Author    : Edward AGAPOV (eap)
25
26 #include "SMESH_MAT2d.hxx"
27
28 #include <list>
29
30 #include <BRepAdaptor_CompCurve.hxx>
31 #include <BRepBuilderAPI_MakeEdge.hxx>
32 #include <BRepBuilderAPI_MakeVertex.hxx>
33 #include <BRep_Builder.hxx>
34 #include <BRep_Tool.hxx>
35 #include <Bnd_B2d.hxx>
36 //#include <GCPnts_AbscissaPoint.hxx>
37 #include <GCPnts_TangentialDeflection.hxx>
38 // #include <GCPnts_UniformAbscissa.hxx>
39 // #include <GCPnts_UniformDeflection.hxx>
40 #include <Geom2d_Curve.hxx>
41 //#include <GeomAdaptor_Curve.hxx>
42 #include <Geom2dAdaptor_Curve.hxx>
43 #include <Geom_Curve.hxx>
44 #include <Geom_Surface.hxx>
45 #include <TopExp.hxx>
46 #include <TopoDS_Vertex.hxx>
47 #include <TopoDS_Wire.hxx>
48
49 #ifdef _DEBUG_
50 //#define _MYDEBUG_
51 #include "SMESH_File.hxx"
52 #include "SMESH_Comment.hxx"
53 #endif
54
55 using namespace std;
56 using boost::polygon::x;
57 using boost::polygon::y;
58 using SMESH_MAT2d::TVD;
59 using SMESH_MAT2d::TVDEdge;
60 using SMESH_MAT2d::TVDCell;
61 using SMESH_MAT2d::TVDVertex;
62
63 namespace
64 {
65   // Input data for construct_voronoi()
66   // -------------------------------------------------------------------------------------
67
68   struct InPoint
69   {
70     int _a, _b;    // coordinates
71     double _param; // param on EDGE
72     InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
73     InPoint() : _a(0), _b(0), _param(0) {}
74
75     // working data
76     list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order
77
78     size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
79     bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
80     bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
81                                                            Abs( _b - v->y() ) < 1. ); }
82   };
83   // -------------------------------------------------------------------------------------
84
85   struct InSegment
86   {
87     InPoint * _p0;
88     InPoint * _p1;
89
90     // working data
91     size_t                 _geomEdgeInd; // EDGE index within the FACE
92     const TVDCell*         _cell;
93     list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
94
95     InSegment( InPoint * p0, InPoint * p1, size_t iE)
96       : _p0(p0), _p1(p1), _geomEdgeInd(iE), _cell(0) {}
97     InSegment() : _p0(0), _p1(0), _geomEdgeInd(0), _cell(0) {}
98
99     const InPoint& point0() const { return *_p0; }
100     const InPoint& point1() const { return *_p1; }
101
102     inline bool isConnected( const TVDEdge* edge );
103
104     inline bool isExternal( const TVDEdge* edge );
105
106     static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); }
107
108     static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); }
109   };
110
111   // check  if a TVDEdge begins at my end or ends at my start
112   inline bool InSegment::isConnected( const TVDEdge* edge )
113   {
114     return (( edge->vertex0() && edge->vertex1() )
115             &&
116             ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
117               Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
118              (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
119               Abs( edge->vertex1()->y() - _p0->_b ) < 1.  )));
120   }
121
122   // check if a MA TVDEdge is outside of a domain
123   inline bool InSegment::isExternal( const TVDEdge* edge )
124   {
125     double dot = // x1*x2 + y1*y2;   (x1,y1) - internal normal of InSegment
126       ( _p0->_b - _p1->_b ) * ( 0.5 * ( edge->vertex0()->x() + edge->vertex1()->x() ) - _p0->_a ) +
127       ( _p1->_a - _p0->_a ) * ( 0.5 * ( edge->vertex0()->y() + edge->vertex1()->y() ) - _p0->_b );
128     return dot < 0.;
129   }
130
131   // // -------------------------------------------------------------------------------------
132   // const size_t theExternMA = 111; // to mark external MA edges
133
134   // bool isExternal( const TVDEdge* edge )
135   // {
136   //   return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA );
137   // }
138
139   // // mark external MA edges
140   // void markExternalEdges( const TVDEdge* edge )
141   // {
142   //   if ( isExternal( edge ))
143   //     return;
144   //   SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge );
145   //   SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() );
146   //   if ( edge->is_primary() && edge->vertex1() )
147   //   {
148   //     const TVDVertex * v = edge->vertex1();
149   //     edge = v->incident_edge();
150   //     do {
151   //       markExternalEdges( edge );
152   //       edge = edge->rot_next();
153   //     } while ( edge != v->incident_edge() );
154   //   }
155   // }
156
157   // -------------------------------------------------------------------------------------
158 #ifdef _MYDEBUG_
159   // writes segments into a txt file readable by voronoi_visualizer
160   void inSegmentsToFile( vector< InSegment>& inSegments)
161   {
162     if ( inSegments.size() > 1000 )
163       return;
164     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
165     const char*     user = getenv("USER");
166     if ( !user || strcmp( user, "eap" )) return;
167     SMESH_File file(fileName, false );
168     file.remove();
169     file.openForWriting();
170     SMESH_Comment text;
171     text << "0\n"; // nb points
172     text << inSegments.size() << "\n"; // nb segments
173     for ( size_t i = 0; i < inSegments.size(); ++i )
174     {
175       text << inSegments[i]._p0->_a << " "
176            << inSegments[i]._p0->_b << " "
177            << inSegments[i]._p1->_a << " "
178            << inSegments[i]._p1->_b << "\n";
179     }
180     text << "\n";
181     file.write( text.c_str(), text.size() );
182     cout << "Write " << fileName << endl;
183   }
184   void dumpEdge( const TVDEdge* edge )
185   {
186     cout << "*Edge_" << edge;
187     if ( !edge->vertex0() )
188       cout << " ( INF, INF";
189     else
190       cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y();
191     if ( !edge->vertex1() )
192       cout << ") -> ( INF, INF";
193     else
194       cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
195     cout << ")\t cell=" << edge->cell()
196          << " iBnd=" << edge->color()
197          << " twin=" << edge->twin()
198          << " twin_cell=" << edge->twin()->cell()
199          << " prev=" << edge->prev() << " next=" << edge->next()
200          << ( edge->is_primary() ? " MA " : " SCND" )
201          << ( edge->is_linear() ? " LIN " : " CURV" )
202          << endl;
203   }
204   void dumpCell( const TVDCell* cell )
205   {
206     cout << "**Cell_" << cell << " GEOM=" << cell->color() << " ";
207     cout << ( cell->contains_segment() ? " SEG " : " PNT " );
208     if ( cell-> is_degenerate() )
209       cout << " degen ";
210     else
211     {
212       cout << endl;
213       const TVDEdge* edge = cell->incident_edge();
214       size_t i = 0;
215       do {
216         edge = edge->next();
217         cout << "   - " << ++i << " ";
218         dumpEdge( edge );
219       } while (edge != cell->incident_edge());
220     }
221   }
222 #else
223   #define inSegmentsToFile(arg) {}
224   //void dumpEdge( const TVDEdge* edge ) {}
225   //void dumpCell( const TVDCell* cell ) {}
226 #endif
227 }
228 // -------------------------------------------------------------------------------------
229
230 namespace boost {
231   namespace polygon {
232
233     template <>
234     struct geometry_concept<InPoint> {
235       typedef point_concept type;
236     };
237     template <>
238     struct point_traits<InPoint> {
239       typedef int coordinate_type;
240
241       static inline coordinate_type get(const InPoint& point, orientation_2d orient) {
242         return (orient == HORIZONTAL) ? point._a : point._b;
243       }
244     };
245
246     template <>
247     struct geometry_concept<InSegment> {
248       typedef segment_concept type;
249     };
250
251     template <>
252     struct segment_traits<InSegment> {
253       typedef int coordinate_type;
254       typedef InPoint point_type;
255
256       static inline point_type get(const InSegment& segment, direction_1d dir) {
257         return *(dir.to_int() ? segment._p1 : segment._p0);
258       }
259     };
260   }  // namespace polygon
261 } // namespace boost
262   // -------------------------------------------------------------------------------------
263
264 namespace
265 {
266   const int    theNoBrachID = 0;
267   double       theScale[2]; // scale used in bndSegsToMesh()
268   const size_t theNoEdgeID = std::numeric_limits<size_t>::max() / 1000;
269
270   // -------------------------------------------------------------------------------------
271   /*!
272    * \brief Intermediate DS to create InPoint's
273    */
274   struct UVU
275   {
276     gp_Pnt2d _uv;
277     double   _u;
278     UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {}
279     InPoint getInPoint( double scale[2] )
280     {
281       return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u );
282     }
283   };
284   // -------------------------------------------------------------------------------------
285   /*!
286    * \brief Segment of EDGE, used to create BndPoints
287    */
288   struct BndSeg
289   {
290     InSegment*       _inSeg;
291     const TVDEdge*   _edge;
292     double           _uLast;
293     BndSeg*          _prev; // previous BndSeg in FACE boundary
294     int              _branchID; // negative ID means reverse direction
295
296     BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
297       _inSeg(seg), _edge(edge), _uLast(u), _prev(0), _branchID( theNoBrachID ) {}
298
299     void setIndexToEdge( size_t id )
300     {
301       SMESH_MAT2d::Branch::setBndSegment( id, _edge );
302     }
303
304     int branchID() const { return Abs( _branchID ); }
305
306     size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
307
308     static BndSeg* getBndSegOfEdge( const TVDEdge*              edge,
309                                     vector< vector< BndSeg > >& bndSegsPerEdge )
310     {
311       BndSeg* seg = 0;
312       if ( edge )
313       {
314         size_t oppSegIndex  = SMESH_MAT2d::Branch::getBndSegment( edge );
315         size_t oppEdgeIndex = SMESH_MAT2d::Branch::getGeomEdge  ( edge );
316         if ( oppEdgeIndex < bndSegsPerEdge.size() &&
317              oppSegIndex < bndSegsPerEdge[ oppEdgeIndex ].size() )
318         {
319           seg = & bndSegsPerEdge[ oppEdgeIndex ][ oppSegIndex ];
320         }
321       }
322       return seg;
323     }
324
325     void setBranch( int branchID, vector< vector< BndSeg > >& bndSegsPerEdge )
326     {
327       _branchID = branchID;
328
329       // pass branch to an opposite BndSeg
330       if ( _edge )
331         if ( BndSeg* oppSeg = getBndSegOfEdge( _edge->twin(), bndSegsPerEdge ))
332         {
333           if ( oppSeg->_branchID == theNoBrachID )
334             oppSeg->_branchID = -branchID;
335         }
336     }
337     bool hasOppositeEdge()
338     {
339       if ( !_edge ) return false;
340       return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != theNoEdgeID );
341     }
342
343     // check a next segment in CCW order
344     bool isSameBranch( const BndSeg& seg2 )
345     {
346       if ( !_edge || !seg2._edge )
347         return true;
348
349       if ( _edge->twin() == seg2._edge )
350         return true;
351
352       const TVDCell* cell1 = this->_edge->twin()->cell();
353       const TVDCell* cell2 = seg2. _edge->twin()->cell();
354       if ( cell1 == cell2 )
355         return true;
356
357       const TVDEdge* edgeMedium1 = this->_edge->twin()->next();
358       const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev();
359
360       if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() )
361       {
362         if ( edgeMedium1->twin() == edgeMedium2 )
363           return true;
364         // edgeMedium's are edges whose twin()->cell is built on an end point of inSegment
365         // and is located between cell1 and cell2
366         if ( edgeMedium1->twin() == edgeMedium2->twin() ) // is this possible???
367           return true;
368         if ( edgeMedium1->twin() == edgeMedium2->twin()->next() &&
369              edgeMedium1->twin()->cell()->contains_point() )
370           return true;
371       }
372       else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
373       {
374         if ( edgeMedium1->twin() == edgeMedium2 &&
375              SMESH_MAT2d::Branch::getGeomEdge( edgeMedium1 ) ==
376              SMESH_MAT2d::Branch::getGeomEdge( edgeMedium2 ))
377           // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
378           return true;
379       }
380
381       return false;
382     }
383   }; // struct BndSeg
384
385   // -------------------------------------------------------------------------------------
386   /*!
387    * \brief Iterator 
388    */
389   struct BranchIterator
390   {
391     int                                 _i, _size;
392     const std::vector<const TVDEdge*> & _edges;
393     bool                                _closed;
394
395     BranchIterator(const std::vector<const TVDEdge*> & edges, int i )
396       :_i( i ), _size( edges.size() ), _edges( edges )
397     {
398       _closed = ( edges[0]->vertex1() == edges.back()->vertex0() || // closed branch
399                   edges[0]->vertex0() == edges.back()->vertex1() );
400     }
401     const TVDEdge* operator++() { ++_i; return edge(); }
402     const TVDEdge* operator--() { --_i; return edge(); }
403     bool operator<( const BranchIterator& other ) { return _i < other._i; }
404     BranchIterator& operator=( const BranchIterator& other ) { _i = other._i; return *this; }
405     void set(int i) { _i = i; }
406     int  index() const { return _i; }
407     int  indexMod() const { return ( _i + _size ) % _size; }
408     const TVDEdge* edge() const {
409       return _closed ? _edges[ indexMod() ] : ( _i < 0 || _i >= _size ) ? 0 : _edges[ _i ];
410     }
411     const TVDEdge* edgePrev() { --_i; const TVDEdge* e = edge(); ++_i; return e; }
412     const TVDEdge* edgeNext() { ++_i; const TVDEdge* e = edge(); --_i; return e; }
413   };
414
415   //================================================================================
416   /*!
417    * \brief debug: to visually check found MA edges
418    */
419   //================================================================================
420
421   void bndSegsToMesh( const vector< vector< BndSeg > >& bndSegsPerEdge )
422   {
423 #ifdef _MYDEBUG_
424     if ( !getenv("bndSegsToMesh")) return;
425     map< const TVDVertex *, int > v2Node;
426     map< const TVDVertex *, int >::iterator v2n;
427     set< const TVDEdge* > addedEdges;
428
429     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
430     SMESH_File file(fileName, false );
431     file.remove();
432     file.openForWriting();
433     SMESH_Comment text;
434     text << "import salome, SMESH\n";
435     text << "salome.salome_init()\n";
436     text << "from salome.smesh import smeshBuilder\n";
437     text << "smesh = smeshBuilder.New()\n";
438     text << "m=smesh.Mesh()\n";
439     for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
440     {
441       const vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
442       for ( size_t i = 0; i < bndSegs.size(); ++i )
443       {
444         if ( !bndSegs[i]._edge )
445           text << "# E=" << iE << " i=" << i << " NULL edge\n";
446         else if ( !bndSegs[i]._edge->vertex0() ||
447                   !bndSegs[i]._edge->vertex1() )
448           text << "# E=" << iE << " i=" << i << " INFINITE edge\n";
449         else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
450                   addedEdges.insert( bndSegs[i]._edge->twin() ).second )
451         {
452           v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
453           size_t n0 = v2n->second;
454           if ( n0 == v2Node.size() )
455             text << "n" << n0 << " = m.AddNode( "
456                  << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
457                  << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
458
459           v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
460           size_t n1 = v2n->second;
461           if ( n1 == v2Node.size() )
462             text << "n" << n1 << " = m.AddNode( "
463                  << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
464                  << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
465
466           text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
467         }
468       }
469     }
470     text << "\n";
471     file.write( text.c_str(), text.size() );
472     cout << fileName << endl;
473 #endif
474   }
475
476   //================================================================================
477   /*!
478    * \brief Computes length of a TVDEdge
479    */
480   //================================================================================
481
482   double length( const TVDEdge* edge )
483   {
484     gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(),
485              edge->vertex0()->y() - edge->vertex1()->y() );
486     return d.Modulus();
487   }
488
489   //================================================================================
490   /*!
491    * \brief Compute scale to have the same 2d proportions as in 3d
492    */
493   //================================================================================
494
495   void computeProportionScale( const TopoDS_Face& face,
496                                const Bnd_B2d&     uvBox,
497                                double             scale[2])
498   {
499     scale[0] = scale[1] = 1.;
500     if ( uvBox.IsVoid() ) return;
501
502     TopLoc_Location loc;
503     Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
504
505     const int nbDiv = 30;
506     gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
507     gp_XY uvMid = 0.5 * ( uvMin + uvMax );
508     double du = ( uvMax.X() - uvMin.X() ) / nbDiv;
509     double dv = ( uvMax.Y() - uvMin.Y() ) / nbDiv;
510
511     double uLen3d = 0, vLen3d = 0;
512     gp_Pnt uPrevP = surface->Value( uvMin.X(), uvMid.Y() );
513     gp_Pnt vPrevP = surface->Value( uvMid.X(), uvMin.Y() );
514     for (int i = 1; i <= nbDiv; i++)
515     {
516       double u = uvMin.X() + du * i;
517       double v = uvMin.Y() + dv * i;
518       gp_Pnt uP = surface->Value( u, uvMid.Y() );
519       gp_Pnt vP = surface->Value( uvMid.X(), v );
520       uLen3d += uP.Distance( uPrevP );
521       vLen3d += vP.Distance( vPrevP );
522       uPrevP = uP;
523       vPrevP = vP;
524     }
525     scale[0] = uLen3d / ( uvMax.X() - uvMin.X() );
526     scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() );
527   }
528
529   //================================================================================
530   /*!
531    * \brief Fill input data for construct_voronoi()
532    */
533   //================================================================================
534
535   bool makeInputData(const TopoDS_Face&                face,
536                      const std::vector< TopoDS_Edge >& edges,
537                      const double                      minSegLen,
538                      vector< InPoint >&                inPoints,
539                      vector< InSegment>&               inSegments,
540                      double                            scale[2])
541   {
542     const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization
543     TopLoc_Location loc;
544
545     // discretize the EDGEs to get 2d points and segments
546
547     vector< vector< UVU > > uvuVec( edges.size() );
548     Bnd_B2d uvBox;
549     for ( size_t iE = 0; iE < edges.size(); ++iE )
550     {
551       vector< UVU > & points = uvuVec[ iE ];
552
553       double f,l;
554       Handle(Geom_Curve)   c3d = BRep_Tool::Curve         ( edges[ iE ], loc, f, l );
555       Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface( edges[ iE ], face, f, l );
556       if ( c2d.IsNull() ) return false;
557
558       points.push_back( UVU( c2d->Value( f ), f ));
559       uvBox.Add( points.back()._uv );
560
561       Geom2dAdaptor_Curve c2dAdaptor (c2d, f,l );
562       double curDeflect = 0.3; //0.01;  //Curvature deflection
563       double angDeflect = 0.2; // 0.09; //Angular deflection
564
565       GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect);
566       // if ( discret.NbPoints() > 2 )
567       // {
568       //   cout << endl;
569       //   do
570       //   {
571       //     discret.Initialize( c2dAdaptor, 100, curDeflect );
572       //     cout << "C " << curDeflect << " " << discret.NbPoints() << endl;
573       //     curDeflect *= 1.5;
574       //   }
575       //   while ( discret.NbPoints() > 5 );
576       //   cout << endl;
577       //   do
578       //   {
579       //     discret.Initialize( c2dAdaptor, angDeflect, 100 );
580       //     cout << "A " << angDeflect << " " << discret.NbPoints() << endl;
581       //     angDeflect *= 1.5;
582       //   }
583       //   while ( discret.NbPoints() > 5 );
584       // }
585       gp_Pnt p, pPrev;
586       if ( !c3d.IsNull() )
587         pPrev = c3d->Value( f );
588       if ( discret.NbPoints() > 2 )
589         for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
590         {
591           double u = discret.Parameter(i);
592           if ( !c3d.IsNull() )
593           {
594             p = c3d->Value( u );
595             int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
596             double dU = ( u - points.back()._u ) / nbDiv;
597             for ( int iD = 1; iD < nbDiv; ++iD )
598             {
599               double uD = points.back()._u + dU;
600               points.push_back( UVU( c2d->Value( uD ), uD ));
601             }
602             pPrev = p;
603           }
604           points.push_back( UVU( c2d->Value( u ), u ));
605           uvBox.Add( points.back()._uv );
606         }
607       // if ( !c3d.IsNull() )
608       // {
609       //   vector<double> params;
610       //   GeomAdaptor_Curve c3dAdaptor( c3d,f,l );
611       //   if ( useDefl )
612       //   {
613       //     const double deflection = minSegLen * 0.1;
614       //     GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true );
615       //     if ( !discret.IsDone() )
616       //       return false;
617       //     int nbP = discret.NbPoints();
618       //     for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
619       //       params.push_back( discret.Parameter(i) );
620       //   }
621       //   else
622       //   {
623       //     double   eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor );
624       //     int     nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef ));
625       //     double segLen = eLen / nbSeg;
626       //     GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l );
627       //     int nbP = Min( discret.NbPoints(), nbSeg + 1 );
628       //     for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
629       //       params.push_back( discret.Parameter(i) );
630       //   }
631       //   for ( size_t i = 0; i < params.size(); ++i )
632       //   {
633       //     points.push_back( UVU( c2d->Value( params[i] ), params[i] ));
634       //     uvBox.Add( points.back()._uv );
635       //   }
636       // }
637       if ( points.size() < 2 )
638       {
639         points.push_back( UVU( c2d->Value( l ), l ));
640         uvBox.Add( points.back()._uv );
641       }
642       if ( edges[ iE ].Orientation() == TopAbs_REVERSED )
643         std::reverse( points.begin(), points.end() );
644     }
645
646     // make connected EDGEs have same UV at shared VERTEX
647     TopoDS_Vertex vShared;
648     for ( size_t iE = 0; iE < edges.size(); ++iE )
649     {
650       size_t iE2 = (iE+1) % edges.size();
651       if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared )) // FACE with several WIREs?
652         for ( size_t i = 1; i < edges.size(); ++i )
653         {
654           iE2 = (iE2+1) % edges.size();
655           if ( iE != iE2 &&
656                TopExp::CommonVertex( edges[iE], edges[iE2], vShared ) &&
657                vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
658             break;
659         }
660       if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
661         continue;
662         //return false;
663       vector< UVU > & points1 = uvuVec[ iE ];
664       vector< UVU > & points2 = uvuVec[ iE2 ];
665       gp_Pnt2d & uv1 = points1.back() ._uv;
666       gp_Pnt2d & uv2 = points2.front()._uv;
667       uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() );
668     }
669
670     // get scale to have the same 2d proportions as in 3d
671     computeProportionScale( face, uvBox, scale );
672
673     // make 'scale' such that to have coordinates precise enough when converted to int
674
675     gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
676     uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
677     uvMin.ChangeCoord(2) = uvMin.Y() * scale[1];
678     uvMax.ChangeCoord(1) = uvMax.X() * scale[0];
679     uvMax.ChangeCoord(2) = uvMax.Y() * scale[1];
680     double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
681                        Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
682     int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
683     const double precision = Min( 1e-5, minSegLen * 1e-2 );
684     double preciScale = Min( vMax[iMax] / precision,
685                              std::numeric_limits<int>::max() / vMax[iMax] );
686     preciScale /= scale[iMax];
687     double roundedScale = 10; // to ease debug
688     while ( roundedScale * 10 < preciScale )
689       roundedScale *= 10.;
690     scale[0] *= roundedScale;
691     scale[1] *= roundedScale;
692
693     // create input points and segments
694
695     inPoints.clear();
696     inSegments.clear();
697     size_t nbPnt = 0;
698     for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
699       nbPnt += uvuVec[ iE ].size();
700     inPoints.resize( nbPnt );
701     inSegments.reserve( nbPnt );
702
703     size_t iP = 0;
704     if ( face.Orientation() == TopAbs_REVERSED )
705     {
706       for ( int iE = uvuVec.size()-1; iE >= 0; --iE )
707       {
708         vector< UVU > & points = uvuVec[ iE ];
709         inPoints[ iP++ ] = points.back().getInPoint( scale );
710         for ( size_t i = points.size()-1; i >= 1; --i )
711         {
712           inPoints[ iP++ ] = points[i-1].getInPoint( scale );
713           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
714           if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
715             return false; // too short segment
716         }
717       }
718     }
719     else
720     {
721       for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
722       {
723         vector< UVU > & points = uvuVec[ iE ];
724         inPoints[ iP++ ] = points[0].getInPoint( scale );
725         for ( size_t i = 1; i < points.size(); ++i )
726         {
727           inPoints[ iP++ ] = points[i].getInPoint( scale );
728           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
729           if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
730             return false; // too short segment
731         }
732       }
733     }
734     // debug
735     theScale[0] = scale[0];
736     theScale[1] = scale[1];
737
738     return true;
739   }
740
741   //================================================================================
742   /*!
743    * \brief Update a branch joined to another one
744    */
745   //================================================================================
746
747   void updateJoinedBranch( vector< const TVDEdge* > &   branchEdges,
748                            const size_t                 newID,
749                            vector< vector< BndSeg > > & bndSegs,
750                            const bool                   reverse)
751   {
752     BndSeg *seg1, *seg2;
753     if ( reverse )
754     {
755       for ( size_t i = 0; i < branchEdges.size(); ++i )
756       {
757         if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i],         bndSegs )) &&
758             ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
759         {
760           seg1->_branchID /= seg1->branchID();
761           seg2->_branchID /= seg2->branchID();
762           seg1->_branchID *= -newID;
763           seg2->_branchID *= -newID;
764           branchEdges[i] = branchEdges[i]->twin();
765         }
766       }
767       std::reverse( branchEdges.begin(), branchEdges.end() );
768     }
769     else
770     {
771       for ( size_t i = 0; i < branchEdges.size(); ++i )
772       {
773         if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i],         bndSegs )) &&
774             ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
775         {
776           seg1->_branchID /= seg1->branchID();
777           seg2->_branchID /= seg2->branchID();
778           seg1->_branchID *= newID;
779           seg2->_branchID *= newID;
780         }
781       }
782     }
783   }
784
785   //================================================================================
786   /*!
787    * \brief Create MA branches and FACE boundary data
788    *  \param [in] vd - voronoi diagram of \a inSegments
789    *  \param [in] inPoints - FACE boundary points
790    *  \param [in,out] inSegments - FACE boundary segments
791    *  \param [out] branch - MA branches to fill
792    *  \param [out] branchEnd - ends of MA branches to fill
793    *  \param [out] boundary - FACE boundary to fill
794    */
795   //================================================================================
796
797   void makeMA( const TVD&                               vd,
798                const bool                               ignoreCorners,
799                vector< InPoint >&                       inPoints,
800                vector< InSegment > &                    inSegments,
801                vector< SMESH_MAT2d::Branch >&           branch,
802                vector< const SMESH_MAT2d::BranchEnd* >& branchPnt,
803                SMESH_MAT2d::Boundary&                   boundary )
804   {
805     // Associate MA cells with geom EDGEs
806     for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
807     {
808       const TVDCell* cell = &(*it);
809       if ( cell->is_degenerate() )
810       {
811         std::cerr << "SMESH_MAT2d: encounter degenerate voronoi_cell. Invalid input data?"
812                   << std::endl;
813         return;
814       }
815       if ( cell->contains_segment() )
816       {
817         InSegment& seg = inSegments[ cell->source_index() ];
818         seg._cell = cell;
819         seg.setGeomEdgeToCell( cell, seg._geomEdgeInd );
820       }
821       else
822       {
823         InSegment::setGeomEdgeToCell( cell, theNoEdgeID );
824       }
825     }
826
827     vector< bool > inPntChecked( inPoints.size(), false );
828
829     // Find MA edges of each inSegment
830
831     for ( size_t i = 0; i < inSegments.size(); ++i )
832     {
833       InSegment& inSeg = inSegments[i];
834
835       // get edges around the cell lying on MA
836       bool hasSecondary = false;
837       const TVDEdge* edge = inSeg._cell->incident_edge();
838       do {
839         edge = edge->next(); // Returns the CCW next edge within the cell.
840         if ( edge->is_primary() && !inSeg.isExternal( edge ) )
841           inSeg._edges.push_back( edge ); // edge equidistant from two InSegments
842         else
843           hasSecondary = true;
844       } while (edge != inSeg._cell->incident_edge());
845
846       // there can be several continuous MA edges but maEdges can begin in the middle of
847       // a chain of continuous MA edges. Make the chain continuous.
848       list< const TVDEdge* >& maEdges = inSeg._edges;
849       if ( maEdges.empty() )
850         continue;
851       if ( hasSecondary )
852         while ( maEdges.back()->next() == maEdges.front() )
853           maEdges.splice( maEdges.end(), maEdges, maEdges.begin() );
854
855       // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE
856       list< const TVDEdge* >::iterator e = maEdges.begin();
857       while ( e != maEdges.end() )
858       {
859         const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge
860         size_t         geoE2 = InSegment::getGeomEdge( cell2 );
861         bool        toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e ));
862         if ( toRemove )
863           e = maEdges.erase( e );
864         else
865           ++e;
866       }
867       if ( maEdges.empty() )
868         continue;
869
870       // add MA edges corresponding to concave InPoints
871       for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg
872       {
873         InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 );
874         size_t    pInd = inPnt.index( inPoints );
875         if ( inPntChecked[ pInd ] )
876           continue;
877         if ( pInd > 0 &&
878              inPntChecked[ pInd-1 ] &&
879              inPoints[ pInd-1 ] == inPnt )
880           continue;
881         inPntChecked[ pInd ] = true;
882
883         const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
884         if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
885           continue;
886         const TVDEdge* edge =  // a secondary TVDEdge connecting inPnt and maE
887           is2nd ? maE->prev() : maE->next();
888         while ( inSeg.isConnected( edge ))
889         {
890           if ( edge->is_primary() ) break; // this should not happen
891           const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
892           if ( inSeg.getGeomEdge( edge2->cell() ) != theNoEdgeID )
893             break; // cell of an InSegment
894           bool hasInfinite = false;
895           list< const TVDEdge* > pointEdges;
896           edge = edge2;
897           do
898           {
899             edge = edge->next(); // Returns the CCW next edge within the cell.
900             if ( edge->is_infinite() )
901               hasInfinite = true;
902             else if ( edge->is_primary() && !inSeg.isExternal( edge ))
903               pointEdges.push_back( edge );
904           }
905           while ( edge != edge2 && !hasInfinite );
906
907           if ( hasInfinite || pointEdges.empty() )
908             break;
909           inPnt._edges.splice( inPnt._edges.end(), pointEdges );
910           inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd );
911
912           edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next();
913         }
914       } // add MA edges corresponding to concave InPoints
915
916     } // loop on inSegments to find corresponding MA edges
917
918
919     // -------------------------------------------
920     // Create Branches and BndPoints for each EDGE
921     // -------------------------------------------
922
923     if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ )
924     {
925       inPntChecked[0] = false; // do not use the 1st point twice
926       //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), theNoEdgeID );
927       inPoints[0]._edges.clear();
928     }
929
930     // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
931
932     vector< vector< BndSeg > > bndSegsPerEdge( boundary.nbEdges() ); // all BndSeg's
933     {
934       vector< BndSeg > bndSegs; // bndSeg's of a current EDGE
935       size_t prevGeomEdge = theNoEdgeID;
936
937       list< const TVDEdge* >::reverse_iterator e;
938       for ( size_t i = 0; i < inSegments.size(); ++i )
939       {
940         InSegment& inSeg = inSegments[i];
941
942         if ( inSeg._geomEdgeInd != prevGeomEdge )
943         {
944           if ( !bndSegs.empty() )
945             bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
946           prevGeomEdge = inSeg._geomEdgeInd;
947         }
948
949         // segments around 1st concave point
950         size_t ip0 = inSeg._p0->index( inPoints );
951         if ( inPntChecked[ ip0 ] )
952           for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
953             bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
954         inPntChecked[ ip0 ] = false;
955
956         // segments of InSegment's
957         const size_t nbMaEdges = inSeg._edges.size();
958         switch ( nbMaEdges ) {
959         case 0: // "around" circle center
960           bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
961         case 1:
962           bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
963         default:
964           gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
965                           inSeg._p1->_b - inSeg._p0->_b );
966           const double inSegLen2 = inSegDir.SquareModulus();
967           e = inSeg._edges.rbegin();
968           for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
969           {
970             gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
971                         (*e)->vertex0()->y() - inSeg._p0->_b );
972             double r = toMA * inSegDir / inSegLen2;
973             double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
974             bndSegs.push_back( BndSeg( &inSeg, *e, u ));
975           }
976           bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
977         }
978         // segments around 2nd concave point
979         size_t ip1 = inSeg._p1->index( inPoints );
980         if ( inPntChecked[ ip1 ] )
981           for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
982             bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
983         inPntChecked[ ip1 ] = false;
984       }
985       if ( !bndSegs.empty() )
986         bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
987     }
988
989     // prepare to MA branch search
990     for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
991     {
992       // 1) make TVDEdge's know it's BndSeg to enable passing branchID to
993       // an opposite BndSeg in BndSeg::setBranch(); geom EDGE ID is known from TVDCell
994       // 2) connect bndSegs via BndSeg::_prev
995
996       vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
997       if ( bndSegs.empty() ) continue;
998
999       for ( size_t i = 1; i < bndSegs.size(); ++i )
1000       {
1001         bndSegs[i]._prev = & bndSegs[i-1];
1002         bndSegs[i].setIndexToEdge( i );
1003       }
1004       // look for the last bndSeg of previous EDGE to set bndSegs[0]._prev
1005       const InPoint& p0 = bndSegs[0]._inSeg->point0();
1006       for ( size_t iE2 = 0; iE2 < bndSegsPerEdge.size(); ++iE2 )
1007         if ( p0 == bndSegsPerEdge[ iE2 ].back()._inSeg->point1() )
1008         {
1009           bndSegs[0]._prev = & bndSegsPerEdge[ iE2 ].back();
1010           break;
1011         }
1012       bndSegs[0].setIndexToEdge( 0 );
1013     }
1014
1015     bndSegsToMesh( bndSegsPerEdge ); // debug: visually check found MA edges
1016
1017
1018     // Find TVDEdge's of Branches and associate them with bndSegs
1019
1020     vector< vector<const TVDEdge*> > branchEdges;
1021     branchEdges.reserve( boundary.nbEdges() * 4 );
1022
1023     map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
1024
1025     int branchID = 1; // we code orientation as branchID sign
1026     branchEdges.resize( branchID );
1027
1028     vector< std::pair< int, const TVDVertex* > > branchesToCheckEnd;
1029
1030     for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
1031     {
1032       vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
1033       for ( size_t i = 0; i < bndSegs.size(); ++i )
1034       {
1035         if ( bndSegs[i].branchID() )
1036         {
1037           if ( bndSegs[i]._prev &&
1038                bndSegs[i]._branchID == -bndSegs[i]._prev->_branchID &&
1039                bndSegs[i]._edge )
1040           {
1041             SMESH_MAT2d::BranchEndType type =
1042               ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
1043                 SMESH_MAT2d::BE_ON_VERTEX :
1044                 SMESH_MAT2d::BE_END );
1045             endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
1046           }
1047           continue;
1048         }
1049         if ( !bndSegs[i]._prev &&
1050              !bndSegs[i].hasOppositeEdge() )
1051           continue;
1052
1053         if ( !bndSegs[i]._prev ||
1054              !bndSegs[i]._prev->isSameBranch( bndSegs[i] ))
1055         {
1056           branchEdges.resize(( branchID = branchEdges.size()) + 1 );
1057           if ( bndSegs[i]._edge && bndSegs[i]._prev )
1058           {
1059             endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_BRANCH_POINT ));
1060             if ( bndSegs[i]._prev->_branchID < 0 )
1061               // 0023404: a branch-point is inside a branch
1062               branchesToCheckEnd.push_back( make_pair( bndSegs[i]._prev->branchID(),
1063                                                        bndSegs[i]._edge->vertex1() ));
1064           }
1065         }
1066         else if ( bndSegs[i]._prev->_branchID )
1067         {
1068           branchID = bndSegs[i]._prev->_branchID;  // with sign
1069         }
1070         else if ( bndSegs[i]._edge ) // 1st bndSeg of a WIRE
1071         {
1072           branchEdges.resize(( branchID = branchEdges.size()) + 1 );
1073           if ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ))
1074           {
1075             if ( bndSegs[i]._inSeg->point0() == bndSegs[i]._edge->vertex1() )
1076               endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_ON_VERTEX ));
1077             else
1078               endType.insert( make_pair( bndSegs[i]._edge->vertex0(), SMESH_MAT2d::BE_ON_VERTEX ));
1079           }
1080         }
1081
1082         bndSegs[i].setBranch( branchID, bndSegsPerEdge ); // set to i-th and to the opposite bndSeg
1083         if ( bndSegs[i].hasOppositeEdge() )
1084           branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
1085       }
1086     }
1087
1088     if ( !ignoreCorners && !branchesToCheckEnd.empty() )
1089     {
1090       // split branches having branch-point inside
1091       // (a branch-point was not detected since another branch is joined at the opposite side)
1092       for ( size_t i = 0; i < branchesToCheckEnd.size(); ++i )
1093       {
1094         vector<const TVDEdge*> & branch = branchEdges[ branchesToCheckEnd[i].first ];
1095         const TVDVertex*    branchPoint = branchesToCheckEnd[i].second;
1096         if ( branch.front()->vertex1() == branchPoint ||
1097              branch.back ()->vertex0() == branchPoint )
1098           continue; // OK - branchPoint is at a branch end
1099
1100         // find a MA edge where another branch begins
1101         size_t iE;
1102         for ( iE = 0; iE < branch.size(); ++iE )
1103           if ( branch[iE]->vertex1() == branchPoint )
1104             break;
1105         if ( iE < branch.size() )
1106         {
1107           // split the branch
1108           branchEdges.resize(( branchID = branchEdges.size()) + 1 );
1109           vector<const TVDEdge*> & branch2 = branchEdges[ branchID ];
1110           branch2.assign( branch.begin()+iE, branch.end() );
1111           branch.resize( iE );
1112           for ( iE = 0; iE < branch2.size(); ++iE )
1113             if ( BndSeg* bs = BndSeg::getBndSegOfEdge( branch2[iE], bndSegsPerEdge ))
1114               bs->setBranch( branchID, bndSegsPerEdge );
1115         }
1116       }
1117     }
1118
1119     // join the 1st and the last branch edges if it is the same branch
1120     // if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
1121     //      bndSegs.back().isSameBranch( bndSegs.front() ))
1122     // {
1123     //   vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
1124     //   vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID()  ];
1125     //   br1.insert( br1.begin(), br2.begin(), br2.end() );
1126     //   br2.clear();
1127     // }
1128
1129     // remove branches ending at BE_ON_VERTEX and BE_END
1130
1131     vector<bool> isBranchRemoved( branchEdges.size(), false );
1132
1133     std::set< SMESH_MAT2d::BranchEndType > endTypeToRm;
1134     endTypeToRm.insert( SMESH_MAT2d::BE_ON_VERTEX );
1135     endTypeToRm.insert( SMESH_MAT2d::BE_END );
1136
1137     if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
1138     {
1139       // find branches to remove
1140       map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et;
1141       for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
1142       {
1143         if ( branchEdges[iB].empty() )
1144           continue;
1145         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
1146         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
1147         v2et = endType.find( v0 );
1148         if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
1149           isBranchRemoved[ iB ] = true;
1150         v2et = endType.find( v1 );
1151         if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
1152           isBranchRemoved[ iB ] = true;
1153       }
1154       // try to join not removed branches into one
1155       for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
1156       {
1157         if ( branchEdges[iB].empty() || isBranchRemoved[iB] )
1158           continue;
1159         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
1160         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
1161         v2et = endType.find( v0 );
1162         if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
1163           v0 = 0;
1164         v2et = endType.find( v1 );
1165         if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
1166           v1 = 0;
1167         if ( !v0 && !v1 )
1168           continue;
1169
1170         for ( int isV0 = 0; isV0 < 2; ++isV0 )
1171         {
1172           const TVDVertex* v = isV0 ? v0 : v1;
1173           size_t iBrToJoin = 0;
1174           for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
1175           {
1176             if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
1177               continue;
1178             const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
1179             const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
1180             if ( v == v02 || v == v12 )
1181             {
1182               if ( iBrToJoin > 0 )
1183               {
1184                 iBrToJoin = 0;
1185                 break; // more than 2 not removed branches meat at a TVDVertex
1186               }
1187               iBrToJoin = iB2;
1188             }
1189           }
1190           if ( iBrToJoin > 0 )
1191           {
1192             vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
1193             const TVDVertex* v02 = branch[0]->vertex1();
1194             const TVDVertex* v12 = branch.back()->vertex0();
1195             updateJoinedBranch( branch, iB, bndSegsPerEdge, /*reverse=*/(v0 == v02 || v1 == v12 ));
1196             if ( v0 == v02 || v0 == v12 )
1197               branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
1198             else
1199               branchEdges[iB].insert( branchEdges[iB].end(),   branch.begin(), branch.end() );
1200             branch.clear();
1201           }
1202         }
1203       } // loop on branchEdges
1204     } // if ( ignoreCorners )
1205
1206     // associate branchIDs and the input branch vector (arg)
1207     vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 );
1208     int nbBranches = 0;
1209     for ( size_t i = 0; i < branchEdges.size(); ++i )
1210     {
1211       nbBranches += ( !branchEdges[i].empty() );
1212     }
1213     branch.resize( nbBranches );
1214     size_t iBr = 0;
1215     for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed
1216     {
1217       if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ])
1218         branchByID[ brID ] = & branch[ iBr++ ];
1219     }
1220     for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed
1221     {
1222       if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ])
1223         branchByID[ brID ] = & branch[ iBr++ ];
1224     }
1225
1226     // Fill in BndPoints of each EDGE of the boundary
1227
1228     //size_t iSeg = 0;
1229     int edgeInd = -1, dInd = 0;
1230     for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
1231     {
1232       vector< BndSeg >&          bndSegs = bndSegsPerEdge[ iE ];
1233       SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( iE );
1234
1235       // make TVDEdge know an index of bndSegs within BndPoints
1236       for ( size_t i = 0; i < bndSegs.size(); ++i )
1237         if ( bndSegs[i]._edge )
1238           SMESH_MAT2d::Branch::setBndSegment( i, bndSegs[i]._edge );
1239
1240       // parameters on EDGE
1241
1242       bndPoints._params.reserve( bndSegs.size() + 1 );
1243       bndPoints._params.push_back( bndSegs[ 0 ]._inSeg->_p0->_param );
1244
1245       for ( size_t i = 0; i < bndSegs.size(); ++i )
1246         bndPoints._params.push_back( bndSegs[ i ]._uLast );
1247
1248       // MA edges
1249
1250       bndPoints._maEdges.reserve( bndSegs.size() );
1251
1252       for ( size_t i = 0; i < bndSegs.size(); ++i )
1253       {
1254         const size_t              brID = bndSegs[ i ].branchID();
1255         const SMESH_MAT2d::Branch*  br = branchByID[ brID ];
1256
1257         if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
1258         {
1259           edgeInd += dInd;
1260
1261           if (( edgeInd < 0 ||
1262                 edgeInd >= (int) branchEdges[ brID ].size() ) ||
1263               ( branchEdges[ brID ][ edgeInd ]         != bndSegs[ i ]._edge &&
1264                 branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge ))
1265           {
1266             if ( bndSegs[ i ]._branchID < 0 )
1267             {
1268               dInd = -1;
1269               for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd )
1270                 if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge )
1271                   break;
1272             }
1273             else // bndSegs[ i ]._branchID > 0
1274             {
1275               dInd = +1;
1276               for ( edgeInd = 0; edgeInd < (int)branchEdges[ brID ].size(); ++edgeInd )
1277                 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
1278                   break;
1279             }
1280           }
1281         }
1282         else
1283         {
1284           // no MA edge, bndSeg corresponds to an end point of a branch
1285           if ( bndPoints._maEdges.empty() )
1286             edgeInd = 0;
1287           else
1288             edgeInd = branchEdges[ brID ].size();
1289           dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
1290         }
1291
1292         bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
1293
1294       } // loop on bndSegs of an EDGE
1295     } // loop on all bndSegs to construct Boundary
1296
1297     // Initialize branches
1298
1299     // find a not removed branch
1300     size_t iBrNorRemoved = 0;
1301     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1302       if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] )
1303       {
1304         iBrNorRemoved = brID;
1305         break;
1306       }
1307     // fill the branches with MA edges
1308     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1309       if ( !branchEdges[brID].empty() )
1310       {
1311         branchByID[ brID ]->init( branchEdges[brID], & boundary, endType );
1312       }
1313     // mark removed branches
1314     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1315       if ( isBranchRemoved[brID] && iBrNorRemoved > 0 )
1316       {
1317         SMESH_MAT2d::Branch* branch     = branchByID[ brID ];
1318         SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ];
1319         bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT );
1320         const TVDVertex* branchVextex = 
1321           is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex;
1322         SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex );
1323         branch->setRemoved( bp );
1324       }
1325     // set branches to branch ends
1326     for ( size_t i = 0; i < branch.size(); ++i )
1327       if ( !branch[i].isRemoved() )
1328         branch[i].setBranchesToEnds( branch );
1329
1330     // fill branchPnt arg
1331     map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
1332     for ( size_t i = 0; i < branch.size(); ++i )
1333     {
1334       if ( branch[i].getEnd(0)->_branches.size() > 2 )
1335         v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
1336       if ( branch[i].getEnd(1)->_branches.size() > 2 )
1337         v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
1338     }
1339     branchPnt.resize( v2end.size() );
1340     map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
1341     for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
1342       branchPnt[ i ] = v2e->second;
1343
1344   } // makeMA()
1345
1346 } // namespace
1347
1348 //================================================================================
1349 /*!
1350  * \brief MedialAxis constructor
1351  *  \param [in] face - a face to create MA for
1352  *  \param [in] edges - edges of the face (possibly not all) on the order they
1353  *              encounter in the face boundary.
1354  *  \param [in] minSegLen - minimal length of a mesh segment used to discretize
1355  *              the edges. It is used to define precision of MA approximation
1356  */
1357 //================================================================================
1358
1359 SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face&                face,
1360                                     const std::vector< TopoDS_Edge >& edges,
1361                                     const double                      minSegLen,
1362                                     const bool                        ignoreCorners):
1363   _face( face ), _boundary( edges.size() )
1364 {
1365   // input to construct_voronoi()
1366   vector< InPoint >  inPoints;
1367   vector< InSegment> inSegments;
1368   if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
1369     return;
1370
1371   inSegmentsToFile( inSegments );
1372
1373   // build voronoi diagram
1374   construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
1375
1376   // make MA data
1377   makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary );
1378
1379   // count valid branches
1380   _nbBranches = _branch.size();
1381   for ( size_t i = 0; i < _branch.size(); ++i )
1382     if ( _branch[i].isRemoved() )
1383       --_nbBranches;
1384 }
1385
1386 //================================================================================
1387 /*!
1388  * \brief Returns the i-th branch
1389  */
1390 //================================================================================
1391
1392 const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const
1393 {
1394   return i < _nbBranches ? &_branch[i] : 0;
1395 }
1396
1397 //================================================================================
1398 /*!
1399  * \brief Return UVs of ends of MA edges of a branch
1400  */
1401 //================================================================================
1402
1403 void SMESH_MAT2d::MedialAxis::getPoints( const Branch*         branch,
1404                                          std::vector< gp_XY >& points) const
1405 {
1406   branch->getPoints( points, _scale );
1407 }
1408
1409 //================================================================================
1410 /*!
1411  * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
1412  *  \param [in] iEdge - index of geom EDGE within a vector passed at MA construction
1413  *  \param [in] u - parameter of the point on EDGE curve
1414  *  \param [out] p - the found BranchPoint
1415  *  \return bool - is OK
1416  */
1417 //================================================================================
1418
1419 bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
1420                                             double            u,
1421                                             BranchPoint&      p ) const
1422 {
1423   if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1424     return false;
1425
1426   const BndPoints& points = _pointsPerEdge[ iEdge ];
1427   const bool  edgeReverse = ( points._params[0] > points._params.back() );
1428
1429   if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
1430     u = edgeReverse ? points._params.back() : points._params[0];
1431   else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
1432     u = edgeReverse ? points._params[0] : points._params.back();
1433
1434   double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
1435   int    i = int( r * double( points._maEdges.size()-1 ));
1436   if ( edgeReverse )
1437   {
1438     while ( points._params[i  ] < u ) --i;
1439     while ( points._params[i+1] > u ) ++i;
1440   }
1441   else
1442   {
1443     while ( points._params[i  ] > u ) --i;
1444     while ( points._params[i+1] < u ) ++i;
1445   }
1446
1447   if ( points._params[i] == points._params[i+1] ) // coincident points at some end
1448   {
1449     int di = ( points._params[0] == points._params[i] ) ? +1 : -1;
1450     while ( points._params[i] == points._params[i+1] )
1451       i += di;
1452     if ( i < 0 || i+1 >= (int)points._params.size() )
1453       i = 0;
1454   }
1455
1456   double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
1457
1458   if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
1459   {
1460     if ( i < (int)points._maEdges.size() / 2 ) // near 1st point
1461     {
1462       while ( i < (int)points._maEdges.size()-1 && !points._maEdges[ i ].second )
1463         ++i;
1464       edgeParam = edgeReverse;
1465     }
1466     else // near last point
1467     {
1468       while ( i > 0 && !points._maEdges[ i ].second )
1469         --i;
1470       edgeParam = !edgeReverse;
1471     }
1472   }
1473   const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
1474   bool maReverse = ( maE.second < 0 );
1475
1476   p._branch    = maE.first;
1477   p._iEdge     = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
1478   p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam;
1479
1480   return true;
1481 }
1482
1483 //================================================================================
1484 /*!
1485  * \brief Returns a BranchPoint corresponding to a given BoundaryPoint on a geom EDGE
1486  *  \param [in] bp - the BoundaryPoint
1487  *  \param [out] p - the found BranchPoint
1488  *  \return bool - is OK
1489  */
1490 //================================================================================
1491
1492 bool SMESH_MAT2d::Boundary::getBranchPoint( const BoundaryPoint& bp,
1493                                             BranchPoint&         p ) const
1494 {
1495   return getBranchPoint( bp._edgeIndex, bp._param, p );
1496 }
1497
1498 //================================================================================
1499 /*!
1500  * \brief Check if a given boundary segment is a null-length segment on a concave
1501  *        boundary corner.
1502  *  \param [in] iEdge - index of a geom EDGE
1503  *  \param [in] iSeg - index of a boundary segment
1504  *  \return bool - true if the segment is on concave corner
1505  */
1506 //================================================================================
1507
1508 bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
1509 {
1510   if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1511     return false;
1512
1513   const BndPoints& points = _pointsPerEdge[ iEdge ];
1514   if ( points._params.size() <= iSeg+1 )
1515     return false;
1516
1517   return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20;
1518 }
1519
1520 //================================================================================
1521 /*!
1522  * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end
1523  */
1524 //================================================================================
1525
1526 bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
1527 {
1528   if ( bp._edgeIndex >= _pointsPerEdge.size() )
1529     return false;
1530
1531   const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
1532   if ( Abs( bp._param - points._params[0]) < Abs( points._params.back() - bp._param ))
1533     bp._param = points._params[0];
1534   else 
1535     bp._param = points._params.back();
1536
1537   return true;
1538 }
1539
1540 //================================================================================
1541 /*!
1542  * \brief Creates a 3d curve corresponding to a Branch
1543  *  \param [in] branch - the Branch
1544  *  \return Adaptor3d_Curve* - the new curve the caller is to delete
1545  */
1546 //================================================================================
1547
1548 Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
1549 {
1550   Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
1551   if ( surface.IsNull() )
1552     return 0;
1553
1554   vector< gp_XY > uv;
1555   branch.getPoints( uv, _scale );
1556   if ( uv.size() < 2 )
1557     return 0;
1558
1559   vector< TopoDS_Vertex > vertex( uv.size() );
1560   for ( size_t i = 0; i < uv.size(); ++i )
1561     vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
1562
1563   TopoDS_Wire aWire;
1564   BRep_Builder aBuilder;
1565   aBuilder.MakeWire(aWire);
1566   for ( size_t i = 1; i < vertex.size(); ++i )
1567   {
1568     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
1569     aBuilder.Add( aWire, edge );
1570   }
1571
1572   // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
1573   //   aWire.Closed(true); // issue 0021141
1574
1575   return new BRepAdaptor_CompCurve( aWire );
1576 }
1577
1578 //================================================================================
1579 /*!
1580  * \brief Copy points of an EDGE
1581  */
1582 //================================================================================
1583
1584 void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                 maEdges,
1585                                 const Boundary*                         boundary,
1586                                 map< const TVDVertex*, BranchEndType >& endType )
1587 {
1588   if ( maEdges.empty() ) return;
1589
1590   _boundary = boundary;
1591   _maEdges.swap( maEdges );
1592
1593
1594   _params.reserve( _maEdges.size() + 1 );
1595   _params.push_back( 0. );
1596   for ( size_t i = 0; i < _maEdges.size(); ++i )
1597     _params.push_back( _params.back() + length( _maEdges[i] ));
1598   
1599   for ( size_t i = 1; i < _params.size(); ++i )
1600     _params[i] /= _params.back();
1601
1602
1603   _endPoint1._vertex = _maEdges.front()->vertex1();
1604   _endPoint2._vertex = _maEdges.back ()->vertex0();
1605
1606   if ( endType.count( _endPoint1._vertex ))
1607     _endPoint1._type = endType[ _endPoint1._vertex ];
1608   if ( endType.count( _endPoint2._vertex ))
1609     _endPoint2._type = endType[ _endPoint2._vertex ];
1610 }
1611
1612 //================================================================================
1613 /*!
1614  * \brief fills BranchEnd::_branches of its ends
1615  */
1616 //================================================================================
1617
1618 void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
1619 {
1620   for ( size_t i = 0; i < branches.size(); ++i )
1621   {
1622     if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
1623          this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
1624       this->_endPoint1._branches.push_back( &branches[i] );
1625
1626     if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
1627          this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
1628       this->_endPoint2._branches.push_back( &branches[i] );
1629   }
1630 }
1631
1632 //================================================================================
1633 /*!
1634  * \brief returns a BranchPoint corresponding to a TVDVertex
1635  */
1636 //================================================================================
1637
1638 SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const
1639 {
1640   BranchPoint p;
1641   p._branch = this;
1642   p._iEdge  = 0;
1643
1644   if ( vertex == _maEdges[0]->vertex1() )
1645   {
1646     p._edgeParam = 0;
1647   }
1648   else
1649   {
1650     for ( ; p._iEdge < _maEdges.size(); ++p._iEdge )
1651       if ( vertex == _maEdges[ p._iEdge ]->vertex0() )
1652       {
1653         p._edgeParam = _params[ p._iEdge ];
1654         break;
1655       }
1656   }
1657   return p;
1658 }
1659
1660 //================================================================================
1661 /*!
1662  * \brief Sets a proxy point for a removed branch
1663  *  \param [in] proxyPoint - a point of another branch to which all points of this
1664  *         branch are mapped
1665  */
1666 //================================================================================
1667
1668 void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint )
1669 {
1670   _proxyPoint = proxyPoint;
1671 }
1672
1673 //================================================================================
1674 /*!
1675  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1676  *  \param [in] param - [0;1] normalized param on the Branch
1677  *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1678  *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1679  *  \return bool - true if the BoundaryPoint's found
1680  */
1681 //================================================================================
1682
1683 bool SMESH_MAT2d::Branch::getBoundaryPoints(double         param,
1684                                             BoundaryPoint& bp1,
1685                                             BoundaryPoint& bp2 ) const
1686 {
1687   if ( param < _params[0] || param > _params.back() )
1688     return false;
1689
1690   // look for an index of a MA edge by param
1691   double ip = param * _params.size();
1692   size_t  i = size_t( Min( int( _maEdges.size()-1), int( ip )));
1693
1694   while ( param < _params[i  ] ) --i;
1695   while ( param > _params[i+1] ) ++i;
1696
1697   double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
1698
1699   return getBoundaryPoints( i, r, bp1, bp2 );
1700 }
1701
1702 //================================================================================
1703 /*!
1704  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1705  *  \param [in] iMAEdge - index of a MA edge within this Branch
1706  *  \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
1707  *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1708  *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1709  *  \return bool - true if the BoundaryPoint's found
1710  */
1711 //================================================================================
1712
1713 bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t    iMAEdge,
1714                                             double         maEdgeParam,
1715                                             BoundaryPoint& bp1,
1716                                             BoundaryPoint& bp2 ) const
1717 {
1718   if ( isRemoved() )
1719     return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 );
1720
1721   if ( iMAEdge > _maEdges.size() )
1722     return false;
1723   if ( iMAEdge == _maEdges.size() )
1724     iMAEdge = _maEdges.size() - 1;
1725
1726   size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
1727   size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
1728   size_t iSeg1  = getBndSegment( _maEdges[ iMAEdge ] );
1729   size_t iSeg2  = getBndSegment( _maEdges[ iMAEdge ]->twin() );
1730
1731   return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
1732            _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
1733 }
1734
1735 //================================================================================
1736 /*!
1737  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1738  */
1739 //================================================================================
1740
1741 bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
1742                                             BoundaryPoint&     bp1,
1743                                             BoundaryPoint&     bp2 ) const
1744 {
1745   return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
1746 }
1747
1748 //================================================================================
1749 /*!
1750  * \brief Return a parameter of a BranchPoint normalized within this Branch
1751  */
1752 //================================================================================
1753
1754 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
1755 {
1756   if ( this != p._branch && p._branch )
1757     return p._branch->getParameter( p, u );
1758
1759   if ( isRemoved() )
1760     return _proxyPoint._branch->getParameter( _proxyPoint, u );
1761
1762   if ( p._iEdge > _params.size()-1 )
1763     return false;
1764   if ( p._iEdge == _params.size()-1 )
1765     return ( u = 1. );
1766
1767   u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
1768         _params[ p._iEdge+1 ] * p._edgeParam );
1769
1770   return true;
1771 }
1772
1773 //================================================================================
1774 /*!
1775  * \brief Check type of both ends
1776  */
1777 //================================================================================
1778
1779 bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
1780 {
1781   return ( _endPoint1._type == type || _endPoint2._type == type );
1782 }
1783
1784 //================================================================================
1785 /*!
1786  * \brief Returns MA points
1787  *  \param [out] points - the 2d points
1788  *  \param [in] scale - the scale that was used to scale the 2d space of MA
1789  */
1790 //================================================================================
1791
1792 void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
1793                                      const double          scale[2]) const
1794 {
1795   points.resize( _maEdges.size() + 1 );
1796
1797   points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
1798                       _maEdges[0]->vertex1()->y() / scale[1] );
1799
1800   for ( size_t i = 0; i < _maEdges.size(); ++i )
1801     points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
1802                           _maEdges[i]->vertex0()->y() / scale[1] );
1803 }
1804
1805 //================================================================================
1806 /*!
1807  * \brief Return indices of EDGEs equidistant from this branch
1808  */
1809 //================================================================================
1810
1811 void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
1812                                         std::vector< std::size_t >& edgeIDs2 ) const
1813 {
1814   edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1815   edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1816
1817   for ( size_t i = 1; i < _maEdges.size(); ++i )
1818   {
1819     size_t ie1 = getGeomEdge( _maEdges[i] );
1820     size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1821
1822     if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
1823     if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
1824   }
1825 }
1826
1827 //================================================================================
1828 /*!
1829  * \brief Looks for a BranchPoint position around a concave VERTEX
1830  */
1831 //================================================================================
1832
1833 bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&   edgeIDs1,
1834                                                    std::vector< std::size_t >&   edgeIDs2,
1835                                                    std::vector< BranchPoint >&   divPoints,
1836                                                    const vector<const TVDEdge*>& maEdges,
1837                                                    const vector<const TVDEdge*>& maEdgesTwin,
1838                                                    int &                         i) const
1839 {
1840   // if there is a concave vertex between EDGEs
1841   // then position of a dividing BranchPoint is undefined, it is somewhere
1842   // on an arc-shaped part of the Branch around the concave vertex.
1843   // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
1844   // of the arc if there is no opposite VERTEX.
1845   // All null-length segments around a VERTEX belong to one of EDGEs.
1846
1847   BranchPoint divisionPnt;
1848   divisionPnt._branch = this;
1849
1850   BranchIterator iCur( maEdges, i );
1851
1852   size_t ie1 = getGeomEdge( maEdges    [i] );
1853   size_t ie2 = getGeomEdge( maEdgesTwin[i] );
1854
1855   size_t iSeg1  = getBndSegment( iCur.edgePrev() );
1856   size_t iSeg2  = getBndSegment( iCur.edge() );
1857   bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
1858   bool isConcaNext = _boundary->isConcaveSegment( ie1,             iSeg2 );
1859   if ( !isConcaNext && !isConcaPrev )
1860     return false;
1861
1862   bool isConcaveV = false;
1863
1864   const TVDEdge* maE;
1865   BranchIterator iPrev( maEdges, i ), iNext( maEdges, i );
1866    --iPrev;
1867   if ( isConcaNext ) // all null-length segments follow
1868   {
1869     // look for a VERTEX of the opposite EDGE
1870     // iNext - next after all null-length segments
1871     while (( maE = ++iNext ))
1872     {
1873       iSeg2 = getBndSegment( maE );
1874       if ( !_boundary->isConcaveSegment( ie1, iSeg2 ))
1875         break;
1876     }
1877     bool vertexFound = false;
1878     for ( ++iCur; iCur < iNext; ++iCur )
1879     {
1880       ie2 = getGeomEdge( maEdgesTwin[ iCur.indexMod() ] );
1881       if ( ie2 != edgeIDs2.back() )
1882       {
1883         // opposite VERTEX found
1884         divisionPnt._iEdge = iCur.indexMod();
1885         divisionPnt._edgeParam = 0;
1886         divPoints.push_back( divisionPnt );
1887         edgeIDs1.push_back( ie1 );
1888         edgeIDs2.push_back( ie2 );
1889         vertexFound = true;
1890       }
1891     }
1892     if ( vertexFound )
1893     {
1894       --iNext;
1895       iPrev = iNext; // not to add a BP in the moddle
1896       i = iNext.indexMod();
1897       isConcaveV = true;
1898     }
1899   }
1900   else if ( isConcaPrev )
1901   {
1902     // all null-length segments passed, find their beginning
1903     while (( maE = iPrev.edgePrev() ))
1904     {
1905       iSeg1 = getBndSegment( maE );
1906       if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
1907         --iPrev;
1908       else
1909         break;
1910     }
1911   }
1912
1913   if ( iPrev.index() < i-1 || iNext.index() > i )
1914   {
1915     // no VERTEX on the opposite EDGE, put the Branch Point in the middle
1916     divisionPnt._iEdge = iPrev.indexMod();
1917     ++iPrev;
1918     double   par1 = _params[ iPrev.indexMod() ], par2 = _params[ iNext.indexMod() ];
1919     double midPar = 0.5 * ( par1 + par2 );
1920     for ( ; _params[ iPrev.indexMod() ] < midPar; ++iPrev )
1921       divisionPnt._iEdge = iPrev.indexMod();
1922     divisionPnt._edgeParam =
1923       ( _params[ iPrev.indexMod() ] - midPar ) /
1924       ( _params[ iPrev.indexMod() ] - _params[ divisionPnt._iEdge ] );
1925     divPoints.push_back( divisionPnt );
1926     isConcaveV = true;
1927   }
1928
1929   return isConcaveV;
1930 }
1931
1932 //================================================================================
1933 /*!
1934  * \brief Return indices of opposite parts of EDGEs equidistant from this branch
1935  *  \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
1936  *  \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
1937  *  \param [out] divPoints - BranchPoint's located between two successive unique
1938  *         pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
1939  *         of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
1940  *         than number of \a edgeIDs
1941  */
1942 //================================================================================
1943
1944 void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
1945                                                 std::vector< std::size_t >& edgeIDs2,
1946                                                 std::vector< BranchPoint >& divPoints) const
1947 {
1948   edgeIDs1.clear();
1949   edgeIDs2.clear();
1950   divPoints.clear();
1951
1952   std::vector<const TVDEdge*> twins( _maEdges.size() );
1953   for ( size_t i = 0; i < _maEdges.size(); ++i )
1954     twins[i] = _maEdges[i]->twin();
1955
1956   BranchIterator maIter ( _maEdges, 0 );
1957   BranchIterator twIter ( twins, 0 );
1958   // size_t lastConcaE1 = _boundary.nbEdges();
1959   // size_t lastConcaE2 = _boundary.nbEdges();
1960
1961   // if ( maIter._closed ) // closed branch
1962   // {
1963   //   edgeIDs1.push_back( getGeomEdge( _maEdges.back() ));
1964   //   edgeIDs2.push_back( getGeomEdge( _maEdges.back()->twin() ));
1965   // }
1966   // else
1967   {
1968     edgeIDs1.push_back( getGeomEdge( maIter.edge() ));
1969     edgeIDs2.push_back( getGeomEdge( twIter.edge() ));
1970   }
1971
1972   BranchPoint divisionPnt;
1973   divisionPnt._branch = this;
1974
1975   for ( ++maIter, ++twIter; maIter.index() < (int)_maEdges.size(); ++maIter, ++twIter )
1976   {
1977     size_t ie1 = getGeomEdge( maIter.edge() );
1978     size_t ie2 = getGeomEdge( twIter.edge() );
1979
1980     bool otherE1 = ( edgeIDs1.back() != ie1 );
1981     bool otherE2 = ( edgeIDs2.back() != ie2 );
1982
1983     if ( !otherE1 && !otherE2 && maIter._closed )
1984     {
1985       int iSegPrev1 = getBndSegment( maIter.edgePrev() );
1986       int iSegCur1  = getBndSegment( maIter.edge() );
1987       otherE1 = Abs( iSegPrev1 - iSegCur1 ) != 1;
1988       int iSegPrev2 = getBndSegment( twIter.edgePrev() );
1989       int iSegCur2  = getBndSegment( twIter.edge() );
1990       otherE2 = Abs( iSegPrev2 - iSegCur2 ) != 1;
1991     }
1992
1993     if ( otherE1 || otherE2 )
1994     {
1995       bool isConcaveV = false;
1996       if ( otherE1 && !otherE2 )
1997       {
1998         isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints,
1999                                               _maEdges, twins, maIter._i );
2000       }
2001       if ( !otherE1 && otherE2 )
2002       {
2003         isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints,
2004                                               twins, _maEdges, maIter._i );
2005       }
2006
2007       if ( isConcaveV )
2008       {
2009         ie1 = getGeomEdge( maIter.edge() );
2010         ie2 = getGeomEdge( twIter.edge() );
2011       }
2012       if ( !isConcaveV || otherE1 || otherE2 )
2013       {
2014         edgeIDs1.push_back( ie1 );
2015         edgeIDs2.push_back( ie2 );
2016       }
2017       if ( divPoints.size() < edgeIDs1.size() - 1 )
2018       {
2019         divisionPnt._iEdge = maIter.index();
2020         divisionPnt._edgeParam = 0;
2021         divPoints.push_back( divisionPnt );
2022       }
2023
2024     } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
2025   } // loop on _maEdges
2026 }
2027
2028 //================================================================================
2029 /*!
2030  * \brief Store data of boundary segments in TVDEdge
2031  */
2032 //================================================================================
2033
2034 void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
2035 {
2036   if ( maEdge ) maEdge->cell()->color( geomIndex );
2037 }
2038 std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
2039 {
2040   return maEdge ? maEdge->cell()->color() : std::string::npos;
2041 }
2042 void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
2043 {
2044   if ( maEdge ) maEdge->color( segIndex );
2045 }
2046 std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
2047 {
2048   return maEdge ? maEdge->color() : std::string::npos;
2049 }
2050
2051 //================================================================================
2052 /*!
2053  * \brief Returns a boundary point on a given EDGE
2054  *  \param [in] iEdge - index of the EDGE within MedialAxis
2055  *  \param [in] iSeg - index of a boundary segment within this Branch
2056  *  \param [in] u - [0;1] normalized param within \a iSeg-th segment
2057  *  \param [out] bp - the found BoundaryPoint
2058  *  \return bool - true if the BoundaryPoint is found
2059  */
2060 //================================================================================
2061
2062 bool SMESH_MAT2d::Boundary::getPoint( std::size_t    iEdge,
2063                                       std::size_t    iSeg,
2064                                       double         u,
2065                                       BoundaryPoint& bp ) const
2066 {
2067   if ( iEdge >= _pointsPerEdge.size() )
2068     return false;
2069   if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
2070     return false;
2071
2072   // This method is called by Branch that can have an opposite orientation,
2073   // hence u is inverted depending on orientation coded as a sign of _maEdge index
2074   bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
2075   if ( isReverse )
2076     u = 1. - u;
2077
2078   double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
2079   double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
2080
2081   bp._param = p0 * ( 1. - u ) + p1 * u;
2082   bp._edgeIndex = iEdge;
2083
2084   return true;
2085 }
2086