src/SMESHUtils/SMESH_MAT2d.cxx

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