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;
  70     double _param;
  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
 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       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         vector< double > len;
 802         len.push_back(0);
 803         for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e )
 804           len.push_back( len.back() + length( *e ));
 805
 806         e = inSeg._edges.rbegin();
 807         for ( size_t l = 1; l < len.size(); ++e, ++l )
 808         {
 809           double dl = len[l] / len.back();
 810           double u  = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param;
 811           bndSegs.push_back( BndSeg( &inSeg, *e, u ));
 812         }
 813       }
 814       // segments around 2nd concave point
 815       size_t ip1 = inSeg._p1->index( inPoints );
 816       if ( inPntChecked[ ip1 ] )
 817         for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
 818           bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
 819       inPntChecked[ ip1 ] = false;
 820     }
 821
 822     // make TVDEdge's know it's BndSeg to enable passing branchID to
 823     // an opposite BndSeg in BndSeg::setBranch()
 824     for ( size_t i = 0; i < bndSegs.size(); ++i )
 825       bndSegs[i].setIndexToEdge( i );
 826
 827
 828     // Find TVDEdge's of Branches and associate them with bndSegs
 829
 830     vector< vector<const TVDEdge*> > branchEdges;
 831     branchEdges.reserve( boundary.nbEdges() * 4 );
 832
 833     map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
 834
 835     int branchID = 1; // we code orientation as branchID sign
 836     branchEdges.resize( branchID + 1 );
 837
 838     size_t i1st = 0;
 839     while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
 840       ++i1st;
 841     bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and the opposite bndSeg
 842     branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
 843
 844     for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
 845     {
 846       if ( bndSegs[i].branchID() )
 847       {
 848         branchID = bndSegs[i]._branchID; // with sign
 849
 850         if ( bndSegs[i]._branchID == -bndSegs[i-1]._branchID &&
 851              bndSegs[i]._edge )
 852         {
 853           SMESH_MAT2d::BranchEndType type =
 854             ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
 855               SMESH_MAT2d::BE_ON_VERTEX :
 856               SMESH_MAT2d::BE_END );
 857           endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
 858         }
 859         continue;
 860       }
 861       if ( !bndSegs[i-1].isSameBranch( bndSegs[i] ))
 862       {
 863         branchEdges.resize(( branchID = branchEdges.size()) + 1 );
 864         if ( bndSegs[i]._edge )
 865           endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
 866                                      SMESH_MAT2d::BE_BRANCH_POINT ));
 867       }
 868       bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
 869       if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
 870         branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
 871     }
 872     // define BranchEndType of the first TVDVertex
 873     if ( bndSegs.front()._branchID == -bndSegs.back()._branchID )
 874     {
 875       if ( bndSegs[0]._edge )
 876       {
 877         SMESH_MAT2d::BranchEndType type =
 878           ( bndSegs[0]._inSeg->isConnected( bndSegs[0]._edge ) ?
 879             SMESH_MAT2d::BE_ON_VERTEX :
 880             SMESH_MAT2d::BE_END );
 881         endType.insert( make_pair( bndSegs[0]._edge->vertex1(), type ));
 882       }
 883       else if ( bndSegs.back()._edge )
 884       {
 885         SMESH_MAT2d::BranchEndType type =
 886           ( bndSegs.back()._inSeg->isConnected( bndSegs.back()._edge ) ?
 887             SMESH_MAT2d::BE_ON_VERTEX :
 888             SMESH_MAT2d::BE_END );
 889         endType.insert( make_pair( bndSegs.back()._edge->vertex0(), type ));
 890       }
 891     }
 892     // join the 1st and the last branch edges if it is the same branch
 893     if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
 894          bndSegs.back().isSameBranch( bndSegs.front() ))
 895     {
 896       vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
 897       vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID()  ];
 898       br1.insert( br1.begin(), br2.begin(), br2.end() );
 899       br2.clear();
 900     }
 901
 902     // remove branches ending at BE_ON_VERTEX
 903
 904     vector<bool> isBranchRemoved( branchEdges.size(), false );
 905
 906     if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
 907     {
 908       // find branches to remove
 909       map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et;
 910       for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
 911       {
 912         if ( branchEdges[iB].empty() )
 913           continue;
 914         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
 915         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
 916         v2et = endType.find( v0 );
 917         if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
 918           isBranchRemoved[ iB ] = true;
 919         v2et = endType.find( v1 );
 920         if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
 921           isBranchRemoved[ iB ] = true;
 922       }
 923       // try to join not removed branches into one
 924       for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
 925       {
 926         if ( branchEdges[iB].empty() || isBranchRemoved[iB] )
 927           continue;
 928         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
 929         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
 930         v2et = endType.find( v0 );
 931         if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
 932           v0 = 0;
 933         v2et = endType.find( v1 );
 934         if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
 935           v1 = 0;
 936         if ( !v0 && !v1 )
 937           continue;
 938
 939         size_t iBrToJoin = 0;
 940         for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
 941         {
 942           if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
 943             continue;
 944           const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
 945           const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
 946           if ( v0 == v02 || v0 == v12 || v1 == v02 || v1 == v12 )
 947           {
 948             if ( iBrToJoin > 0 )
 949             {
 950               iBrToJoin = 0;
 951               break; // more than 2 not removed branches meat at a TVDVertex
 952             }
 953             iBrToJoin = iB2;
 954           }
 955         }
 956         if ( iBrToJoin > 0 )
 957         {
 958           vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
 959           const TVDVertex* v02 = branch[0]->vertex1();
 960           const TVDVertex* v12 = branch.back()->vertex0();
 961           updateJoinedBranch( branch, iB, bndSegs, /*reverse=*/(v0 == v02 || v1 == v12 ));
 962           if ( v0 == v02 || v0 == v12 )
 963             branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
 964           else
 965             branchEdges[iB].insert( branchEdges[iB].end(),   branch.begin(), branch.end() );
 966           branch.clear();
 967         }
 968       } // loop on branchEdges
 969     } // if ( ignoreCorners )
 970
 971     // associate branchIDs and the input branch vector (arg)
 972     vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 );
 973     int nbBranches = 0;
 974     for ( size_t i = 0; i < branchEdges.size(); ++i )
 975     {
 976       nbBranches += ( !branchEdges[i].empty() );
 977     }
 978     branch.resize( nbBranches );
 979     size_t iBr = 0;
 980     for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed
 981     {
 982       if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ])
 983         branchByID[ brID ] = & branch[ iBr++ ];
 984     }
 985     for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed
 986     {
 987       if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ])
 988         branchByID[ brID ] = & branch[ iBr++ ];
 989     }
 990
 991     // Fill in BndPoints of each EDGE of the boundary
 992
 993     size_t iSeg = 0;
 994     int edgeInd = -1, dInd = 0;
 995     while ( iSeg < bndSegs.size() )
 996     {
 997       const size_t                geomID = bndSegs[ iSeg ].geomEdge();
 998       SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( geomID );
 999
1000       size_t nbSegs = 0;
1001       for ( size_t i = iSeg; i < bndSegs.size() && geomID == bndSegs[ i ].geomEdge(); ++i )
1002         ++nbSegs;
1003       size_t iSegEnd = iSeg + nbSegs;
1004
1005       // make TVDEdge know an index of bndSegs within BndPoints
1006       for ( size_t i = iSeg; i < iSegEnd; ++i )
1007         if ( bndSegs[i]._edge )
1008           SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge );
1009
1010       // parameters on EDGE
1011
1012       bndPoints._params.reserve( nbSegs + 1 );
1013       bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param );
1014
1015       for ( size_t i = iSeg; i < iSegEnd; ++i )
1016         bndPoints._params.push_back( bndSegs[ i ]._uLast );
1017
1018       // MA edges
1019
1020       bndPoints._maEdges.reserve( nbSegs );
1021
1022       for ( size_t i = iSeg; i < iSegEnd; ++i )
1023       {
1024         const size_t              brID = bndSegs[ i ].branchID();
1025         const SMESH_MAT2d::Branch*  br = branchByID[ brID ];
1026
1027         if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
1028         {
1029           edgeInd += dInd;
1030
1031           if (( edgeInd < 0 ||
1032                 edgeInd >= (int) branchEdges[ brID ].size() ) ||
1033               ( branchEdges[ brID ][ edgeInd ]         != bndSegs[ i ]._edge &&
1034                 branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge ))
1035           {
1036             if ( bndSegs[ i ]._branchID < 0 )
1037             {
1038               dInd = -1;
1039               for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd )
1040                 if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge )
1041                   break;
1042             }
1043             else // bndSegs[ i ]._branchID > 0
1044             {
1045               dInd = +1;
1046               for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
1047                 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
1048                   break;
1049             }
1050           }
1051         }
1052         else
1053         {
1054           // no MA edge, bndSeg corresponds to an end point of a branch
1055           if ( bndPoints._maEdges.empty() )
1056           {
1057             // should not get here according to algo design???
1058             edgeInd = 0;
1059           }
1060           else
1061           {
1062             edgeInd = branchEdges[ brID ].size();
1063             dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
1064           }
1065         }
1066
1067         bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
1068
1069       } // loop on bndSegs of an EDGE
1070
1071       iSeg = iSegEnd;
1072
1073     } // loop on all bndSegs
1074
1075     // Initialize branches
1076
1077     // find a not removed branch
1078     size_t iBrNorRemoved = 0;
1079     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1080       if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] )
1081       {
1082         iBrNorRemoved = brID;
1083         break;
1084       }
1085     // fill the branches with MA edges
1086     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1087       if ( !branchEdges[brID].empty() )
1088       {
1089         branchByID[ brID ]->init( branchEdges[brID], & boundary, endType );
1090       }
1091     // mark removed branches
1092     for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1093       if ( isBranchRemoved[brID] && iBrNorRemoved > 0 )
1094       {
1095         SMESH_MAT2d::Branch* branch     = branchByID[ brID ];
1096         SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ];
1097         bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT );
1098         const TVDVertex* branchVextex =
1099           is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex;
1100         SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex );
1101         branch->setRemoved( bp );
1102       }
1103     // set branches to branch ends
1104     for ( size_t i = 0; i < branch.size(); ++i )
1105       if ( !branch[i].isRemoved() )
1106         branch[i].setBranchesToEnds( branch );
1107
1108     // fill branchPnt arg
1109     map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
1110     for ( size_t i = 0; i < branch.size(); ++i )
1111     {
1112       if ( branch[i].getEnd(0)->_branches.size() > 2 )
1113         v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
1114       if ( branch[i].getEnd(1)->_branches.size() > 2 )
1115         v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
1116     }
1117     branchPnt.resize( v2end.size() );
1118     map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
1119     for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
1120       branchPnt[ i ] = v2e->second;
1121
1122   } // makeMA()
1123
1124 } // namespace
1125
1126 //================================================================================
1127 /*!
1128  * \brief MedialAxis constructor
1129  *  \param [in] face - a face to create MA for
1130  *  \param [in] edges - edges of the face (possibly not all) on the order they
1131  *              encounter in the face boundary.
1132  *  \param [in] minSegLen - minimal length of a mesh segment used to discretize
1133  *              the edges. It is used to define precision of MA approximation
1134  */
1135 //================================================================================
1136
1137 SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face&                face,
1138                                     const std::vector< TopoDS_Edge >& edges,
1139                                     const double                      minSegLen,
1140                                     const bool                        ignoreCorners):
1141   _face( face ), _boundary( edges.size() )
1142 {
1143   // input to construct_voronoi()
1144   vector< InPoint >  inPoints;
1145   vector< InSegment> inSegments;
1146   if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
1147     return;
1148
1149   inSegmentsToFile( inSegments );
1150
1151   // build voronoi diagram
1152   construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
1153
1154   // make MA data
1155   makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary );
1156
1157   // count valid branches
1158   _nbBranches = _branch.size();
1159   for ( size_t i = 0; i < _branch.size(); ++i )
1160     if ( _branch[i].isRemoved() )
1161       --_nbBranches;
1162 }
1163
1164 //================================================================================
1165 /*!
1166  * \brief Returns the i-th branch
1167  */
1168 //================================================================================
1169
1170 const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const
1171 {
1172   return i < _nbBranches ? &_branch[i] : 0;
1173 }
1174
1175 //================================================================================
1176 /*!
1177  * \brief Return UVs of ends of MA edges of a branch
1178  */
1179 //================================================================================
1180
1181 void SMESH_MAT2d::MedialAxis::getPoints( const Branch*         branch,
1182                                          std::vector< gp_XY >& points) const
1183 {
1184   branch->getPoints( points, _scale );
1185 }
1186
1187 //================================================================================
1188 /*!
1189  * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
1190  *  \param [in] iGeomEdge - index of geom EDGE within a vector passed at MA construction
1191  *  \param [in] u - parameter of the point on EDGE curve
1192  *  \param [out] p - the found BranchPoint
1193  *  \return bool - is OK
1194  */
1195 //================================================================================
1196
1197 bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
1198                                             double            u,
1199                                             BranchPoint&      p ) const
1200 {
1201   if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1202     return false;
1203
1204   const BndPoints& points = _pointsPerEdge[ iEdge ];
1205   const bool  edgeReverse = ( points._params[0] > points._params.back() );
1206
1207   if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
1208     u = edgeReverse ? points._params.back() : points._params[0];
1209   else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
1210     u = edgeReverse ? points._params[0] : points._params.back();
1211
1212   double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
1213   int    i = int( r * double( points._maEdges.size()-1 ));
1214   if ( edgeReverse )
1215   {
1216     while ( points._params[i  ] < u ) --i;
1217     while ( points._params[i+1] > u ) ++i;
1218   }
1219   else
1220   {
1221     while ( points._params[i  ] > u ) --i;
1222     while ( points._params[i+1] < u ) ++i;
1223   }
1224
1225   double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
1226
1227   if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
1228   {
1229     if ( i < points._maEdges.size() / 2 ) // near 1st point
1230     {
1231       while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second )
1232         ++i;
1233       edgeParam = edgeReverse;
1234     }
1235     else // near last point
1236     {
1237       while ( i > 0 && !points._maEdges[ i ].second )
1238         --i;
1239       edgeParam = !edgeReverse;
1240     }
1241   }
1242   const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
1243   bool maReverse = ( maE.second < 0 );
1244
1245   p._branch = maE.first;
1246   p._iEdge  = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
1247   p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam;
1248
1249   return true;
1250 }
1251
1252 //================================================================================
1253 /*!
1254  * \brief Check if a given boundary segment is a null-length segment on a concave
1255  *        boundary corner.
1256  *  \param [in] iEdge - index of a geom EDGE
1257  *  \param [in] iSeg - index of a boundary segment
1258  *  \return bool - true if the segment is on concave corner
1259  */
1260 //================================================================================
1261
1262 bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
1263 {
1264   if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1265     return false;
1266
1267   const BndPoints& points = _pointsPerEdge[ iEdge ];
1268   if ( points._params.size() <= iSeg+1 )
1269     return false;
1270
1271   return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20;
1272 }
1273
1274 //================================================================================
1275 /*!
1276  * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end
1277  */
1278 //================================================================================
1279
1280 bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
1281 {
1282   if ( bp._edgeIndex >= _pointsPerEdge.size() )
1283     return false;
1284
1285   const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
1286   if ( bp._param - points._params[0] < points._params.back() - bp._param )
1287     bp._param = points._params[0];
1288   else
1289     bp._param = points._params.back();
1290
1291   return true;
1292 }
1293
1294 //================================================================================
1295 /*!
1296  * \brief Creates a 3d curve corresponding to a Branch
1297  *  \param [in] branch - the Branch
1298  *  \return Adaptor3d_Curve* - the new curve the caller is to delete
1299  */
1300 //================================================================================
1301
1302 Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
1303 {
1304   Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
1305   if ( surface.IsNull() )
1306     return 0;
1307
1308   vector< gp_XY > uv;
1309   branch.getPoints( uv, _scale );
1310   if ( uv.size() < 2 )
1311     return 0;
1312
1313   vector< TopoDS_Vertex > vertex( uv.size() );
1314   for ( size_t i = 0; i < uv.size(); ++i )
1315     vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
1316
1317   TopoDS_Wire aWire;
1318   BRep_Builder aBuilder;
1319   aBuilder.MakeWire(aWire);
1320   for ( size_t i = 1; i < vertex.size(); ++i )
1321   {
1322     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
1323     aBuilder.Add( aWire, edge );
1324   }
1325
1326   // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
1327   //   aWire.Closed(true); // issue 0021141
1328
1329   return new BRepAdaptor_CompCurve( aWire );
1330 }
1331
1332 //================================================================================
1333 /*!
1334  * \brief Copy points of an EDGE
1335  */
1336 //================================================================================
1337
1338 void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                maEdges,
1339                                 const Boundary*                        boundary,
1340                                 map< const TVDVertex*, BranchEndType > endType )
1341 {
1342   if ( maEdges.empty() ) return;
1343
1344   _boundary = boundary;
1345   _maEdges.swap( maEdges );
1346
1347
1348   _params.reserve( _maEdges.size() + 1 );
1349   _params.push_back( 0. );
1350   for ( size_t i = 0; i < _maEdges.size(); ++i )
1351     _params.push_back( _params.back() + length( _maEdges[i] ));
1352
1353   for ( size_t i = 1; i < _params.size(); ++i )
1354     _params[i] /= _params.back();
1355
1356
1357   _endPoint1._vertex = _maEdges.front()->vertex1();
1358   _endPoint2._vertex = _maEdges.back ()->vertex0();
1359
1360   if ( endType.count( _endPoint1._vertex ))
1361     _endPoint1._type = endType[ _endPoint1._vertex ];
1362   if ( endType.count( _endPoint2._vertex ))
1363     _endPoint2._type = endType[ _endPoint2._vertex ];
1364 }
1365
1366 //================================================================================
1367 /*!
1368  * \brief fills BranchEnd::_branches of its ends
1369  */
1370 //================================================================================
1371
1372 void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
1373 {
1374   for ( size_t i = 0; i < branches.size(); ++i )
1375   {
1376     if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
1377          this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
1378       this->_endPoint1._branches.push_back( &branches[i] );
1379
1380     if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
1381          this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
1382       this->_endPoint2._branches.push_back( &branches[i] );
1383   }
1384 }
1385
1386 //================================================================================
1387 /*!
1388  * \brief returns a BranchPoint corresponding to a TVDVertex
1389  */
1390 //================================================================================
1391
1392 SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const
1393 {
1394   BranchPoint p;
1395   p._branch = this;
1396   p._iEdge  = 0;
1397
1398   if ( vertex == _maEdges[0]->vertex1() )
1399   {
1400     p._edgeParam = 0;
1401   }
1402   else
1403   {
1404     for ( ; p._iEdge < _maEdges.size(); ++p._iEdge )
1405       if ( vertex == _maEdges[ p._iEdge ]->vertex0() )
1406       {
1407         p._edgeParam = _params[ p._iEdge ];
1408         break;
1409       }
1410   }
1411   return p;
1412 }
1413
1414 //================================================================================
1415 /*!
1416  * \brief Sets a proxy point for a removed branch
1417  *  \param [in] proxyPoint - a point of another branch to which all points of this
1418  *         branch are mapped
1419  */
1420 //================================================================================
1421
1422 void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint )
1423 {
1424   _proxyPoint = proxyPoint;
1425 }
1426
1427 //================================================================================
1428 /*!
1429  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1430  *  \param [in] param - [0;1] normalized param on the Branch
1431  *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1432  *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1433  *  \return bool - true if the BoundaryPoint's found
1434  */
1435 //================================================================================
1436
1437 bool SMESH_MAT2d::Branch::getBoundaryPoints(double         param,
1438                                             BoundaryPoint& bp1,
1439                                             BoundaryPoint& bp2 ) const
1440 {
1441   if ( param < _params[0] || param > _params.back() )
1442     return false;
1443
1444   // look for an index of a MA edge by param
1445   double ip = param * _params.size();
1446   size_t  i = size_t( Min( int( _maEdges.size()-1), int( ip )));
1447
1448   while ( param < _params[i  ] ) --i;
1449   while ( param > _params[i+1] ) ++i;
1450
1451   double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
1452
1453   return getBoundaryPoints( i, r, bp1, bp2 );
1454 }
1455
1456 //================================================================================
1457 /*!
1458  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1459  *  \param [in] iMAEdge - index of a MA edge within this Branch
1460  *  \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
1461  *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1462  *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1463  *  \return bool - true if the BoundaryPoint's found
1464  */
1465 //================================================================================
1466
1467 bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t    iMAEdge,
1468                                             double         maEdgeParam,
1469                                             BoundaryPoint& bp1,
1470                                             BoundaryPoint& bp2 ) const
1471 {
1472   if ( isRemoved() )
1473     return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 );
1474
1475   if ( iMAEdge > _maEdges.size() )
1476     return false;
1477   if ( iMAEdge == _maEdges.size() )
1478     iMAEdge = _maEdges.size() - 1;
1479
1480   size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
1481   size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
1482   size_t iSeg1  = getBndSegment( _maEdges[ iMAEdge ] );
1483   size_t iSeg2  = getBndSegment( _maEdges[ iMAEdge ]->twin() );
1484
1485   return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
1486            _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
1487 }
1488
1489 //================================================================================
1490 /*!
1491  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1492  */
1493 //================================================================================
1494
1495 bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
1496                                             BoundaryPoint&     bp1,
1497                                             BoundaryPoint&     bp2 ) const
1498 {
1499   return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
1500 }
1501
1502 //================================================================================
1503 /*!
1504  * \brief Return a parameter of a BranchPoint normalized within this Branch
1505  */
1506 //================================================================================
1507
1508 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
1509 {
1510   if ( p._iEdge > _params.size()-1 )
1511     return false;
1512   if ( p._iEdge == _params.size()-1 )
1513     return u = 1.;
1514
1515   u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
1516         _params[ p._iEdge+1 ] * p._edgeParam );
1517
1518   return true;
1519 }
1520
1521 //================================================================================
1522 /*!
1523  * \brief Check type of both ends
1524  */
1525 //================================================================================
1526
1527 bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
1528 {
1529   return ( _endPoint1._type == type || _endPoint2._type == type );
1530 }
1531
1532 //================================================================================
1533 /*!
1534  * \brief Returns MA points
1535  *  \param [out] points - the 2d points
1536  *  \param [in] scale - the scale that was used to scale the 2d space of MA
1537  */
1538 //================================================================================
1539
1540 void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
1541                                      const double          scale[2]) const
1542 {
1543   points.resize( _maEdges.size() + 1 );
1544
1545   points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
1546                       _maEdges[0]->vertex1()->y() / scale[1] );
1547
1548   for ( size_t i = 0; i < _maEdges.size(); ++i )
1549     points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
1550                           _maEdges[i]->vertex0()->y() / scale[1] );
1551 }
1552
1553 //================================================================================
1554 /*!
1555  * \brief Return indices of EDGEs equidistant from this branch
1556  */
1557 //================================================================================
1558
1559 void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
1560                                         std::vector< std::size_t >& edgeIDs2 ) const
1561 {
1562   edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1563   edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1564
1565   for ( size_t i = 1; i < _maEdges.size(); ++i )
1566   {
1567     size_t ie1 = getGeomEdge( _maEdges[i] );
1568     size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1569
1570     if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
1571     if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
1572   }
1573 }
1574
1575 //================================================================================
1576 /*!
1577  * \brief Looks for a BranchPoint position around a concave VERTEX
1578  */
1579 //================================================================================
1580
1581 bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&   edgeIDs1,
1582                                                    std::vector< std::size_t >&   edgeIDs2,
1583                                                    std::vector< BranchPoint >&   divPoints,
1584                                                    const vector<const TVDEdge*>& maEdges,
1585                                                    const vector<const TVDEdge*>& maEdgesTwin,
1586                                                    size_t &                    i) const
1587 {
1588   // if there is a concave vertex between EDGEs
1589   // then position of a dividing BranchPoint is undefined, it is somewhere
1590   // on an arc-shaped part of the Branch around the concave vertex.
1591   // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
1592   // of the arc if there is no opposite VERTEX.
1593   // All null-length segments around a VERTEX belong to one of EDGEs.
1594
1595   BranchPoint divisionPnt;
1596   divisionPnt._branch = this;
1597
1598   size_t ie1 = getGeomEdge( maEdges    [i] );
1599   size_t ie2 = getGeomEdge( maEdgesTwin[i] );
1600
1601   size_t iSeg1  = getBndSegment( maEdges[ i-1 ] );
1602   size_t iSeg2  = getBndSegment( maEdges[ i ] );
1603   bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
1604   bool isConcaNext = _boundary->isConcaveSegment( ie1,             iSeg2 );
1605   if ( !isConcaNext && !isConcaPrev )
1606     return false;
1607
1608   bool isConcaveV = false;
1609
1610   int iPrev = i-1, iNext = i;
1611   if ( isConcaNext ) // all null-length segments follow
1612   {
1613     // look for a VERTEX of the opposite EDGE
1614     ++iNext; // end of null-length segments
1615     while ( iNext < maEdges.size() )
1616     {
1617       iSeg2 = getBndSegment( maEdges[ iNext ] );
1618       if ( _boundary->isConcaveSegment( ie1, iSeg2 ))
1619         ++iNext;
1620       else
1621         break;
1622     }
1623     bool vertexFound = false;
1624     for ( size_t iE = i+1; iE < iNext; ++iE )
1625     {
1626       ie2 = getGeomEdge( maEdgesTwin[iE] );
1627       if ( ie2 != edgeIDs2.back() )
1628       {
1629         // opposite VERTEX found
1630         divisionPnt._iEdge = iE;
1631         divisionPnt._edgeParam = 0;
1632         divPoints.push_back( divisionPnt );
1633         edgeIDs1.push_back( ie1 );
1634         edgeIDs2.push_back( ie2 );
1635         vertexFound = true;
1636         }
1637     }
1638     if ( vertexFound )
1639     {
1640       iPrev = i = --iNext; // not to add a BP in the moddle
1641       isConcaveV = true;
1642     }
1643   }
1644   else if ( isConcaPrev )
1645   {
1646     // all null-length segments passed, find their beginning
1647     while ( iPrev-1 >= 0 )
1648     {
1649       iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
1650       if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
1651         --iPrev;
1652       else
1653         break;
1654     }
1655   }
1656
1657   if ( iPrev < i-1 || iNext > i )
1658   {
1659     // no VERTEX on the opposite EDGE, put the Branch Point in the middle
1660     double par1 = _params[ iPrev+1 ], par2 = _params[ iNext ];
1661     double midPar = 0.5 * ( par1 + par2 );
1662     divisionPnt._iEdge = iPrev;
1663     while ( _params[ divisionPnt._iEdge + 1 ] < midPar )
1664       ++divisionPnt._iEdge;
1665     divisionPnt._edgeParam =
1666       ( _params[ divisionPnt._iEdge + 1 ] - midPar ) /
1667       ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] );
1668     divPoints.push_back( divisionPnt );
1669     isConcaveV = true;
1670   }
1671
1672   return isConcaveV;
1673 }
1674
1675 //================================================================================
1676 /*!
1677  * \brief Return indices of opposite parts of EDGEs equidistant from this branch
1678  *  \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
1679  *  \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
1680  *  \param [out] divPoints - BranchPoint's located between two successive unique
1681  *         pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
1682  *         of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
1683  *         than number of \a edgeIDs
1684  */
1685 //================================================================================
1686
1687 void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
1688                                                 std::vector< std::size_t >& edgeIDs2,
1689                                                 std::vector< BranchPoint >& divPoints) const
1690 {
1691   edgeIDs1.clear();
1692   edgeIDs2.clear();
1693   divPoints.clear();
1694
1695   edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1696   edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1697
1698   std::vector<const TVDEdge*> twins( _maEdges.size() );
1699   for ( size_t i = 0; i < _maEdges.size(); ++i )
1700     twins[i] = _maEdges[i]->twin();
1701
1702   // size_t lastConcaE1 = _boundary.nbEdges();
1703   // size_t lastConcaE2 = _boundary.nbEdges();
1704
1705   BranchPoint divisionPnt;
1706   divisionPnt._branch = this;
1707
1708   for ( size_t i = 0; i < _maEdges.size(); ++i )
1709   {
1710     size_t ie1 = getGeomEdge( _maEdges[i] );
1711     size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1712
1713     if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
1714     {
1715       bool isConcaveV = false;
1716       if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 )
1717       {
1718         isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i );
1719       }
1720       if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 )
1721       {
1722         isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i );
1723       }
1724
1725       if ( isConcaveV )
1726       {
1727         ie1 = getGeomEdge( _maEdges[i] );
1728         ie2 = getGeomEdge( _maEdges[i]->twin() );
1729       }
1730       if (( !isConcaveV ) ||
1731           ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ))
1732       {
1733         edgeIDs1.push_back( ie1 );
1734         edgeIDs2.push_back( ie2 );
1735       }
1736       if ( divPoints.size() < edgeIDs1.size() - 1 )
1737       {
1738         divisionPnt._iEdge = i;
1739         divisionPnt._edgeParam = 0;
1740         divPoints.push_back( divisionPnt );
1741       }
1742
1743     } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
1744   } // loop on _maEdges
1745 }
1746
1747 //================================================================================
1748 /*!
1749  * \brief Store data of boundary segments in TVDEdge
1750  */
1751 //================================================================================
1752
1753 void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
1754 {
1755   if ( maEdge ) maEdge->cell()->color( geomIndex );
1756 }
1757 std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
1758 {
1759   return maEdge ? maEdge->cell()->color() : std::string::npos;
1760 }
1761 void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
1762 {
1763   if ( maEdge ) maEdge->color( segIndex );
1764 }
1765 std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
1766 {
1767   return maEdge ? maEdge->color() : std::string::npos;
1768 }
1769
1770 //================================================================================
1771 /*!
1772  * \brief Returns a boundary point on a given EDGE
1773  *  \param [in] iEdge - index of the EDGE within MedialAxis
1774  *  \param [in] iSeg - index of a boundary segment within this Branch
1775  *  \param [in] u - [0;1] normalized param within \a iSeg-th segment
1776  *  \param [out] bp - the found BoundaryPoint
1777  *  \return bool - true if the BoundaryPoint is found
1778  */
1779 //================================================================================
1780
1781 bool SMESH_MAT2d::Boundary::getPoint( std::size_t    iEdge,
1782                                       std::size_t    iSeg,
1783                                       double         u,
1784                                       BoundaryPoint& bp ) const
1785 {
1786   if ( iEdge >= _pointsPerEdge.size() )
1787     return false;
1788   if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
1789     return false;
1790
1791   // This method is called by Branch that can have an opposite orientation,
1792   // hence u is inverted depending on orientation coded as a sign of _maEdge index
1793   bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
1794   if ( isReverse )
1795     u = 1. - u;
1796
1797   double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
1798   double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
1799
1800   bp._param = p0 * ( 1. - u ) + p1 * u;
1801   bp._edgeIndex = iEdge;
1802
1803   return true;
1804 }
1805