+namespace
+{
+ /*!
+ * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
+ */
+ struct PolyVertex
+ {
+ SMESH_TNodeXYZ _nxyz;
+ gp_XY _xy;
+ PolyVertex* _prev;
+ PolyVertex* _next;
+
+ void SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v )
+ {
+ _nxyz.Set( n );
+ _next = &v;
+ v._prev = this;
+ }
+ PolyVertex* Delete()
+ {
+ _prev->_next = _next;
+ _next->_prev = _prev;
+ return _next;
+ }
+ void GetTriaNodes( const SMDS_MeshNode** nodes) const
+ {
+ nodes[0] = _prev->_nxyz._node;
+ nodes[1] = this->_nxyz._node;
+ nodes[2] = _next->_nxyz._node;
+ }
+
+ inline static double Area( const PolyVertex* v0, const PolyVertex* v1, const PolyVertex* v2 )
+ {
+ gp_XY vPrev = v0->_xy - v1->_xy;
+ gp_XY vNext = v2->_xy - v1->_xy;
+ return vNext ^ vPrev;
+ }
+ double TriaArea() const { return Area( _prev, this, _next ); }
+
+ bool IsInsideTria( const PolyVertex* v )
+ {
+ gp_XY p = _prev->_xy - v->_xy;
+ gp_XY t = this->_xy - v->_xy;
+ gp_XY n = _next->_xy - v->_xy;
+ const double tol = -1e-12;
+ return (( p ^ t ) >= tol &&
+ ( t ^ n ) >= tol &&
+ ( n ^ p ) >= tol );
+ // return ( Area( _prev, this, v ) > 0 &&
+ // Area( this, _next, v ) > 0 &&
+ // Area( _next, _prev, v ) > 0 );
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Triangulate a polygon. Assure correct orientation for concave polygons
+ */
+ //================================================================================
+
+ bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes )
+ {
+ // connect nodes into a ring
+ std::vector< PolyVertex > pv( nbNodes );
+ for ( size_t i = 1; i < nbNodes; ++i )
+ pv[i-1].SetNodeAndNext( nodes[i-1], pv[i] );
+ pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], pv[0] );
+
+ // get a polygon normal
+ gp_XYZ normal(0,0,0), p0,v01,v02;
+ p0 = pv[0]._nxyz;
+ v01 = pv[1]._nxyz - p0;
+ for ( size_t i = 2; i < nbNodes; ++i )
+ {
+ v02 = pv[i]._nxyz - p0;
+ normal += v01 ^ v02;
+ v01 = v02;
+ }
+ // project nodes to the found plane
+ gp_Ax2 axes;
+ try {
+ axes = gp_Ax2( p0, normal, v01 );
+ }
+ catch ( Standard_Failure ) {
+ return false;
+ }
+ for ( size_t i = 0; i < nbNodes; ++i )
+ {
+ gp_XYZ p = pv[i]._nxyz - p0;
+ pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
+ pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
+ }
+
+ // in a loop, find triangles with positive area and having no vertices inside
+ int iN = 0, nbTria = nbNodes - 2;
+ nodes.reserve( nbTria * 3 );
+ const double minArea = 1e-6;
+ PolyVertex* v = &pv[0], *vi;
+ int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
+ while ( nbBadTria < nbVertices )
+ {
+ if (( isGoodTria = v->TriaArea() > minArea ))
+ {
+ for ( vi = v->_next->_next;
+ vi != v->_prev;
+ vi = vi->_next )
+ {
+ if ( v->IsInsideTria( vi ))
+ break;
+ }
+ isGoodTria = ( vi == v->_prev );
+ }
+ if ( isGoodTria )
+ {
+ v->GetTriaNodes( &nodes[ iN ] );
+ iN += 3;
+ v = v->Delete();
+ if ( --nbVertices == 3 )
+ {
+ // last triangle remains
+ v->GetTriaNodes( &nodes[ iN ] );
+ return true;
+ }
+ nbBadTria = 0;
+ }
+ else
+ {
+ v = v->_next;
+ ++nbBadTria;
+ }
+ }
+
+ // the polygon is invalid; add triangles with positive area
+ nbBadTria = 0;
+ while ( nbBadTria < nbVertices )
+ {
+ isGoodTria = v->TriaArea() > minArea;
+ if ( isGoodTria )
+ {
+ v->GetTriaNodes( &nodes[ iN ] );
+ iN += 3;
+ v = v->Delete();
+ if ( --nbVertices == 3 )
+ {
+ // last triangle remains
+ v->GetTriaNodes( &nodes[ iN ] );
+ return true;
+ }
+ nbBadTria = 0;
+ }
+ else
+ {
+ v = v->_next;
+ ++nbBadTria;
+ }
+ }
+
+ // add all the rest triangles
+ while ( nbVertices >= 3 )
+ {
+ v->GetTriaNodes( &nodes[ iN ] );
+ iN += 3;
+ v = v->Delete();
+ --nbVertices;
+ }
+
+ return true;
+
+ } // triangulate()
+} // namespace
+