src/SMESHUtils/SMESH_MAT2d.cxx

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