src/SMESHUtils/SMESH_MAT2d.cxx

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