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