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