src/SMESHUtils/SMESH_MAT2d.cxx

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