src/SMESHUtils/SMESH_MAT2d.cxx

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