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