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