1 // Copyright (C) 2007-2015 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>
50 #include "SMESH_File.hxx"
51 #include "SMESH_Comment.hxx"
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;
64 // Input data for construct_voronoi()
65 // -------------------------------------------------------------------------------------
69 int _a, _b; // coordinates
70 double _param; // param on EDGE
71 InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
72 InPoint() : _a(0), _b(0), _param(0) {}
75 list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order
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; }
80 // -------------------------------------------------------------------------------------
88 size_t _geomEdgeInd; // EDGE index within the FACE
90 list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
92 InSegment( InPoint * p0, InPoint * p1, size_t iE)
93 : _p0(p0), _p1(p1), _geomEdgeInd(iE) {}
94 InSegment() : _p0(0), _p1(0), _geomEdgeInd(0) {}
96 inline bool isConnected( const TVDEdge* edge );
98 inline bool isExternal( const TVDEdge* edge );
100 static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); }
102 static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); }
105 // check if a TVDEdge begins at my end or ends at my start
106 inline bool InSegment::isConnected( const TVDEdge* edge )
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. ));
114 // check if a MA TVDEdge is outside of a domain
115 inline bool InSegment::isExternal( const TVDEdge* edge )
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 );
123 // // -------------------------------------------------------------------------------------
124 // const size_t theExternMA = 111; // to mark external MA edges
126 // bool isExternal( const TVDEdge* edge )
128 // return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA );
131 // // mark external MA edges
132 // void markExternalEdges( const TVDEdge* edge )
134 // if ( isExternal( edge ))
136 // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge );
137 // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() );
138 // if ( edge->is_primary() && edge->vertex1() )
140 // const TVDVertex * v = edge->vertex1();
141 // edge = v->incident_edge();
143 // markExternalEdges( edge );
144 // edge = edge->rot_next();
145 // } while ( edge != v->incident_edge() );
149 // -------------------------------------------------------------------------------------
151 // writes segments into a txt file readable by voronoi_visualizer
152 void inSegmentsToFile( vector< InSegment>& inSegments)
154 if ( inSegments.size() > 1000 )
156 const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
157 SMESH_File file(fileName, false );
158 file.openForWriting();
160 text << "0\n"; // nb points
161 text << inSegments.size() << "\n"; // nb segments
162 for ( size_t i = 0; i < inSegments.size(); ++i )
164 text << inSegments[i]._p0->_a << " "
165 << inSegments[i]._p0->_b << " "
166 << inSegments[i]._p1->_a << " "
167 << inSegments[i]._p1->_b << "\n";
170 file.write( text.c_str(), text.size() );
171 cout << "Write " << fileName << endl;
173 void dumpEdge( const TVDEdge* edge )
175 cout << "*Edge_" << edge;
176 if ( !edge->vertex0() )
177 cout << " ( INF, INF";
179 cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y();
180 if ( !edge->vertex1() )
181 cout << ") -> ( INF, INF";
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" )
193 void dumpCell( const TVDCell* cell )
195 cout << "**Cell_" << cell << " GEOM=" << cell->color() << " ";
196 cout << ( cell->contains_segment() ? " SEG " : " PNT " );
197 if ( cell-> is_degenerate() )
202 const TVDEdge* edge = cell->incident_edge();
206 cout << " - " << ++i << " ";
208 } while (edge != cell->incident_edge());
212 void inSegmentsToFile( vector< InSegment>& inSegments) {}
213 void dumpEdge( const TVDEdge* edge ) {}
214 void dumpCell( const TVDCell* cell ) {}
217 // -------------------------------------------------------------------------------------
223 struct geometry_concept<InPoint> {
224 typedef point_concept type;
227 struct point_traits<InPoint> {
228 typedef int coordinate_type;
230 static inline coordinate_type get(const InPoint& point, orientation_2d orient) {
231 return (orient == HORIZONTAL) ? point._a : point._b;
236 struct geometry_concept<InSegment> {
237 typedef segment_concept type;
241 struct segment_traits<InSegment> {
242 typedef int coordinate_type;
243 typedef InPoint point_type;
245 static inline point_type get(const InSegment& segment, direction_1d dir) {
246 return *(dir.to_int() ? segment._p1 : segment._p0);
249 } // namespace polygon
251 // -------------------------------------------------------------------------------------
255 const int theNoBrachID = 0; // std::numeric_limits<int>::max();
257 // -------------------------------------------------------------------------------------
259 * \brief Intermediate DS to create InPoint's
265 UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {}
266 InPoint getInPoint( double scale[2] )
268 return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u );
271 // -------------------------------------------------------------------------------------
273 * \brief A segment on EDGE, used to create BndPoints
278 const TVDEdge* _edge;
280 int _branchID; // negative ID means reverse direction
282 BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
283 _inSeg(seg), _edge(edge), _uLast(u), _branchID( theNoBrachID ) {}
285 void setIndexToEdge( size_t id )
287 SMESH_MAT2d::Branch::setBndSegment( id, _edge );
290 int branchID() const { return Abs( _branchID ); }
292 size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
294 void setBranch( int branchID, vector< BndSeg >& bndSegs )
296 _branchID = branchID;
298 if ( _edge ) // pass branch to an opposite BndSeg
300 size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
301 if ( oppSegIndex < bndSegs.size() && bndSegs[ oppSegIndex ]._branchID == theNoBrachID )
302 bndSegs[ oppSegIndex ]._branchID = -branchID;
305 bool hasOppositeEdge( const size_t noEdgeID )
307 if ( !_edge ) return false;
308 return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != noEdgeID );
311 // check a next segment in CW order
312 bool isSameBranch( const BndSeg& seg2 )
314 if ( !_edge || !seg2._edge )
317 const TVDCell* cell1 = this->_edge->twin()->cell();
318 const TVDCell* cell2 = seg2. _edge->twin()->cell();
319 if ( cell1 == cell2 )
322 const TVDEdge* edgeMedium1 = this->_edge->twin()->next();
323 const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev();
325 if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() )
327 if ( edgeMedium1->twin() == edgeMedium2 )
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???
333 if ( edgeMedium1->twin() == edgeMedium2->twin()->next() &&
334 edgeMedium1->twin()->cell()->contains_point() )
337 else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
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
350 //================================================================================
352 * \brief Computes length of of TVDEdge
354 //================================================================================
356 double length( const TVDEdge* edge )
358 gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(),
359 edge->vertex0()->y() - edge->vertex1()->y() );
363 //================================================================================
365 * \brief Compute scale to have the same 2d proportions as in 3d
367 //================================================================================
369 void computeProportionScale( const TopoDS_Face& face,
370 const Bnd_B2d& uvBox,
373 scale[0] = scale[1] = 1.;
374 if ( uvBox.IsVoid() ) return;
377 Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
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;
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++)
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 );
399 scale[0] = uLen3d / ( uvMax.X() - uvMin.X() );
400 scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() );
403 //================================================================================
405 * \brief Fill input data for construct_voronoi()
407 //================================================================================
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,
416 const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization
419 // discretize the EDGEs to get 2d points and segments
421 vector< vector< UVU > > uvuVec( edges.size() );
423 for ( size_t iE = 0; iE < edges.size(); ++iE )
425 vector< UVU > & points = uvuVec[ iE ];
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;
432 points.push_back( UVU( c2d->Value( f ), f ));
433 uvBox.Add( points.back()._uv );
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
439 GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect);
440 // if ( discret.NbPoints() > 2 )
445 // discret.Initialize( c2dAdaptor, 100, curDeflect );
446 // cout << "C " << curDeflect << " " << discret.NbPoints() << endl;
447 // curDeflect *= 1.5;
449 // while ( discret.NbPoints() > 5 );
453 // discret.Initialize( c2dAdaptor, angDeflect, 100 );
454 // cout << "A " << angDeflect << " " << discret.NbPoints() << endl;
455 // angDeflect *= 1.5;
457 // while ( discret.NbPoints() > 5 );
461 pPrev = c3d->Value( f );
462 if ( discret.NbPoints() > 2 )
463 for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
465 double u = discret.Parameter(i);
469 int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
470 double dU = ( u - points.back()._u ) / nbDiv;
471 for ( int iD = 1; iD < nbDiv; ++iD )
473 double uD = points.back()._u + dU;
474 points.push_back( UVU( c2d->Value( uD ), uD ));
478 points.push_back( UVU( c2d->Value( u ), u ));
479 uvBox.Add( points.back()._uv );
481 // if ( !c3d.IsNull() )
483 // vector<double> params;
484 // GeomAdaptor_Curve c3dAdaptor( c3d,f,l );
487 // const double deflection = minSegLen * 0.1;
488 // GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true );
489 // if ( !discret.IsDone() )
491 // int nbP = discret.NbPoints();
492 // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
493 // params.push_back( discret.Parameter(i) );
497 // double eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor );
498 // int nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef ));
499 // double segLen = eLen / nbSeg;
500 // GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l );
501 // int nbP = Min( discret.NbPoints(), nbSeg + 1 );
502 // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
503 // params.push_back( discret.Parameter(i) );
505 // for ( size_t i = 0; i < params.size(); ++i )
507 // points.push_back( UVU( c2d->Value( params[i] ), params[i] ));
508 // uvBox.Add( points.back()._uv );
511 if ( points.size() < 2 )
513 points.push_back( UVU( c2d->Value( l ), l ));
514 uvBox.Add( points.back()._uv );
516 if ( edges[ iE ].Orientation() == TopAbs_REVERSED )
517 std::reverse( points.begin(), points.end() );
520 // make connected EDGEs have same UV at shared VERTEX
521 TopoDS_Vertex vShared;
522 for ( size_t iE = 0; iE < edges.size(); ++iE )
524 size_t iE2 = (iE+1) % edges.size();
525 if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared ))
527 if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
529 vector< UVU > & points1 = uvuVec[ iE ];
530 vector< UVU > & points2 = uvuVec[ iE2 ];
531 gp_Pnt2d & uv1 = points1.back() ._uv;
532 gp_Pnt2d & uv2 = points2.front()._uv;
533 uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() );
536 // get scale to have the same 2d proportions as in 3d
537 computeProportionScale( face, uvBox, scale );
539 // make scale to have coordinates precise enough when converted to int
541 gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
542 uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
543 uvMin.ChangeCoord(2) = uvMin.Y() * scale[1];
544 uvMax.ChangeCoord(1) = uvMax.X() * scale[0];
545 uvMax.ChangeCoord(2) = uvMax.Y() * scale[1];
546 double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
547 Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
548 int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
549 const double precision = 1e-5;
550 double preciScale = Min( vMax[iMax] / precision,
551 std::numeric_limits<int>::max() / vMax[iMax] );
552 preciScale /= scale[iMax];
553 double roundedScale = 10; // to ease debug
554 while ( roundedScale * 10 < preciScale )
556 scale[0] *= roundedScale;
557 scale[1] *= roundedScale;
559 // create input points and segments
564 for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
565 nbPnt += uvuVec[ iE ].size();
566 inPoints.resize( nbPnt );
567 inSegments.reserve( nbPnt );
570 if ( face.Orientation() == TopAbs_REVERSED )
572 for ( int iE = uvuVec.size()-1; iE >= 0; --iE )
574 vector< UVU > & points = uvuVec[ iE ];
575 inPoints[ iP++ ] = points.back().getInPoint( scale );
576 for ( size_t i = points.size()-1; i >= 1; --i )
578 inPoints[ iP++ ] = points[i-1].getInPoint( scale );
579 inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
585 for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
587 vector< UVU > & points = uvuVec[ iE ];
588 inPoints[ iP++ ] = points[0].getInPoint( scale );
589 for ( size_t i = 1; i < points.size(); ++i )
591 inPoints[ iP++ ] = points[i].getInPoint( scale );
592 inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
599 //================================================================================
601 * \brief Update a branch joined to another one
603 //================================================================================
605 void updateJoinedBranch( vector< const TVDEdge* > & branchEdges,
607 vector< BndSeg > & bndSegs,
612 for ( size_t i = 0; i < branchEdges.size(); ++i )
614 size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
615 size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
616 bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
617 bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
618 bndSegs[ seg1 ]._branchID *= -newID;
619 bndSegs[ seg2 ]._branchID *= -newID;
620 branchEdges[i] = branchEdges[i]->twin();
622 std::reverse( branchEdges.begin(), branchEdges.end() );
626 for ( size_t i = 0; i < branchEdges.size(); ++i )
628 size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
629 size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
630 bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
631 bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
632 bndSegs[ seg1 ]._branchID *= newID;
633 bndSegs[ seg2 ]._branchID *= newID;
638 //================================================================================
640 * \brief Create MA branches and FACE boundary data
641 * \param [in] vd - voronoi diagram of \a inSegments
642 * \param [in] inPoints - FACE boundary points
643 * \param [in,out] inSegments - FACE boundary segments
644 * \param [out] branch - MA branches to fill
645 * \param [out] branchEnd - ends of MA branches to fill
646 * \param [out] boundary - FACE boundary to fill
648 //================================================================================
650 void makeMA( const TVD& vd,
651 const bool ignoreCorners,
652 vector< InPoint >& inPoints,
653 vector< InSegment > & inSegments,
654 vector< SMESH_MAT2d::Branch >& branch,
655 vector< const SMESH_MAT2d::BranchEnd* >& branchPnt,
656 SMESH_MAT2d::Boundary& boundary )
658 const size_t noEdgeID = inSegments.size() + 1; // ID of non-existent geom EDGE
660 // Associate MA cells with inSegments
661 for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
663 const TVDCell* cell = &(*it);
664 if ( cell->contains_segment() )
666 InSegment& seg = inSegments[ cell->source_index() ];
668 seg.setGeomEdgeToCell( cell, seg._geomEdgeInd );
672 InSegment::setGeomEdgeToCell( cell, noEdgeID );
676 vector< bool > inPntChecked( inPoints.size(), false );
678 // Find MA edges of each inSegment
680 for ( size_t i = 0; i < inSegments.size(); ++i )
682 InSegment& inSeg = inSegments[i];
684 // get edges around the cell lying on MA
685 bool hasSecondary = false;
686 const TVDEdge* edge = inSeg._cell->incident_edge();
688 edge = edge->next(); // Returns the CCW next edge within the cell.
689 if ( edge->is_primary() && !inSeg.isExternal( edge ) )
690 inSeg._edges.push_back( edge ); // edge equidistant from two InSegments
693 } while (edge != inSeg._cell->incident_edge());
695 // there can be several continuous MA edges but maEdges can begin in the middle of
696 // a chain of continuous MA edges. Make the chain continuous.
697 list< const TVDEdge* >& maEdges = inSeg._edges;
698 if ( maEdges.empty() )
701 while ( maEdges.back()->next() == maEdges.front() )
702 maEdges.splice( maEdges.end(), maEdges, maEdges.begin() );
704 // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE
705 list< const TVDEdge* >::iterator e = maEdges.begin();
706 while ( e != maEdges.end() )
708 const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge
709 size_t geoE2 = InSegment::getGeomEdge( cell2 );
710 bool toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e ));
712 e = maEdges.erase( e );
716 if ( maEdges.empty() )
719 // add MA edges corresponding to concave InPoints
720 for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg
722 InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 );
723 size_t pInd = inPnt.index( inPoints );
724 if ( inPntChecked[ pInd ] )
727 inPntChecked[ pInd-1 ] &&
728 inPoints[ pInd-1 ] == inPnt )
730 inPntChecked[ pInd ] = true;
732 const TVDEdge* edge = // a TVDEdge passing through an end of inSeg
733 is2nd ? maEdges.front()->prev() : maEdges.back()->next();
736 if ( edge->is_primary() ) break; // this should not happen
737 const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
738 if ( inSeg.getGeomEdge( edge2->cell() ) != noEdgeID )
739 break; // cell of an InSegment
740 bool hasInfinite = false;
741 list< const TVDEdge* > pointEdges;
745 edge = edge->next(); // Returns the CCW next edge within the cell.
746 if ( edge->is_infinite() )
748 else if ( edge->is_primary() && !inSeg.isExternal( edge ))
749 pointEdges.push_back( edge );
751 while ( edge != edge2 && !hasInfinite );
753 if ( hasInfinite || pointEdges.empty() )
755 inPnt._edges.splice( inPnt._edges.end(), pointEdges );
756 inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd );
758 edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next();
760 } // add MA edges corresponding to concave InPoints
762 } // loop on inSegments to find corresponding MA edges
765 // -------------------------------------------
766 // Create Branches and BndPoints for each EDGE
767 // -------------------------------------------
769 if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ )
771 inPntChecked[0] = false; // do not use the 1st point twice
772 //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), noEdgeID );
773 inPoints[0]._edges.clear();
776 // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
778 vector< BndSeg > bndSegs;
779 bndSegs.reserve( inSegments.size() * 3 );
781 list< const TVDEdge* >::reverse_iterator e;
782 for ( size_t i = 0; i < inSegments.size(); ++i )
784 InSegment& inSeg = inSegments[i];
786 // segments around 1st concave point
787 size_t ip0 = inSeg._p0->index( inPoints );
788 if ( inPntChecked[ ip0 ] )
789 for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
790 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
791 inPntChecked[ ip0 ] = false;
793 // segments of InSegment's
794 const size_t nbMaEdges = inSeg._edges.size();
795 switch ( nbMaEdges ) {
796 case 0: // "around" circle center
797 bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
799 bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
801 gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
802 inSeg._p1->_b - inSeg._p0->_b );
803 const double inSegLen2 = inSegDir.SquareModulus();
804 e = inSeg._edges.rbegin();
805 for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
807 gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
808 (*e)->vertex0()->y() - inSeg._p0->_b );
809 double r = toMA * inSegDir / inSegLen2;
810 double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
811 bndSegs.push_back( BndSeg( &inSeg, *e, u ));
813 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
815 // segments around 2nd concave point
816 size_t ip1 = inSeg._p1->index( inPoints );
817 if ( inPntChecked[ ip1 ] )
818 for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
819 bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
820 inPntChecked[ ip1 ] = false;
823 // make TVDEdge's know it's BndSeg to enable passing branchID to
824 // an opposite BndSeg in BndSeg::setBranch()
825 for ( size_t i = 0; i < bndSegs.size(); ++i )
826 bndSegs[i].setIndexToEdge( i );
829 // Find TVDEdge's of Branches and associate them with bndSegs
831 vector< vector<const TVDEdge*> > branchEdges;
832 branchEdges.reserve( boundary.nbEdges() * 4 );
834 map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
836 int branchID = 1; // we code orientation as branchID sign
837 branchEdges.resize( branchID + 1 );
840 while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
842 bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and the opposite bndSeg
843 branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
845 for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
847 if ( bndSegs[i].branchID() )
849 branchID = bndSegs[i]._branchID; // with sign
851 if ( bndSegs[i]._branchID == -bndSegs[i-1]._branchID &&
854 SMESH_MAT2d::BranchEndType type =
855 ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
856 SMESH_MAT2d::BE_ON_VERTEX :
857 SMESH_MAT2d::BE_END );
858 endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
862 if ( !bndSegs[i-1].isSameBranch( bndSegs[i] ))
864 branchEdges.resize(( branchID = branchEdges.size()) + 1 );
865 if ( bndSegs[i]._edge )
866 endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
867 SMESH_MAT2d::BE_BRANCH_POINT ));
869 bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
870 if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
871 branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
873 // define BranchEndType of the first TVDVertex
874 if ( bndSegs.front()._branchID == -bndSegs.back()._branchID )
876 if ( bndSegs[0]._edge )
878 SMESH_MAT2d::BranchEndType type =
879 ( bndSegs[0]._inSeg->isConnected( bndSegs[0]._edge ) ?
880 SMESH_MAT2d::BE_ON_VERTEX :
881 SMESH_MAT2d::BE_END );
882 endType.insert( make_pair( bndSegs[0]._edge->vertex1(), type ));
884 else if ( bndSegs.back()._edge )
886 SMESH_MAT2d::BranchEndType type =
887 ( bndSegs.back()._inSeg->isConnected( bndSegs.back()._edge ) ?
888 SMESH_MAT2d::BE_ON_VERTEX :
889 SMESH_MAT2d::BE_END );
890 endType.insert( make_pair( bndSegs.back()._edge->vertex0(), type ));
893 // join the 1st and the last branch edges if it is the same branch
894 if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
895 bndSegs.back().isSameBranch( bndSegs.front() ))
897 vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
898 vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID() ];
899 br1.insert( br1.begin(), br2.begin(), br2.end() );
903 // remove branches ending at BE_ON_VERTEX
905 vector<bool> isBranchRemoved( branchEdges.size(), false );
907 if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
909 // find branches to remove
910 map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et;
911 for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
913 if ( branchEdges[iB].empty() )
915 const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
916 const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
917 v2et = endType.find( v0 );
918 if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
919 isBranchRemoved[ iB ] = true;
920 v2et = endType.find( v1 );
921 if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
922 isBranchRemoved[ iB ] = true;
924 // try to join not removed branches into one
925 for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
927 if ( branchEdges[iB].empty() || isBranchRemoved[iB] )
929 const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
930 const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
931 v2et = endType.find( v0 );
932 if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
934 v2et = endType.find( v1 );
935 if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
940 size_t iBrToJoin = 0;
941 for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
943 if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
945 const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
946 const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
947 if ( v0 == v02 || v0 == v12 || v1 == v02 || v1 == v12 )
952 break; // more than 2 not removed branches meat at a TVDVertex
959 vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
960 const TVDVertex* v02 = branch[0]->vertex1();
961 const TVDVertex* v12 = branch.back()->vertex0();
962 updateJoinedBranch( branch, iB, bndSegs, /*reverse=*/(v0 == v02 || v1 == v12 ));
963 if ( v0 == v02 || v0 == v12 )
964 branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
966 branchEdges[iB].insert( branchEdges[iB].end(), branch.begin(), branch.end() );
969 } // loop on branchEdges
970 } // if ( ignoreCorners )
972 // associate branchIDs and the input branch vector (arg)
973 vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 );
975 for ( size_t i = 0; i < branchEdges.size(); ++i )
977 nbBranches += ( !branchEdges[i].empty() );
979 branch.resize( nbBranches );
981 for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed
983 if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ])
984 branchByID[ brID ] = & branch[ iBr++ ];
986 for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed
988 if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ])
989 branchByID[ brID ] = & branch[ iBr++ ];
992 // Fill in BndPoints of each EDGE of the boundary
995 int edgeInd = -1, dInd = 0;
996 while ( iSeg < bndSegs.size() )
998 const size_t geomID = bndSegs[ iSeg ].geomEdge();
999 SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( geomID );
1002 for ( size_t i = iSeg; i < bndSegs.size() && geomID == bndSegs[ i ].geomEdge(); ++i )
1004 size_t iSegEnd = iSeg + nbSegs;
1006 // make TVDEdge know an index of bndSegs within BndPoints
1007 for ( size_t i = iSeg; i < iSegEnd; ++i )
1008 if ( bndSegs[i]._edge )
1009 SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge );
1011 // parameters on EDGE
1013 bndPoints._params.reserve( nbSegs + 1 );
1014 bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param );
1016 for ( size_t i = iSeg; i < iSegEnd; ++i )
1017 bndPoints._params.push_back( bndSegs[ i ]._uLast );
1021 bndPoints._maEdges.reserve( nbSegs );
1023 for ( size_t i = iSeg; i < iSegEnd; ++i )
1025 const size_t brID = bndSegs[ i ].branchID();
1026 const SMESH_MAT2d::Branch* br = branchByID[ brID ];
1028 if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
1032 if (( edgeInd < 0 ||
1033 edgeInd >= (int) branchEdges[ brID ].size() ) ||
1034 ( branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge &&
1035 branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge ))
1037 if ( bndSegs[ i ]._branchID < 0 )
1040 for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd )
1041 if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge )
1044 else // bndSegs[ i ]._branchID > 0
1047 for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
1048 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
1055 // no MA edge, bndSeg corresponds to an end point of a branch
1056 if ( bndPoints._maEdges.empty() )
1059 edgeInd = branchEdges[ brID ].size();
1060 dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
1063 bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
1065 } // loop on bndSegs of an EDGE
1069 } // loop on all bndSegs
1071 // Initialize branches
1073 // find a not removed branch
1074 size_t iBrNorRemoved = 0;
1075 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1076 if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] )
1078 iBrNorRemoved = brID;
1081 // fill the branches with MA edges
1082 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1083 if ( !branchEdges[brID].empty() )
1085 branchByID[ brID ]->init( branchEdges[brID], & boundary, endType );
1087 // mark removed branches
1088 for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
1089 if ( isBranchRemoved[brID] && iBrNorRemoved > 0 )
1091 SMESH_MAT2d::Branch* branch = branchByID[ brID ];
1092 SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ];
1093 bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT );
1094 const TVDVertex* branchVextex =
1095 is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex;
1096 SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex );
1097 branch->setRemoved( bp );
1099 // set branches to branch ends
1100 for ( size_t i = 0; i < branch.size(); ++i )
1101 if ( !branch[i].isRemoved() )
1102 branch[i].setBranchesToEnds( branch );
1104 // fill branchPnt arg
1105 map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
1106 for ( size_t i = 0; i < branch.size(); ++i )
1108 if ( branch[i].getEnd(0)->_branches.size() > 2 )
1109 v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
1110 if ( branch[i].getEnd(1)->_branches.size() > 2 )
1111 v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
1113 branchPnt.resize( v2end.size() );
1114 map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
1115 for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
1116 branchPnt[ i ] = v2e->second;
1122 //================================================================================
1124 * \brief MedialAxis constructor
1125 * \param [in] face - a face to create MA for
1126 * \param [in] edges - edges of the face (possibly not all) on the order they
1127 * encounter in the face boundary.
1128 * \param [in] minSegLen - minimal length of a mesh segment used to discretize
1129 * the edges. It is used to define precision of MA approximation
1131 //================================================================================
1133 SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face,
1134 const std::vector< TopoDS_Edge >& edges,
1135 const double minSegLen,
1136 const bool ignoreCorners):
1137 _face( face ), _boundary( edges.size() )
1139 // input to construct_voronoi()
1140 vector< InPoint > inPoints;
1141 vector< InSegment> inSegments;
1142 if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
1145 inSegmentsToFile( inSegments );
1147 // build voronoi diagram
1148 construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
1151 makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary );
1153 // count valid branches
1154 _nbBranches = _branch.size();
1155 for ( size_t i = 0; i < _branch.size(); ++i )
1156 if ( _branch[i].isRemoved() )
1160 //================================================================================
1162 * \brief Returns the i-th branch
1164 //================================================================================
1166 const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const
1168 return i < _nbBranches ? &_branch[i] : 0;
1171 //================================================================================
1173 * \brief Return UVs of ends of MA edges of a branch
1175 //================================================================================
1177 void SMESH_MAT2d::MedialAxis::getPoints( const Branch* branch,
1178 std::vector< gp_XY >& points) const
1180 branch->getPoints( points, _scale );
1183 //================================================================================
1185 * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
1186 * \param [in] iEdge - index of geom EDGE within a vector passed at MA construction
1187 * \param [in] u - parameter of the point on EDGE curve
1188 * \param [out] p - the found BranchPoint
1189 * \return bool - is OK
1191 //================================================================================
1193 bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
1195 BranchPoint& p ) const
1197 if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1200 const BndPoints& points = _pointsPerEdge[ iEdge ];
1201 const bool edgeReverse = ( points._params[0] > points._params.back() );
1203 if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
1204 u = edgeReverse ? points._params.back() : points._params[0];
1205 else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
1206 u = edgeReverse ? points._params[0] : points._params.back();
1208 double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
1209 int i = int( r * double( points._maEdges.size()-1 ));
1212 while ( points._params[i ] < u ) --i;
1213 while ( points._params[i+1] > u ) ++i;
1217 while ( points._params[i ] > u ) --i;
1218 while ( points._params[i+1] < u ) ++i;
1221 double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
1223 if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
1225 if ( i < points._maEdges.size() / 2 ) // near 1st point
1227 while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second )
1229 edgeParam = edgeReverse;
1231 else // near last point
1233 while ( i > 0 && !points._maEdges[ i ].second )
1235 edgeParam = !edgeReverse;
1238 const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
1239 bool maReverse = ( maE.second < 0 );
1241 p._branch = maE.first;
1242 p._iEdge = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
1243 p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam;
1248 //================================================================================
1250 * \brief Check if a given boundary segment is a null-length segment on a concave
1252 * \param [in] iEdge - index of a geom EDGE
1253 * \param [in] iSeg - index of a boundary segment
1254 * \return bool - true if the segment is on concave corner
1256 //================================================================================
1258 bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
1260 if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
1263 const BndPoints& points = _pointsPerEdge[ iEdge ];
1264 if ( points._params.size() <= iSeg+1 )
1267 return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20;
1270 //================================================================================
1272 * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end
1274 //================================================================================
1276 bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
1278 if ( bp._edgeIndex >= _pointsPerEdge.size() )
1281 const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
1282 if ( bp._param - points._params[0] < points._params.back() - bp._param )
1283 bp._param = points._params[0];
1285 bp._param = points._params.back();
1290 //================================================================================
1292 * \brief Creates a 3d curve corresponding to a Branch
1293 * \param [in] branch - the Branch
1294 * \return Adaptor3d_Curve* - the new curve the caller is to delete
1296 //================================================================================
1298 Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
1300 Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
1301 if ( surface.IsNull() )
1305 branch.getPoints( uv, _scale );
1306 if ( uv.size() < 2 )
1309 vector< TopoDS_Vertex > vertex( uv.size() );
1310 for ( size_t i = 0; i < uv.size(); ++i )
1311 vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
1314 BRep_Builder aBuilder;
1315 aBuilder.MakeWire(aWire);
1316 for ( size_t i = 1; i < vertex.size(); ++i )
1318 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
1319 aBuilder.Add( aWire, edge );
1322 // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
1323 // aWire.Closed(true); // issue 0021141
1325 return new BRepAdaptor_CompCurve( aWire );
1328 //================================================================================
1330 * \brief Copy points of an EDGE
1332 //================================================================================
1334 void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>& maEdges,
1335 const Boundary* boundary,
1336 map< const TVDVertex*, BranchEndType > endType )
1338 if ( maEdges.empty() ) return;
1340 _boundary = boundary;
1341 _maEdges.swap( maEdges );
1344 _params.reserve( _maEdges.size() + 1 );
1345 _params.push_back( 0. );
1346 for ( size_t i = 0; i < _maEdges.size(); ++i )
1347 _params.push_back( _params.back() + length( _maEdges[i] ));
1349 for ( size_t i = 1; i < _params.size(); ++i )
1350 _params[i] /= _params.back();
1353 _endPoint1._vertex = _maEdges.front()->vertex1();
1354 _endPoint2._vertex = _maEdges.back ()->vertex0();
1356 if ( endType.count( _endPoint1._vertex ))
1357 _endPoint1._type = endType[ _endPoint1._vertex ];
1358 if ( endType.count( _endPoint2._vertex ))
1359 _endPoint2._type = endType[ _endPoint2._vertex ];
1362 //================================================================================
1364 * \brief fills BranchEnd::_branches of its ends
1366 //================================================================================
1368 void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
1370 for ( size_t i = 0; i < branches.size(); ++i )
1372 if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
1373 this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
1374 this->_endPoint1._branches.push_back( &branches[i] );
1376 if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
1377 this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
1378 this->_endPoint2._branches.push_back( &branches[i] );
1382 //================================================================================
1384 * \brief returns a BranchPoint corresponding to a TVDVertex
1386 //================================================================================
1388 SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const
1394 if ( vertex == _maEdges[0]->vertex1() )
1400 for ( ; p._iEdge < _maEdges.size(); ++p._iEdge )
1401 if ( vertex == _maEdges[ p._iEdge ]->vertex0() )
1403 p._edgeParam = _params[ p._iEdge ];
1410 //================================================================================
1412 * \brief Sets a proxy point for a removed branch
1413 * \param [in] proxyPoint - a point of another branch to which all points of this
1416 //================================================================================
1418 void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint )
1420 _proxyPoint = proxyPoint;
1423 //================================================================================
1425 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1426 * \param [in] param - [0;1] normalized param on the Branch
1427 * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1428 * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1429 * \return bool - true if the BoundaryPoint's found
1431 //================================================================================
1433 bool SMESH_MAT2d::Branch::getBoundaryPoints(double param,
1435 BoundaryPoint& bp2 ) const
1437 if ( param < _params[0] || param > _params.back() )
1440 // look for an index of a MA edge by param
1441 double ip = param * _params.size();
1442 size_t i = size_t( Min( int( _maEdges.size()-1), int( ip )));
1444 while ( param < _params[i ] ) --i;
1445 while ( param > _params[i+1] ) ++i;
1447 double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
1449 return getBoundaryPoints( i, r, bp1, bp2 );
1452 //================================================================================
1454 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1455 * \param [in] iMAEdge - index of a MA edge within this Branch
1456 * \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
1457 * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
1458 * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
1459 * \return bool - true if the BoundaryPoint's found
1461 //================================================================================
1463 bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t iMAEdge,
1466 BoundaryPoint& bp2 ) const
1469 return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 );
1471 if ( iMAEdge > _maEdges.size() )
1473 if ( iMAEdge == _maEdges.size() )
1474 iMAEdge = _maEdges.size() - 1;
1476 size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
1477 size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
1478 size_t iSeg1 = getBndSegment( _maEdges[ iMAEdge ] );
1479 size_t iSeg2 = getBndSegment( _maEdges[ iMAEdge ]->twin() );
1481 return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
1482 _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
1485 //================================================================================
1487 * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
1489 //================================================================================
1491 bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
1493 BoundaryPoint& bp2 ) const
1495 return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
1498 //================================================================================
1500 * \brief Return a parameter of a BranchPoint normalized within this Branch
1502 //================================================================================
1504 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
1506 if ( this != p._branch && p._branch )
1507 return p._branch->getParameter( p, u );
1510 return _proxyPoint._branch->getParameter( _proxyPoint, u );
1512 if ( p._iEdge > _params.size()-1 )
1514 if ( p._iEdge == _params.size()-1 )
1517 u = ( _params[ p._iEdge ] * ( 1 - p._edgeParam ) +
1518 _params[ p._iEdge+1 ] * p._edgeParam );
1523 //================================================================================
1525 * \brief Check type of both ends
1527 //================================================================================
1529 bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
1531 return ( _endPoint1._type == type || _endPoint2._type == type );
1534 //================================================================================
1536 * \brief Returns MA points
1537 * \param [out] points - the 2d points
1538 * \param [in] scale - the scale that was used to scale the 2d space of MA
1540 //================================================================================
1542 void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
1543 const double scale[2]) const
1545 points.resize( _maEdges.size() + 1 );
1547 points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
1548 _maEdges[0]->vertex1()->y() / scale[1] );
1550 for ( size_t i = 0; i < _maEdges.size(); ++i )
1551 points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
1552 _maEdges[i]->vertex0()->y() / scale[1] );
1555 //================================================================================
1557 * \brief Return indices of EDGEs equidistant from this branch
1559 //================================================================================
1561 void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
1562 std::vector< std::size_t >& edgeIDs2 ) const
1564 edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1565 edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1567 for ( size_t i = 1; i < _maEdges.size(); ++i )
1569 size_t ie1 = getGeomEdge( _maEdges[i] );
1570 size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1572 if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
1573 if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
1577 //================================================================================
1579 * \brief Looks for a BranchPoint position around a concave VERTEX
1581 //================================================================================
1583 bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1,
1584 std::vector< std::size_t >& edgeIDs2,
1585 std::vector< BranchPoint >& divPoints,
1586 const vector<const TVDEdge*>& maEdges,
1587 const vector<const TVDEdge*>& maEdgesTwin,
1590 // if there is a concave vertex between EDGEs
1591 // then position of a dividing BranchPoint is undefined, it is somewhere
1592 // on an arc-shaped part of the Branch around the concave vertex.
1593 // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
1594 // of the arc if there is no opposite VERTEX.
1595 // All null-length segments around a VERTEX belong to one of EDGEs.
1597 BranchPoint divisionPnt;
1598 divisionPnt._branch = this;
1600 size_t ie1 = getGeomEdge( maEdges [i] );
1601 size_t ie2 = getGeomEdge( maEdgesTwin[i] );
1603 size_t iSeg1 = getBndSegment( maEdges[ i-1 ] );
1604 size_t iSeg2 = getBndSegment( maEdges[ i ] );
1605 bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
1606 bool isConcaNext = _boundary->isConcaveSegment( ie1, iSeg2 );
1607 if ( !isConcaNext && !isConcaPrev )
1610 bool isConcaveV = false;
1612 int iPrev = i-1, iNext = i;
1613 if ( isConcaNext ) // all null-length segments follow
1615 // look for a VERTEX of the opposite EDGE
1616 ++iNext; // end of null-length segments
1617 while ( iNext < maEdges.size() )
1619 iSeg2 = getBndSegment( maEdges[ iNext ] );
1620 if ( _boundary->isConcaveSegment( ie1, iSeg2 ))
1625 bool vertexFound = false;
1626 for ( size_t iE = i+1; iE < iNext; ++iE )
1628 ie2 = getGeomEdge( maEdgesTwin[iE] );
1629 if ( ie2 != edgeIDs2.back() )
1631 // opposite VERTEX found
1632 divisionPnt._iEdge = iE;
1633 divisionPnt._edgeParam = 0;
1634 divPoints.push_back( divisionPnt );
1635 edgeIDs1.push_back( ie1 );
1636 edgeIDs2.push_back( ie2 );
1642 iPrev = i = --iNext; // not to add a BP in the moddle
1646 else if ( isConcaPrev )
1648 // all null-length segments passed, find their beginning
1649 while ( iPrev-1 >= 0 )
1651 iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
1652 if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
1659 if ( iPrev < i-1 || iNext > i )
1661 // no VERTEX on the opposite EDGE, put the Branch Point in the middle
1662 double par1 = _params[ iPrev+1 ], par2 = _params[ iNext ];
1663 double midPar = 0.5 * ( par1 + par2 );
1664 divisionPnt._iEdge = iPrev;
1665 while ( _params[ divisionPnt._iEdge + 1 ] < midPar )
1666 ++divisionPnt._iEdge;
1667 divisionPnt._edgeParam =
1668 ( _params[ divisionPnt._iEdge + 1 ] - midPar ) /
1669 ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] );
1670 divPoints.push_back( divisionPnt );
1677 //================================================================================
1679 * \brief Return indices of opposite parts of EDGEs equidistant from this branch
1680 * \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
1681 * \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
1682 * \param [out] divPoints - BranchPoint's located between two successive unique
1683 * pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
1684 * of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
1685 * than number of \a edgeIDs
1687 //================================================================================
1689 void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
1690 std::vector< std::size_t >& edgeIDs2,
1691 std::vector< BranchPoint >& divPoints) const
1697 edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
1698 edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
1700 std::vector<const TVDEdge*> twins( _maEdges.size() );
1701 for ( size_t i = 0; i < _maEdges.size(); ++i )
1702 twins[i] = _maEdges[i]->twin();
1704 // size_t lastConcaE1 = _boundary.nbEdges();
1705 // size_t lastConcaE2 = _boundary.nbEdges();
1707 BranchPoint divisionPnt;
1708 divisionPnt._branch = this;
1710 for ( size_t i = 0; i < _maEdges.size(); ++i )
1712 size_t ie1 = getGeomEdge( _maEdges[i] );
1713 size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
1715 if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
1717 bool isConcaveV = false;
1718 if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 )
1720 isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i );
1722 if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 )
1724 isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i );
1729 ie1 = getGeomEdge( _maEdges[i] );
1730 ie2 = getGeomEdge( _maEdges[i]->twin() );
1732 if (( !isConcaveV ) ||
1733 ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ))
1735 edgeIDs1.push_back( ie1 );
1736 edgeIDs2.push_back( ie2 );
1738 if ( divPoints.size() < edgeIDs1.size() - 1 )
1740 divisionPnt._iEdge = i;
1741 divisionPnt._edgeParam = 0;
1742 divPoints.push_back( divisionPnt );
1745 } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
1746 } // loop on _maEdges
1749 //================================================================================
1751 * \brief Store data of boundary segments in TVDEdge
1753 //================================================================================
1755 void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
1757 if ( maEdge ) maEdge->cell()->color( geomIndex );
1759 std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
1761 return maEdge ? maEdge->cell()->color() : std::string::npos;
1763 void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
1765 if ( maEdge ) maEdge->color( segIndex );
1767 std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
1769 return maEdge ? maEdge->color() : std::string::npos;
1772 //================================================================================
1774 * \brief Returns a boundary point on a given EDGE
1775 * \param [in] iEdge - index of the EDGE within MedialAxis
1776 * \param [in] iSeg - index of a boundary segment within this Branch
1777 * \param [in] u - [0;1] normalized param within \a iSeg-th segment
1778 * \param [out] bp - the found BoundaryPoint
1779 * \return bool - true if the BoundaryPoint is found
1781 //================================================================================
1783 bool SMESH_MAT2d::Boundary::getPoint( std::size_t iEdge,
1786 BoundaryPoint& bp ) const
1788 if ( iEdge >= _pointsPerEdge.size() )
1790 if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
1793 // This method is called by Branch that can have an opposite orientation,
1794 // hence u is inverted depending on orientation coded as a sign of _maEdge index
1795 bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
1799 double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
1800 double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
1802 bp._param = p0 * ( 1. - u ) + p1 * u;
1803 bp._edgeIndex = iEdge;