1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : SMESH_MAT2d.cxx
23 // Created : Thu May 28 17:49:53 2015
24 // Author : Edward AGAPOV (eap)
26 #include "SMESH_MAT2d.hxx"
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>
46 #include <TopoDS_Vertex.hxx>
47 #include <TopoDS_Wire.hxx>
51 #include "SMESH_File.hxx"
52 #include "SMESH_Comment.hxx"
56 using boost::polygon::x;
57 using boost::polygon::y;
58 using SMESH_MAT2d::TVD;
59 using SMESH_MAT2d::TVDEdge;
60 using SMESH_MAT2d::TVDCell;
61 using SMESH_MAT2d::TVDVertex;
65 // Input data for construct_voronoi()
66 // -------------------------------------------------------------------------------------
70 int _a, _b; // coordinates
71 double _param; // param on EDGE
72 InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
73 InPoint() : _a(0), _b(0), _param(0) {}
76 list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order
78 size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
79 bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
80 bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
81 Abs( _b - v->y() ) < 1. ); }
83 // -------------------------------------------------------------------------------------
91 size_t _geomEdgeInd; // EDGE index within the FACE
93 list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
95 InSegment( InPoint * p0, InPoint * p1, size_t iE)
96 : _p0(p0), _p1(p1), _geomEdgeInd(iE), _cell(0) {}
97 InSegment() : _p0(0), _p1(0), _geomEdgeInd(0), _cell(0) {}
99 const InPoint& point0() const { return *_p0; }
100 const InPoint& point1() const { return *_p1; }
102 inline bool isConnected( const TVDEdge* edge );
104 inline bool isExternal( const TVDEdge* edge );
106 static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); }
108 static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); }
111 // check if a TVDEdge begins at my end or ends at my start
112 inline bool InSegment::isConnected( const TVDEdge* edge )
114 return (( edge->vertex0() && edge->vertex1() )
116 ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
117 Abs( edge->vertex0()->y() - _p1->_b ) < 1. ) ||
118 (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
119 Abs( edge->vertex1()->y() - _p0->_b ) < 1. )));
122 // check if a MA TVDEdge is outside of a domain
123 inline bool InSegment::isExternal( const TVDEdge* edge )
125 double dot = // x1*x2 + y1*y2; (x1,y1) - internal normal of InSegment
126 ( _p0->_b - _p1->_b ) * ( 0.5 * ( edge->vertex0()->x() + edge->vertex1()->x() ) - _p0->_a ) +
127 ( _p1->_a - _p0->_a ) * ( 0.5 * ( edge->vertex0()->y() + edge->vertex1()->y() ) - _p0->_b );
131 // // -------------------------------------------------------------------------------------
132 // const size_t theExternMA = 111; // to mark external MA edges
134 // bool isExternal( const TVDEdge* edge )
136 // return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA );
139 // // mark external MA edges
140 // void markExternalEdges( const TVDEdge* edge )
142 // if ( isExternal( edge ))
144 // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge );
145 // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() );
146 // if ( edge->is_primary() && edge->vertex1() )
148 // const TVDVertex * v = edge->vertex1();
149 // edge = v->incident_edge();
151 // markExternalEdges( edge );
152 // edge = edge->rot_next();
153 // } while ( edge != v->incident_edge() );
157 // -------------------------------------------------------------------------------------
159 // writes segments into a txt file readable by voronoi_visualizer
160 void inSegmentsToFile( vector< InSegment>& inSegments)
162 if ( inSegments.size() > 1000 )
164 const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
165 const char* user = getenv("USER");
166 if ( !user || strcmp( user, "eap" )) return;
167 SMESH_File file(fileName, false );
169 file.openForWriting();
171 text << "0\n"; // nb points
172 text << inSegments.size() << "\n"; // nb segments
173 for ( size_t i = 0; i < inSegments.size(); ++i )
175 text << inSegments[i]._p0->_a << " "
176 << inSegments[i]._p0->_b << " "
177 << inSegments[i]._p1->_a << " "
178 << inSegments[i]._p1->_b << "\n";
181 file.write( text.c_str(), text.size() );
182 cout << "Write " << fileName << endl;
184 void dumpEdge( const TVDEdge* edge )
186 cout << "*Edge_" << edge;
187 if ( !edge->vertex0() )
188 cout << " ( INF, INF";
190 cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y();
191 if ( !edge->vertex1() )
192 cout << ") -> ( INF, INF";
194 cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
195 cout << ")\t cell=" << edge->cell()
196 << " iBnd=" << edge->color()
197 << " twin=" << edge->twin()
198 << " twin_cell=" << edge->twin()->cell()
199 << " prev=" << edge->prev() << " next=" << edge->next()
200 << ( edge->is_primary() ? " MA " : " SCND" )
201 << ( edge->is_linear() ? " LIN " : " CURV" )
204 void dumpCell( const TVDCell* cell )
206 cout << "**Cell_" << cell << " GEOM=" << cell->color() << " ";
207 cout << ( cell->contains_segment() ? " SEG " : " PNT " );
208 if ( cell-> is_degenerate() )
213 const TVDEdge* edge = cell->incident_edge();
217 cout << " - " << ++i << " ";
219 } while (edge != cell->incident_edge());
223 #define inSegmentsToFile(arg) {}
224 //void dumpEdge( const TVDEdge* edge ) {}
225 //void dumpCell( const TVDCell* cell ) {}
228 // -------------------------------------------------------------------------------------
234 struct geometry_concept<InPoint> {
235 typedef point_concept type;
238 struct point_traits<InPoint> {
239 typedef int coordinate_type;
241 static inline coordinate_type get(const InPoint& point, orientation_2d orient) {
242 return (orient == HORIZONTAL) ? point._a : point._b;
247 struct geometry_concept<InSegment> {
248 typedef segment_concept type;
252 struct segment_traits<InSegment> {
253 typedef int coordinate_type;
254 typedef InPoint point_type;
256 static inline point_type get(const InSegment& segment, direction_1d dir) {
257 return *(dir.to_int() ? segment._p1 : segment._p0);
260 } // namespace polygon
262 // -------------------------------------------------------------------------------------
266 const int theNoBrachID = 0;
267 double theScale[2]; // scale used in bndSegsToMesh()
268 const size_t theNoEdgeID = std::numeric_limits<size_t>::max() / 1000;
270 // -------------------------------------------------------------------------------------
272 * \brief Intermediate DS to create InPoint's
278 UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {}
279 InPoint getInPoint( double scale[2] )
281 return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u );
284 // -------------------------------------------------------------------------------------
286 * \brief Segment of EDGE, used to create BndPoints
291 const TVDEdge* _edge;
293 BndSeg* _prev; // previous BndSeg in FACE boundary
294 int _branchID; // negative ID means reverse direction
296 BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
297 _inSeg(seg), _edge(edge), _uLast(u), _prev(0), _branchID( theNoBrachID ) {}
299 void setIndexToEdge( size_t id )
301 SMESH_MAT2d::Branch::setBndSegment( id, _edge );
304 int branchID() const { return Abs( _branchID ); }
306 size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
308 static BndSeg* getBndSegOfEdge( const TVDEdge* edge,
309 vector< vector< BndSeg > >& bndSegsPerEdge )
314 size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( edge );
315 size_t oppEdgeIndex = SMESH_MAT2d::Branch::getGeomEdge ( edge );
316 if ( oppEdgeIndex < bndSegsPerEdge.size() &&
317 oppSegIndex < bndSegsPerEdge[ oppEdgeIndex ].size() )
319 seg = & bndSegsPerEdge[ oppEdgeIndex ][ oppSegIndex ];
325 void setBranch( int branchID, vector< vector< BndSeg > >& bndSegsPerEdge )
327 _branchID = branchID;
329 // pass branch to an opposite BndSeg
331 if ( BndSeg* oppSeg = getBndSegOfEdge( _edge->twin(), bndSegsPerEdge ))
333 if ( oppSeg->_branchID == theNoBrachID )
334 oppSeg->_branchID = -branchID;
337 bool hasOppositeEdge()
339 if ( !_edge ) return false;
340 return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != theNoEdgeID );
343 // check a next segment in CW order
344 bool isSameBranch( const BndSeg& seg2 )
346 if ( !_edge || !seg2._edge )
349 if ( _edge->twin() == seg2._edge )
352 const TVDCell* cell1 = this->_edge->twin()->cell();
353 const TVDCell* cell2 = seg2. _edge->twin()->cell();
354 if ( cell1 == cell2 )
357 const TVDEdge* edgeMedium1 = this->_edge->twin()->next();
358 const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev();
360 if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() )
362 if ( edgeMedium1->twin() == edgeMedium2 )
364 // edgeMedium's are edges whose twin()->cell is built on an end point of inSegment
365 // and is located between cell1 and cell2
366 if ( edgeMedium1->twin() == edgeMedium2->twin() ) // is this possible???
368 if ( edgeMedium1->twin() == edgeMedium2->twin()->next() &&
369 edgeMedium1->twin()->cell()->contains_point() )
372 else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
374 if ( edgeMedium1->twin() == edgeMedium2 &&
375 SMESH_MAT2d::Branch::getGeomEdge( edgeMedium1 ) ==
376 SMESH_MAT2d::Branch::getGeomEdge( edgeMedium2 ))
377 // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
385 // -------------------------------------------------------------------------------------
389 struct BranchIterator
392 const std::vector<const TVDEdge*> & _edges;
395 BranchIterator(const std::vector<const TVDEdge*> & edges, int i )
396 :_i( i ), _size( edges.size() ), _edges( edges )
398 _closed = ( edges[0]->vertex1() == edges.back()->vertex0() || // closed branch
399 edges[0]->vertex0() == edges.back()->vertex1() );
401 const TVDEdge* operator++() { ++_i; return edge(); }
402 const TVDEdge* operator--() { --_i; return edge(); }
403 bool operator<( const BranchIterator& other ) { return _i < other._i; }
404 BranchIterator& operator=( const BranchIterator& other ) { _i = other._i; return *this; }
405 void set(int i) { _i = i; }
406 int index() const { return _i; }
407 int indexMod() const { return ( _i + _size ) % _size; }
408 const TVDEdge* edge() const {
409 return _closed ? _edges[ indexMod() ] : ( _i < 0 || _i >= _size ) ? 0 : _edges[ _i ];
411 const TVDEdge* edgePrev() { --_i; const TVDEdge* e = edge(); ++_i; return e; }
412 const TVDEdge* edgeNext() { ++_i; const TVDEdge* e = edge(); --_i; return e; }
415 //================================================================================
417 * \brief debug: to visually check found MA edges
419 //================================================================================
421 void bndSegsToMesh( const vector< vector< BndSeg > >& bndSegsPerEdge )
424 if ( !getenv("bndSegsToMesh")) return;
425 map< const TVDVertex *, int > v2Node;
426 map< const TVDVertex *, int >::iterator v2n;
427 set< const TVDEdge* > addedEdges;
429 const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
430 SMESH_File file(fileName, false );
432 file.openForWriting();
434 text << "import salome, SMESH\n";
435 text << "salome.salome_init()\n";
436 text << "from salome.smesh import smeshBuilder\n";
437 text << "smesh = smeshBuilder.New(salome.myStudy)\n";
438 text << "m=smesh.Mesh()\n";
439 for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
441 const vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
442 for ( size_t i = 0; i < bndSegs.size(); ++i )
444 if ( !bndSegs[i]._edge )
445 text << "# E=" << iE << " i=" << i << " NULL edge\n";
446 else if ( !bndSegs[i]._edge->vertex0() ||
447 !bndSegs[i]._edge->vertex1() )
448 text << "# E=" << iE << " i=" << i << " INFINITE edge\n";
449 else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
450 addedEdges.insert( bndSegs[i]._edge->twin() ).second )
452 v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
453 size_t n0 = v2n->second;
454 if ( n0 == v2Node.size() )
455 text << "n" << n0 << " = m.AddNode( "
456 << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
457 << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
459 v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
460 size_t n1 = v2n->second;
461 if ( n1 == v2Node.size() )
462 text << "n" << n1 << " = m.AddNode( "
463 << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
464 << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
466 text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
471 file.write( text.c_str(), text.size() );
472 cout << "execfile( '" << fileName << "')" << endl;
476 //================================================================================
478 * \brief Computes length of a TVDEdge
480 //================================================================================
482 double length( const TVDEdge* edge )
484 gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(),
485 edge->vertex0()->y() - edge->vertex1()->y() );
489 //================================================================================
491 * \brief Compute scale to have the same 2d proportions as in 3d
493 //================================================================================
495 void computeProportionScale( const TopoDS_Face& face,
496 const Bnd_B2d& uvBox,
499 scale[0] = scale[1] = 1.;
500 if ( uvBox.IsVoid() ) return;
503 Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
505 const int nbDiv = 30;
506 gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
507 gp_XY uvMid = 0.5 * ( uvMin + uvMax );
508 double du = ( uvMax.X() - uvMin.X() ) / nbDiv;
509 double dv = ( uvMax.Y() - uvMin.Y() ) / nbDiv;
511 double uLen3d = 0, vLen3d = 0;
512 gp_Pnt uPrevP = surface->Value( uvMin.X(), uvMid.Y() );
513 gp_Pnt vPrevP = surface->Value( uvMid.X(), uvMin.Y() );
514 for (int i = 1; i <= nbDiv; i++)
516 double u = uvMin.X() + du * i;
517 double v = uvMin.Y() + dv * i;
518 gp_Pnt uP = surface->Value( u, uvMid.Y() );
519 gp_Pnt vP = surface->Value( uvMid.X(), v );
520 uLen3d += uP.Distance( uPrevP );
521 vLen3d += vP.Distance( vPrevP );
525 scale[0] = uLen3d / ( uvMax.X() - uvMin.X() );
526 scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() );
529 //================================================================================
531 * \brief Fill input data for construct_voronoi()
533 //================================================================================
535 bool makeInputData(const TopoDS_Face& face,
536 const std::vector< TopoDS_Edge >& edges,
537 const double minSegLen,
538 vector< InPoint >& inPoints,
539 vector< InSegment>& inSegments,
542 const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization
545 // discretize the EDGEs to get 2d points and segments
547 vector< vector< UVU > > uvuVec( edges.size() );
549 for ( size_t iE = 0; iE < edges.size(); ++iE )
551 vector< UVU > & points = uvuVec[ iE ];
554 Handle(Geom_Curve) c3d = BRep_Tool::Curve ( edges[ iE ], loc, f, l );
555 Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface( edges[ iE ], face, f, l );
556 if ( c2d.IsNull() ) return false;
558 points.push_back( UVU( c2d->Value( f ), f ));
559 uvBox.Add( points.back()._uv );
561 Geom2dAdaptor_Curve c2dAdaptor (c2d, f,l );
562 double curDeflect = 0.3; //0.01; //Curvature deflection
563 double angDeflect = 0.2; // 0.09; //Angular deflection
565 GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect);
566 // if ( discret.NbPoints() > 2 )
571 // discret.Initialize( c2dAdaptor, 100, curDeflect );
572 // cout << "C " << curDeflect << " " << discret.NbPoints() << endl;
573 // curDeflect *= 1.5;
575 // while ( discret.NbPoints() > 5 );
579 // discret.Initialize( c2dAdaptor, angDeflect, 100 );
580 // cout << "A " << angDeflect << " " << discret.NbPoints() << endl;
581 // angDeflect *= 1.5;
583 // while ( discret.NbPoints() > 5 );
587 pPrev = c3d->Value( f );
588 if ( discret.NbPoints() > 2 )
589 for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
591 double u = discret.Parameter(i);
595 int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
596 double dU = ( u - points.back()._u ) / nbDiv;
597 for ( int iD = 1; iD < nbDiv; ++iD )
599 double uD = points.back()._u + dU;
600 points.push_back( UVU( c2d->Value( uD ), uD ));
604 points.push_back( UVU( c2d->Value( u ), u ));
605 uvBox.Add( points.back()._uv );
607 // if ( !c3d.IsNull() )
609 // vector<double> params;
610 // GeomAdaptor_Curve c3dAdaptor( c3d,f,l );
613 // const double deflection = minSegLen * 0.1;
614 // GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true );
615 // if ( !discret.IsDone() )
617 // int nbP = discret.NbPoints();
618 // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
619 // params.push_back( discret.Parameter(i) );
623 // double eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor );
624 // int nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef ));
625 // double segLen = eLen / nbSeg;
626 // GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l );
627 // int nbP = Min( discret.NbPoints(), nbSeg + 1 );
628 // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
629 // params.push_back( discret.Parameter(i) );
631 // for ( size_t i = 0; i < params.size(); ++i )
633 // points.push_back( UVU( c2d->Value( params[i] ), params[i] ));
634 // uvBox.Add( points.back()._uv );
637 if ( points.size() < 2 )
639 points.push_back( UVU( c2d->Value( l ), l ));
640 uvBox.Add( points.back()._uv );
642 if ( edges[ iE ].Orientation() == TopAbs_REVERSED )
643 std::reverse( points.begin(), points.end() );
646 // make connected EDGEs have same UV at shared VERTEX
647 TopoDS_Vertex vShared;
648 for ( size_t iE = 0; iE < edges.size(); ++iE )
650 size_t iE2 = (iE+1) % edges.size();
651 if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared ))
653 if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
655 vector< UVU > & points1 = uvuVec[ iE ];
656 vector< UVU > & points2 = uvuVec[ iE2 ];
657 gp_Pnt2d & uv1 = points1.back() ._uv;
658 gp_Pnt2d & uv2 = points2.front()._uv;
659 uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() );
662 // get scale to have the same 2d proportions as in 3d
663 computeProportionScale( face, uvBox, scale );
665 // make 'scale' such that to have coordinates precise enough when converted to int
667 gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
668 uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
669 uvMin.ChangeCoord(2) = uvMin.Y() * scale[1];
670 uvMax.ChangeCoord(1) = uvMax.X() * scale[0];
671 uvMax.ChangeCoord(2) = uvMax.Y() * scale[1];
672 double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
673 Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
674 int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
675 const double precision = Min( 1e-5, minSegLen * 1e-2 );
676 double preciScale = Min( vMax[iMax] / precision,
677 std::numeric_limits<int>::max() / vMax[iMax] );
678 preciScale /= scale[iMax];
679 double roundedScale = 10; // to ease debug
680 while ( roundedScale * 10 < preciScale )
682 scale[0] *= roundedScale;
683 scale[1] *= roundedScale;
685 // create input points and segments
690 for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
691 nbPnt += uvuVec[ iE ].size();
692 inPoints.resize( nbPnt );
693 inSegments.reserve( nbPnt );
696 if ( face.Orientation() == TopAbs_REVERSED )
698 for ( int iE = uvuVec.size()-1; iE >= 0; --iE )
700 vector< UVU > & points = uvuVec[ iE ];
701 inPoints[ iP++ ] = points.back().getInPoint( scale );
702 for ( size_t i = points.size()-1; i >= 1; --i )
704 inPoints[ iP++ ] = points[i-1].getInPoint( scale );
705 inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
706 if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
707 return false; // too short segment
713 for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
715 vector< UVU > & points = uvuVec[ iE ];
716 inPoints[ iP++ ] = points[0].getInPoint( scale );
717 for ( size_t i = 1; i < points.size(); ++i )
719 inPoints[ iP++ ] = points[i].getInPoint( scale );
720 inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
721 if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
722 return false; // too short segment
727 theScale[0] = scale[0];
728 theScale[1] = scale[1];
733 //================================================================================
735 * \brief Update a branch joined to another one
737 //================================================================================
739 void updateJoinedBranch( vector< const TVDEdge* > & branchEdges,
741 vector< vector< BndSeg > > & bndSegs,
747 for ( size_t i = 0; i < branchEdges.size(); ++i )
749 if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i], bndSegs )) &&
750 ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
752 seg1->_branchID /= seg1->branchID();
753 seg2->_branchID /= seg2->branchID();
754 seg1->_branchID *= -newID;
755 seg2->_branchID *= -newID;
756 branchEdges[i] = branchEdges[i]->twin();
759 std::reverse( branchEdges.begin(), branchEdges.end() );
763 for ( size_t i = 0; i < branchEdges.size(); ++i )
765 if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i], bndSegs )) &&
766 ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
768 seg1->_branchID /= seg1->branchID();
769 seg2->_branchID /= seg2->branchID();
770 seg1->_branchID *= newID;
771 seg2->_branchID *= newID;
777 //================================================================================
779 * \brief Create MA branches and FACE boundary data
780 * \param [in] vd - voronoi diagram of \a inSegments
781 * \param [in] inPoints - FACE boundary points
782 * \param [in,out] inSegments - FACE boundary segments
783 * \param [out] branch - MA branches to fill
784 * \param [out] branchEnd - ends of MA branches to fill
785 * \param [out] boundary - FACE boundary to fill
787 //================================================================================
789 void makeMA( const TVD& vd,
790 const bool ignoreCorners,
791 vector< InPoint >& inPoints,
792 vector< InSegment > & inSegments,
793 vector< SMESH_MAT2d::Branch >& branch,
794 vector< const SMESH_MAT2d::BranchEnd* >& branchPnt,
795 SMESH_MAT2d::Boundary& boundary )
797 // Associate MA cells with geom EDGEs
798 for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
800 const TVDCell* cell = &(*it);
801 if ( cell->contains_segment() )
803 InSegment& seg = inSegments[ cell->source_index() ];
805 seg.setGeomEdgeToCell( cell, seg._geomEdgeInd );
809 InSegment::setGeomEdgeToCell( cell, theNoEdgeID );
813 vector< bool > inPntChecked( inPoints.size(), false );
815 // Find MA edges of each inSegment
817 for ( size_t i = 0; i < inSegments.size(); ++i )
819 InSegment& inSeg = inSegments[i];
821 // get edges around the cell lying on MA
822 bool hasSecondary = false;
823 const TVDEdge* edge = inSeg._cell->incident_edge();
825 edge = edge->next(); // Returns the CCW next edge within the cell.
826 if ( edge->is_primary() && !inSeg.isExternal( edge ) )
827 inSeg._edges.push_back( edge ); // edge equidistant from two InSegments
830 } while (edge != inSeg._cell->incident_edge());
832 // there can be several continuous MA edges but maEdges can begin in the middle of
833 // a chain of continuous MA edges. Make the chain continuous.
834 list< const TVDEdge* >& maEdges = inSeg._edges;
835 if ( maEdges.empty() )
838 while ( maEdges.back()->next() == maEdges.front() )
839 maEdges.splice( maEdges.end(), maEdges, maEdges.begin() );
841 // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE
842 list< const TVDEdge* >::iterator e = maEdges.begin();
843 while ( e != maEdges.end() )
845 const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge
846 size_t geoE2 = InSegment::getGeomEdge( cell2 );
847 bool toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e ));
849 e = maEdges.erase( e );
853 if ( maEdges.empty() )
856 // add MA edges corresponding to concave InPoints
857 for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg
859 InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 );
860 size_t pInd = inPnt.index( inPoints );
861 if ( inPntChecked[ pInd ] )
864 inPntChecked[ pInd-1 ] &&
865 inPoints[ pInd-1 ] == inPnt )
867 inPntChecked[ pInd ] = true;
869 const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
870 if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
872 const TVDEdge* edge = // a secondary TVDEdge connecting inPnt and maE
873 is2nd ? maE->prev() : maE->next();
874 while ( inSeg.isConnected( edge ))
876 if ( edge->is_primary() ) break; // this should not happen
877 const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
878 if ( inSeg.getGeomEdge( edge2->cell() ) != theNoEdgeID )
879 break; // cell of an InSegment
880 bool hasInfinite = false;
881 list< const TVDEdge* > pointEdges;
885 edge = edge->next(); // Returns the CCW next edge within the cell.
886 if ( edge->is_infinite() )
888 else if ( edge->is_primary() && !inSeg.isExternal( edge ))
889 pointEdges.push_back( edge );
891 while ( edge != edge2 && !hasInfinite );
893 if ( hasInfinite || pointEdges.empty() )
895 inPnt._edges.splice( inPnt._edges.end(), pointEdges );
896 inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd );
898 edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next();
900 } // add MA edges corresponding to concave InPoints
902 } // loop on inSegments to find corresponding MA edges
905 // -------------------------------------------
906 // Create Branches and BndPoints for each EDGE
907 // -------------------------------------------
909 if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ )
911 inPntChecked[0] = false; // do not use the 1st point twice
912 //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), theNoEdgeID );
913 inPoints[0]._edges.clear();
916 // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
918 vector< vector< BndSeg > > bndSegsPerEdge( boundary.nbEdges() ); // all BndSeg's
920 vector< BndSeg > bndSegs; // bndSeg's of a current EDGE
921 size_t prevGeomEdge = theNoEdgeID;
923 list< const TVDEdge* >::reverse_iterator e;
924 for ( size_t i = 0; i < inSegments.size(); ++i )
926 InSegment& inSeg = inSegments[i];
928 if ( inSeg._geomEdgeInd != prevGeomEdge )
930 if ( !bndSegs.empty() )
931 bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
932 prevGeomEdge = inSeg._geomEdgeInd;
935 // segments around 1st concave point
936 size_t ip0 = inSeg._p0->index( inPoints );
937 if ( inPntChecked[ ip0 ] )
938 for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
939 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
940 inPntChecked[ ip0 ] = false;
942 // segments of InSegment's
943 const size_t nbMaEdges = inSeg._edges.size();
944 switch ( nbMaEdges ) {
945 case 0: // "around" circle center
946 bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
948 bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
950 gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
951 inSeg._p1->_b - inSeg._p0->_b );
952 const double inSegLen2 = inSegDir.SquareModulus();
953 e = inSeg._edges.rbegin();
954 for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
956 gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
957 (*e)->vertex0()->y() - inSeg._p0->_b );
958 double r = toMA * inSegDir / inSegLen2;
959 double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
960 bndSegs.push_back( BndSeg( &inSeg, *e, u ));
962 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
964 // segments around 2nd concave point
965 size_t ip1 = inSeg._p1->index( inPoints );
966 if ( inPntChecked[ ip1 ] )
967 for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
968 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
969 inPntChecked[ ip1 ] = false;
971 if ( !bndSegs.empty() )
972 bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
975 // prepare to MA branch search
976 for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
978 // 1) make TVDEdge's know it's BndSeg to enable passing branchID to
979 // an opposite BndSeg in BndSeg::setBranch(); geom EDGE ID is known from TVDCell
980 // 2) connect bndSegs via BndSeg::_prev
982 vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
983 if ( bndSegs.empty() ) continue;
985 for ( size_t i = 1; i < bndSegs.size(); ++i )
987 bndSegs[i]._prev = & bndSegs[i-1];
988 bndSegs[i].setIndexToEdge( i );
990 // look for the last bndSeg of previous EDGE to set bndSegs[0]._prev
991 const InPoint& p0 = bndSegs[0]._inSeg->point0();
992 for ( size_t iE2 = 0; iE2 < bndSegsPerEdge.size(); ++iE2 )
993 if ( p0 == bndSegsPerEdge[ iE2 ].back()._inSeg->point1() )
995 bndSegs[0]._prev = & bndSegsPerEdge[ iE2 ].back();
998 bndSegs[0].setIndexToEdge( 0 );
1001 bndSegsToMesh( bndSegsPerEdge ); // debug: visually check found MA edges
1004 // Find TVDEdge's of Branches and associate them with bndSegs
1006 vector< vector<const TVDEdge*> > branchEdges;
1007 branchEdges.reserve( boundary.nbEdges() * 4 );
1009 map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
1011 int branchID = 1; // we code orientation as branchID sign
1012 branchEdges.resize( branchID );
1014 for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
1016 vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
1017 for ( size_t i = 0; i < bndSegs.size(); ++i )
1019 if ( bndSegs[i].branchID() )
1021 if ( bndSegs[i]._prev &&
1022 bndSegs[i]._branchID == -bndSegs[i]._prev->_branchID &&
1025 SMESH_MAT2d::BranchEndType type =
1026 ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
1027 SMESH_MAT2d::BE_ON_VERTEX :
1028 SMESH_MAT2d::BE_END );
1029 endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
1033 if ( !bndSegs[i]._prev &&
1034 !bndSegs[i].hasOppositeEdge() )
1037 if ( !bndSegs[i]._prev ||
1038 !bndSegs[i]._prev->isSameBranch( bndSegs[i] ))
1040 branchEdges.resize(( branchID = branchEdges.size()) + 1 );
1041 if ( bndSegs[i]._edge && bndSegs[i]._prev )
1042 endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_BRANCH_POINT ));
1044 else if ( bndSegs[i]._prev->_branchID )
1046 branchID = bndSegs[i]._prev->_branchID; // with sign
1048 else if ( bndSegs[i]._edge ) // 1st bndSeg of a WIRE
1050 branchEdges.resize(( branchID = branchEdges.size()) + 1 );
1051 if ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ))
1053 if ( bndSegs[i]._inSeg->point0() == bndSegs[i]._edge->vertex1() )
1054 endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_ON_VERTEX ));
1056 endType.insert( make_pair( bndSegs[i]._edge->vertex0(), SMESH_MAT2d::BE_ON_VERTEX ));
1060 bndSegs[i].setBranch( branchID, bndSegsPerEdge ); // set to i-th and to the opposite bndSeg
1061 if ( bndSegs[i].hasOppositeEdge() )
1062 branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
1066 // join the 1st and the last branch edges if it is the same branch
1067 // if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
1068 // bndSegs.back().isSameBranch( bndSegs.front() ))
1070 // vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
1071 // vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID() ];
1072 // br1.insert( br1.begin(), br2.begin(), br2.end() );
1076 // remove branches ending at BE_ON_VERTEX
1078 vector<bool> isBranchRemoved( branchEdges.size(), false );
1080 if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
1082 // find branches to remove
1083 map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et;
1084 for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
1086 if ( branchEdges[iB].empty() )
1088 const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
1089 const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
1090 v2et = endType.find( v0 );
1091 if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
1092 isBranchRemoved[ iB ] = true;
1093 v2et = endType.find( v1 );
1094 if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
1095 isBranchRemoved[ iB ] = true;
1097 // try to join not removed branches into one
1098 for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
1100 if ( branchEdges[iB].empty() || isBranchRemoved[iB] )
1102 const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
1103 const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
1104 v2et = endType.find( v0 );
1105 if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
1107 v2et = endType.find( v1 );
1108 if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
1113 for ( int isV0 = 0; isV0 < 2; ++isV0 )
1115 const TVDVertex* v = isV0 ? v0 : v1;
1116 size_t iBrToJoin = 0;
1117 for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
1119 if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
1121 const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
1122 const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
1123 if ( v == v02 || v == v12 )
1125 if ( iBrToJoin > 0 )
1128 break; // more than 2 not removed branches meat at a TVDVertex
1133 if ( iBrToJoin > 0 )
1135 vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
1136 const TVDVertex* v02 = branch[0]->vertex1();
1137 const TVDVertex* v12 = branch.back()->vertex0();
1138 updateJoinedBranch( branch, iB, bndSegsPerEdge, /*reverse=*/(v0 == v02 || v1 == v12 ));
1139 if ( v0 == v02 || v0 == v12 )
1140 branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
1142 branchEdges[iB].insert( branchEdges[iB].end(), branch.begin(), branch.end() );
1146 } // loop on branchEdges
1147 } // if ( ignoreCorners )
1149 // associate branchIDs and the input branch vector (arg)
1150 vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 );
1152 for ( size_t i = 0; i < branchEdges.size(); ++i )
1154 nbBranches += ( !branchEdges[i].empty() );
1156 branch.resize( nbBranches );
1158 for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed
1160 if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ])
1161 branchByID[ brID ] = & branch[ iBr++ ];
1163 for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed
1165 if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ])
1166 branchByID[ brID ] = & branch[ iBr++ ];
1169 // Fill in BndPoints of each EDGE of the boundary
1172 int edgeInd = -1, dInd = 0;
1173 for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
1175 vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
1176 SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( iE );
1178 // make TVDEdge know an index of bndSegs within BndPoints
1179 for ( size_t i = 0; i < bndSegs.size(); ++i )
1180 if ( bndSegs[i]._edge )
1181 SMESH_MAT2d::Branch::setBndSegment( i, bndSegs[i]._edge );
1183 // parameters on EDGE
1185 bndPoints._params.reserve( bndSegs.size() + 1 );
1186 bndPoints._params.push_back( bndSegs[ 0 ]._inSeg->_p0->_param );
1188 for ( size_t i = 0; i < bndSegs.size(); ++i )
1189 bndPoints._params.push_back( bndSegs[ i ]._uLast );
1193 bndPoints._maEdges.reserve( bndSegs.size() );
1195 for ( size_t i = 0; i < bndSegs.size(); ++i )
1197 const size_t brID = bndSegs[ i ].branchID();
1198 const SMESH_MAT2d::Branch* br = branchByID[ brID ];
1200 if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
1204 if (( edgeInd < 0 ||
1205 edgeInd >= (int) branchEdges[ brID ].size() ) ||
1206 ( branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge &&
1207 branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge ))
1209 if ( bndSegs[ i ]._branchID < 0 )
1212 for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd )
1213 if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge )
1216 else // bndSegs[ i ]._branchID > 0
1219 for ( edgeInd = 0; edgeInd < (int)branchEdges[ brID ].size(); ++edgeInd )
1220 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
1227 // no MA edge, bndSeg corresponds to an end point of a branch
1228 if ( bndPoints._maEdges.empty() )
1231 edgeInd = branchEdges[ brID ].size();
1232 dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
1235 bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
1237 } // loop on bndSegs of an EDGE
1238 } // loop on all bndSegs to construct Boundary
1240 // Initialize branches
1242 // find a not removed branch
1243 size_t iBrNorRemoved = 0;
1244 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1245 if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] )
1247 iBrNorRemoved = brID;
1250 // fill the branches with MA edges
1251 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1252 if ( !branchEdges[brID].empty() )
1254 branchByID[ brID ]->init( branchEdges[brID], & boundary, endType );
1256 // mark removed branches
1257 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1258 if ( isBranchRemoved[brID] && iBrNorRemoved > 0 )
1260 SMESH_MAT2d::Branch* branch = branchByID[ brID ];
1261 SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ];
1262 bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT );
1263 const TVDVertex* branchVextex =
1264 is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex;
1265 SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex );
1266 branch->setRemoved( bp );
1268 // set branches to branch ends
1269 for ( size_t i = 0; i < branch.size(); ++i )
1270 if ( !branch[i].isRemoved() )
1271 branch[i].setBranchesToEnds( branch );
1273 // fill branchPnt arg
1274 map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
1275 for ( size_t i = 0; i < branch.size(); ++i )
1277 if ( branch[i].getEnd(0)->_branches.size() > 2 )
1278 v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
1279 if ( branch[i].getEnd(1)->_branches.size() > 2 )
1280 v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
1282 branchPnt.resize( v2end.size() );
1283 map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
1284 for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
1285 branchPnt[ i ] = v2e->second;
1291 //================================================================================
1293 * \brief MedialAxis constructor
1294 * \param [in] face - a face to create MA for
1295 * \param [in] edges - edges of the face (possibly not all) on the order they
1296 * encounter in the face boundary.
1297 * \param [in] minSegLen - minimal length of a mesh segment used to discretize
1298 * the edges. It is used to define precision of MA approximation
1300 //================================================================================
1302 SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face,
1303 const std::vector< TopoDS_Edge >& edges,
1304 const double minSegLen,
1305 const bool ignoreCorners):
1306 _face( face ), _boundary( edges.size() )
1308 // input to construct_voronoi()
1309 vector< InPoint > inPoints;
1310 vector< InSegment> inSegments;
1311 if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
1314 inSegmentsToFile( inSegments );
1316 // build voronoi diagram
1317 construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
1320 makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary );
1322 // count valid branches
1323 _nbBranches = _branch.size();
1324 for ( size_t i = 0; i < _branch.size(); ++i )
1325 if ( _branch[i].isRemoved() )
1329 //================================================================================
1331 * \brief Returns the i-th branch
1333 //================================================================================
1335 const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const
1337 return i < _nbBranches ? &_branch[i] : 0;
1340 //================================================================================
1342 * \brief Return UVs of ends of MA edges of a branch
1344 //================================================================================
1346 void SMESH_MAT2d::MedialAxis::getPoints( const Branch* branch,
1347 std::vector< gp_XY >& points) const
1349 branch->getPoints( points, _scale );
1352 //================================================================================
1354 * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
1355 * \param [in] iEdge - index of geom EDGE within a vector passed at MA construction
1356 * \param [in] u - parameter of the point on EDGE curve
1357 * \param [out] p - the found BranchPoint
1358 * \return bool - is OK
1360 //================================================================================
1362 bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
1364 BranchPoint& p ) const
1366 if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1369 const BndPoints& points = _pointsPerEdge[ iEdge ];
1370 const bool edgeReverse = ( points._params[0] > points._params.back() );
1372 if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
1373 u = edgeReverse ? points._params.back() : points._params[0];
1374 else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
1375 u = edgeReverse ? points._params[0] : points._params.back();
1377 double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
1378 int i = int( r * double( points._maEdges.size()-1 ));
1381 while ( points._params[i ] < u ) --i;
1382 while ( points._params[i+1] > u ) ++i;
1386 while ( points._params[i ] > u ) --i;
1387 while ( points._params[i+1] < u ) ++i;
1390 if ( points._params[i] == points._params[i+1] ) // coincident points at some end
1392 int di = ( points._params[0] == points._params[i] ) ? +1 : -1;
1393 while ( points._params[i] == points._params[i+1] )
1395 if ( i < 0 || i+1 >= (int)points._params.size() )
1399 double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
1401 if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
1403 if ( i < (int)points._maEdges.size() / 2 ) // near 1st point
1405 while ( i < (int)points._maEdges.size()-1 && !points._maEdges[ i ].second )
1407 edgeParam = edgeReverse;
1409 else // near last point
1411 while ( i > 0 && !points._maEdges[ i ].second )
1413 edgeParam = !edgeReverse;
1416 const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
1417 bool maReverse = ( maE.second < 0 );
1419 p._branch = maE.first;
1420 p._iEdge = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
1421 p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam;
1426 //================================================================================
1428 * \brief Returns a BranchPoint corresponding to a given BoundaryPoint on a geom EDGE
1429 * \param [in] bp - the BoundaryPoint
1430 * \param [out] p - the found BranchPoint
1431 * \return bool - is OK
1433 //================================================================================
1435 bool SMESH_MAT2d::Boundary::getBranchPoint( const BoundaryPoint& bp,
1436 BranchPoint& p ) const
1438 return getBranchPoint( bp._edgeIndex, bp._param, p );
1441 //================================================================================
1443 * \brief Check if a given boundary segment is a null-length segment on a concave
1445 * \param [in] iEdge - index of a geom EDGE
1446 * \param [in] iSeg - index of a boundary segment
1447 * \return bool - true if the segment is on concave corner
1449 //================================================================================
1451 bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
1453 if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1456 const BndPoints& points = _pointsPerEdge[ iEdge ];
1457 if ( points._params.size() <= iSeg+1 )
1460 return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20;
1463 //================================================================================
1465 * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end
1467 //================================================================================
1469 bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
1471 if ( bp._edgeIndex >= _pointsPerEdge.size() )
1474 const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
1475 if ( Abs( bp._param - points._params[0]) < Abs( points._params.back() - bp._param ))
1476 bp._param = points._params[0];
1478 bp._param = points._params.back();
1483 //================================================================================
1485 * \brief Creates a 3d curve corresponding to a Branch
1486 * \param [in] branch - the Branch
1487 * \return Adaptor3d_Curve* - the new curve the caller is to delete
1489 //================================================================================
1491 Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
1493 Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
1494 if ( surface.IsNull() )
1498 branch.getPoints( uv, _scale );
1499 if ( uv.size() < 2 )
1502 vector< TopoDS_Vertex > vertex( uv.size() );
1503 for ( size_t i = 0; i < uv.size(); ++i )
1504 vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
1507 BRep_Builder aBuilder;
1508 aBuilder.MakeWire(aWire);
1509 for ( size_t i = 1; i < vertex.size(); ++i )
1511 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
1512 aBuilder.Add( aWire, edge );
1515 // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
1516 // aWire.Closed(true); // issue 0021141
1518 return new BRepAdaptor_CompCurve( aWire );
1521 //================================================================================
1523 * \brief Copy points of an EDGE
1525 //================================================================================
1527 void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>& maEdges,
1528 const Boundary* boundary,
1529 map< const TVDVertex*, BranchEndType >& endType )
1531 if ( maEdges.empty() ) return;
1533 _boundary = boundary;
1534 _maEdges.swap( maEdges );
1537 _params.reserve( _maEdges.size() + 1 );
1538 _params.push_back( 0. );
1539 for ( size_t i = 0; i < _maEdges.size(); ++i )
1540 _params.push_back( _params.back() + length( _maEdges[i] ));
1542 for ( size_t i = 1; i < _params.size(); ++i )
1543 _params[i] /= _params.back();
1546 _endPoint1._vertex = _maEdges.front()->vertex1();
1547 _endPoint2._vertex = _maEdges.back ()->vertex0();
1549 if ( endType.count( _endPoint1._vertex ))
1550 _endPoint1._type = endType[ _endPoint1._vertex ];
1551 if ( endType.count( _endPoint2._vertex ))
1552 _endPoint2._type = endType[ _endPoint2._vertex ];
1555 //================================================================================
1557 * \brief fills BranchEnd::_branches of its ends
1559 //================================================================================
1561 void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
1563 for ( size_t i = 0; i < branches.size(); ++i )
1565 if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
1566 this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
1567 this->_endPoint1._branches.push_back( &branches[i] );
1569 if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
1570 this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
1571 this->_endPoint2._branches.push_back( &branches[i] );
1575 //================================================================================
1577 * \brief returns a BranchPoint corresponding to a TVDVertex
1579 //================================================================================
1581 SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const
1587 if ( vertex == _maEdges[0]->vertex1() )
1593 for ( ; p._iEdge < _maEdges.size(); ++p._iEdge )
1594 if ( vertex == _maEdges[ p._iEdge ]->vertex0() )
1596 p._edgeParam = _params[ p._iEdge ];
1603 //================================================================================
1605 * \brief Sets a proxy point for a removed branch
1606 * \param [in] proxyPoint - a point of another branch to which all points of this
1609 //================================================================================
1611 void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint )
1613 _proxyPoint = proxyPoint;
1616 //================================================================================
1618 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1619 * \param [in] param - [0;1] normalized param on the Branch
1620 * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1621 * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1622 * \return bool - true if the BoundaryPoint's found
1624 //================================================================================
1626 bool SMESH_MAT2d::Branch::getBoundaryPoints(double param,
1628 BoundaryPoint& bp2 ) const
1630 if ( param < _params[0] || param > _params.back() )
1633 // look for an index of a MA edge by param
1634 double ip = param * _params.size();
1635 size_t i = size_t( Min( int( _maEdges.size()-1), int( ip )));
1637 while ( param < _params[i ] ) --i;
1638 while ( param > _params[i+1] ) ++i;
1640 double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
1642 return getBoundaryPoints( i, r, bp1, bp2 );
1645 //================================================================================
1647 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1648 * \param [in] iMAEdge - index of a MA edge within this Branch
1649 * \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
1650 * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1651 * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1652 * \return bool - true if the BoundaryPoint's found
1654 //================================================================================
1656 bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t iMAEdge,
1659 BoundaryPoint& bp2 ) const
1662 return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 );
1664 if ( iMAEdge > _maEdges.size() )
1666 if ( iMAEdge == _maEdges.size() )
1667 iMAEdge = _maEdges.size() - 1;
1669 size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
1670 size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
1671 size_t iSeg1 = getBndSegment( _maEdges[ iMAEdge ] );
1672 size_t iSeg2 = getBndSegment( _maEdges[ iMAEdge ]->twin() );
1674 return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
1675 _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
1678 //================================================================================
1680 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1682 //================================================================================
1684 bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
1686 BoundaryPoint& bp2 ) const
1688 return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
1691 //================================================================================
1693 * \brief Return a parameter of a BranchPoint normalized within this Branch
1695 //================================================================================
1697 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
1699 if ( this != p._branch && p._branch )
1700 return p._branch->getParameter( p, u );
1703 return _proxyPoint._branch->getParameter( _proxyPoint, u );
1705 if ( p._iEdge > _params.size()-1 )
1707 if ( p._iEdge == _params.size()-1 )
1710 u = ( _params[ p._iEdge ] * ( 1 - p._edgeParam ) +
1711 _params[ p._iEdge+1 ] * p._edgeParam );
1716 //================================================================================
1718 * \brief Check type of both ends
1720 //================================================================================
1722 bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
1724 return ( _endPoint1._type == type || _endPoint2._type == type );
1727 //================================================================================
1729 * \brief Returns MA points
1730 * \param [out] points - the 2d points
1731 * \param [in] scale - the scale that was used to scale the 2d space of MA
1733 //================================================================================
1735 void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
1736 const double scale[2]) const
1738 points.resize( _maEdges.size() + 1 );
1740 points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
1741 _maEdges[0]->vertex1()->y() / scale[1] );
1743 for ( size_t i = 0; i < _maEdges.size(); ++i )
1744 points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
1745 _maEdges[i]->vertex0()->y() / scale[1] );
1748 //================================================================================
1750 * \brief Return indices of EDGEs equidistant from this branch
1752 //================================================================================
1754 void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
1755 std::vector< std::size_t >& edgeIDs2 ) const
1757 edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1758 edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1760 for ( size_t i = 1; i < _maEdges.size(); ++i )
1762 size_t ie1 = getGeomEdge( _maEdges[i] );
1763 size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1765 if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
1766 if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
1770 //================================================================================
1772 * \brief Looks for a BranchPoint position around a concave VERTEX
1774 //================================================================================
1776 bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1,
1777 std::vector< std::size_t >& edgeIDs2,
1778 std::vector< BranchPoint >& divPoints,
1779 const vector<const TVDEdge*>& maEdges,
1780 const vector<const TVDEdge*>& maEdgesTwin,
1783 // if there is a concave vertex between EDGEs
1784 // then position of a dividing BranchPoint is undefined, it is somewhere
1785 // on an arc-shaped part of the Branch around the concave vertex.
1786 // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
1787 // of the arc if there is no opposite VERTEX.
1788 // All null-length segments around a VERTEX belong to one of EDGEs.
1790 BranchPoint divisionPnt;
1791 divisionPnt._branch = this;
1793 BranchIterator iCur( maEdges, i );
1795 size_t ie1 = getGeomEdge( maEdges [i] );
1796 size_t ie2 = getGeomEdge( maEdgesTwin[i] );
1798 size_t iSeg1 = getBndSegment( iCur.edgePrev() );
1799 size_t iSeg2 = getBndSegment( iCur.edge() );
1800 bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
1801 bool isConcaNext = _boundary->isConcaveSegment( ie1, iSeg2 );
1802 if ( !isConcaNext && !isConcaPrev )
1805 bool isConcaveV = false;
1808 BranchIterator iPrev( maEdges, i ), iNext( maEdges, i );
1810 if ( isConcaNext ) // all null-length segments follow
1812 // look for a VERTEX of the opposite EDGE
1813 // iNext - next after all null-length segments
1814 while (( maE = ++iNext ))
1816 iSeg2 = getBndSegment( maE );
1817 if ( !_boundary->isConcaveSegment( ie1, iSeg2 ))
1820 bool vertexFound = false;
1821 for ( ++iCur; iCur < iNext; ++iCur )
1823 ie2 = getGeomEdge( maEdgesTwin[ iCur.indexMod() ] );
1824 if ( ie2 != edgeIDs2.back() )
1826 // opposite VERTEX found
1827 divisionPnt._iEdge = iCur.indexMod();
1828 divisionPnt._edgeParam = 0;
1829 divPoints.push_back( divisionPnt );
1830 edgeIDs1.push_back( ie1 );
1831 edgeIDs2.push_back( ie2 );
1838 iPrev = iNext; // not to add a BP in the moddle
1839 i = iNext.indexMod();
1843 else if ( isConcaPrev )
1845 // all null-length segments passed, find their beginning
1846 while (( maE = iPrev.edgePrev() ))
1848 iSeg1 = getBndSegment( maE );
1849 if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
1856 if ( iPrev.index() < i-1 || iNext.index() > i )
1858 // no VERTEX on the opposite EDGE, put the Branch Point in the middle
1859 divisionPnt._iEdge = iPrev.indexMod();
1861 double par1 = _params[ iPrev.indexMod() ], par2 = _params[ iNext.indexMod() ];
1862 double midPar = 0.5 * ( par1 + par2 );
1863 for ( ; _params[ iPrev.indexMod() ] < midPar; ++iPrev )
1864 divisionPnt._iEdge = iPrev.indexMod();
1865 divisionPnt._edgeParam =
1866 ( _params[ iPrev.indexMod() ] - midPar ) /
1867 ( _params[ iPrev.indexMod() ] - _params[ divisionPnt._iEdge ] );
1868 divPoints.push_back( divisionPnt );
1875 //================================================================================
1877 * \brief Return indices of opposite parts of EDGEs equidistant from this branch
1878 * \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
1879 * \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
1880 * \param [out] divPoints - BranchPoint's located between two successive unique
1881 * pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
1882 * of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
1883 * than number of \a edgeIDs
1885 //================================================================================
1887 void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
1888 std::vector< std::size_t >& edgeIDs2,
1889 std::vector< BranchPoint >& divPoints) const
1895 std::vector<const TVDEdge*> twins( _maEdges.size() );
1896 for ( size_t i = 0; i < _maEdges.size(); ++i )
1897 twins[i] = _maEdges[i]->twin();
1899 BranchIterator maIter ( _maEdges, 0 );
1900 BranchIterator twIter ( twins, 0 );
1901 // size_t lastConcaE1 = _boundary.nbEdges();
1902 // size_t lastConcaE2 = _boundary.nbEdges();
1904 // if ( maIter._closed ) // closed branch
1906 // edgeIDs1.push_back( getGeomEdge( _maEdges.back() ));
1907 // edgeIDs2.push_back( getGeomEdge( _maEdges.back()->twin() ));
1911 edgeIDs1.push_back( getGeomEdge( maIter.edge() ));
1912 edgeIDs2.push_back( getGeomEdge( twIter.edge() ));
1915 BranchPoint divisionPnt;
1916 divisionPnt._branch = this;
1918 for ( ++maIter, ++twIter; maIter.index() < (int)_maEdges.size(); ++maIter, ++twIter )
1920 size_t ie1 = getGeomEdge( maIter.edge() );
1921 size_t ie2 = getGeomEdge( twIter.edge() );
1923 bool otherE1 = ( edgeIDs1.back() != ie1 );
1924 bool otherE2 = ( edgeIDs2.back() != ie2 );
1926 if ( !otherE1 && !otherE2 && maIter._closed )
1928 int iSegPrev1 = getBndSegment( maIter.edgePrev() );
1929 int iSegCur1 = getBndSegment( maIter.edge() );
1930 otherE1 = Abs( iSegPrev1 - iSegCur1 ) != 1;
1931 int iSegPrev2 = getBndSegment( twIter.edgePrev() );
1932 int iSegCur2 = getBndSegment( twIter.edge() );
1933 otherE2 = Abs( iSegPrev2 - iSegCur2 ) != 1;
1936 if ( otherE1 || otherE2 )
1938 bool isConcaveV = false;
1939 if ( otherE1 && !otherE2 )
1941 isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints,
1942 _maEdges, twins, maIter._i );
1944 if ( !otherE1 && otherE2 )
1946 isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints,
1947 twins, _maEdges, maIter._i );
1952 ie1 = getGeomEdge( maIter.edge() );
1953 ie2 = getGeomEdge( twIter.edge() );
1955 if ( !isConcaveV || otherE1 || otherE2 )
1957 edgeIDs1.push_back( ie1 );
1958 edgeIDs2.push_back( ie2 );
1960 if ( divPoints.size() < edgeIDs1.size() - 1 )
1962 divisionPnt._iEdge = maIter.index();
1963 divisionPnt._edgeParam = 0;
1964 divPoints.push_back( divisionPnt );
1967 } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
1968 } // loop on _maEdges
1971 //================================================================================
1973 * \brief Store data of boundary segments in TVDEdge
1975 //================================================================================
1977 void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
1979 if ( maEdge ) maEdge->cell()->color( geomIndex );
1981 std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
1983 return maEdge ? maEdge->cell()->color() : std::string::npos;
1985 void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
1987 if ( maEdge ) maEdge->color( segIndex );
1989 std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
1991 return maEdge ? maEdge->color() : std::string::npos;
1994 //================================================================================
1996 * \brief Returns a boundary point on a given EDGE
1997 * \param [in] iEdge - index of the EDGE within MedialAxis
1998 * \param [in] iSeg - index of a boundary segment within this Branch
1999 * \param [in] u - [0;1] normalized param within \a iSeg-th segment
2000 * \param [out] bp - the found BoundaryPoint
2001 * \return bool - true if the BoundaryPoint is found
2003 //================================================================================
2005 bool SMESH_MAT2d::Boundary::getPoint( std::size_t iEdge,
2008 BoundaryPoint& bp ) const
2010 if ( iEdge >= _pointsPerEdge.size() )
2012 if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
2015 // This method is called by Branch that can have an opposite orientation,
2016 // hence u is inverted depending on orientation coded as a sign of _maEdge index
2017 bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
2021 double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
2022 double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
2024 bp._param = p0 * ( 1. - u ) + p1 * u;
2025 bp._edgeIndex = iEdge;