1 // Copyright (C) 2007-2020 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_Triangulate.cxx
23 // Created : Thu Jan 18 18:00:13 2018
24 // Author : Edward AGAPOV (eap)
26 // Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx
28 #include "SMESH_MeshAlgos.hxx"
30 #include <Standard_ErrorHandler.hxx>
31 #include <Standard_Failure.hxx>
34 #include <boost/container/flat_set.hpp>
36 using namespace SMESH_MeshAlgos;
40 struct Node // node of a triangle
42 size_t _triaIndex; // triangle index == index of the 1st triangle node in triangulation array
43 size_t _nodeIndex; // node index within triangle [0-2]
45 //! return node index within the node array
46 size_t Index() const { return _triaIndex + _nodeIndex; }
48 //! return local 3-d index [0-2]
49 static size_t ThirdIndex( size_t i1, size_t i2 )
51 size_t i3 = ( i2 + 1 ) % 3;
56 //! return 3-d node index within the node array
57 static size_t ThirdIndex( const Node& n1, const Node& n2 )
59 return n1._triaIndex + ThirdIndex( n1._nodeIndex, n2._nodeIndex );
61 bool operator<(const Node& other) const { return _triaIndex < other._triaIndex; }
63 typedef boost::container::flat_set< Node > TriaNodeSet;
67 * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
69 struct Triangulate::PolyVertex
77 void SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index );
78 void GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const;
79 double TriaArea() const;
80 bool IsInsideTria( const PolyVertex* v );
83 // compare PolyVertex'es by node
84 bool operator()(const PolyVertex* a, const PolyVertex* b) const
86 return ( a->_nxyz.Node() < b->_nxyz.Node() );
88 // set of PolyVertex sorted by mesh node
89 typedef boost::container::flat_set< PolyVertex*, PolyVertex > PVSet;
92 struct Triangulate::Data
94 std::vector< PolyVertex > _pv;
95 std::vector< size_t > _nodeIndex;
96 PolyVertex::PVSet _uniqueNodePV;
99 struct Triangulate::Optimizer
101 std::vector< TriaNodeSet > _nodeUsage; // inclusions of a node in triangles
103 //================================================================================
105 * \brief Optimize triangles by edge swapping
106 * \param [inout] nodes - polygon triangulation, i.e. connectivity of all triangles to optimize
107 * \param [in] points - coordinates of nodes of the input polygon
108 * \param [in] nodeIndices - indices of triangulation nodes within the input polygon
110 //================================================================================
112 void optimize( std::vector< const SMDS_MeshNode*>& nodes,
113 std::vector< PolyVertex > & points,
114 std::vector< size_t > & nodeIndices)
116 // for each node of the polygon, remember triangles using it
117 _nodeUsage.resize( points.size() );
118 for ( size_t i = 0; i < points.size(); ++i ) // clear old data
120 _nodeUsage[ i ].clear();
122 for ( size_t i = 0, iTria = 0; i < nodeIndices.size(); ++iTria )
124 _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 0 });
125 _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 1 });
126 _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 2 });
130 for ( size_t iTria = 0; iTria < nodeIndices.size(); iTria += 3 )
132 double badness1 = computeBadness( nodeIndices[ iTria + 0 ],
133 nodeIndices[ iTria + 1 ],
134 nodeIndices[ iTria + 2 ],
136 for ( size_t i = 0; i < 3; ++i ) // loop on triangle edges to find a neighbor triangle
138 size_t i1 = iTria + i; // node index in nodeIndices
139 size_t i2 = iTria + ( i + 1 ) % 3;
140 size_t ind1 = nodeIndices[ i1 ]; // node index in points
141 size_t ind2 = nodeIndices[ i2 ];
142 TriaNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
143 TriaNodeSet & usage2 = _nodeUsage[ ind2 ];
144 if ( usage1.size() < 2 ||
148 // look for another triangle using two nodes
149 TriaNodeSet::iterator usIt1 = usage1.begin();
150 for ( ; usIt1 != usage1.end(); ++usIt1 )
152 if ( usIt1->_triaIndex == iTria )
153 continue; // current triangle
154 TriaNodeSet::iterator usIt2 = usage2.find( *usIt1 );
155 if ( usIt2 == usage2.end() )
156 continue; // no common _triaIndex in two usages
158 size_t i3 = iTria + ( i + 2 ) % 3;
159 size_t i4 = Node::ThirdIndex( *usIt1, *usIt2 ); // 4th node of quadrangle
160 size_t ind3 = nodeIndices[ i3 ];
161 size_t ind4 = nodeIndices[ i4 ];
163 double badness2 = computeBadness( ind2, ind1, ind4, points );
164 double badness3 = computeBadness( ind1, ind4, ind3, points, /*checkArea=*/true );
165 double badness4 = computeBadness( ind2, ind3, ind4, points, /*checkArea=*/true );
167 if ( Max( badness1, badness2 ) < Max( badness3, badness4 ))
170 // swap edge by modifying nodeIndices
172 nodeIndices[ i2 ] = ind4;
173 _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria });
174 _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
177 nodeIndices[ i1 ] = ind3;
178 _nodeUsage[ ind3 ].insert( *usIt1 );
179 _nodeUsage[ ind1 ].erase ( *usIt1 );
181 --i; // to re-check a current edge
188 // update nodes by updated nodeIndices
189 for ( size_t i = 0; i < nodeIndices.size(); ++i )
190 nodes[ i ] = points[ nodeIndices[ i ]]._nxyz.Node();
195 //================================================================================
197 * \brief Return 1./area. Initially: max cos^2 of triangle angles
199 //================================================================================
201 double computeBadness( size_t i1, size_t i2, size_t i3,
202 std::vector< PolyVertex > & points,
203 bool checkArea = false )
207 points[ i2 ]._prev = & points[ i1 ];
208 points[ i2 ]._next = & points[ i3 ];
209 double a = points[ i2 ].TriaArea();
211 // return std::numeric_limits<double>::max();
217 const gp_XY & p1 = points[ i1 ]._xy;
218 const gp_XY & p2 = points[ i2 ]._xy;
219 const gp_XY & p3 = points[ i3 ]._xy;
220 gp_XY vec[3] = { p2 - p1,
223 double len[3] = { vec[0].SquareModulus(),
224 vec[1].SquareModulus(),
225 vec[2].SquareModulus() };
226 if ( len[0] < gp::Resolution() ||
227 len[1] < gp::Resolution() ||
228 len[2] < gp::Resolution() )
232 for ( int i = 0; i < 3; ++i )
234 int i2 = ( i+1 ) % 3;
235 double dot = -vec[ i ] * vec[ i2 ];
237 maxCos2 = Max( maxCos2, dot * dot / len[ i ] / len[ i2 ] );
243 //================================================================================
245 * \brief Initialization
247 //================================================================================
249 void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
258 //================================================================================
260 * \brief Remove self from a polygon
262 //================================================================================
264 Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
266 _prev->_next = _next;
267 _next->_prev = _prev;
271 //================================================================================
273 * \brief Return nodes of a triangle
275 //================================================================================
277 void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes,
278 size_t* nodeIndices) const
280 nodes[0] = _prev->_nxyz._node;
281 nodes[1] = this->_nxyz._node;
282 nodes[2] = _next->_nxyz._node;
283 nodeIndices[0] = _prev->_index;
284 nodeIndices[1] = this->_index;
285 nodeIndices[2] = _next->_index;
288 //================================================================================
290 * \brief Compute triangle area
292 //================================================================================
294 inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
296 gp_XY vPrev = xy0 - xy1;
297 gp_XY vNext = xy2 - xy1;
298 return vNext ^ vPrev;
301 //================================================================================
303 * \brief Compute triangle area
305 //================================================================================
307 double Triangulate::PolyVertex::TriaArea() const
309 return Area( _prev->_xy, this->_xy, _next->_xy );
312 //================================================================================
314 * \brief Check if a vertex is inside a triangle
316 //================================================================================
318 bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
320 if ( this ->_nxyz == v->_nxyz ||
321 _prev->_nxyz == v->_nxyz ||
322 _next->_nxyz == v->_nxyz )
325 gp_XY p = _prev->_xy - v->_xy;
326 gp_XY t = this->_xy - v->_xy;
327 gp_XY n = _next->_xy - v->_xy;
328 const double tol = -1e-7;
329 return (( p ^ t ) >= tol &&
332 // return ( Area( _prev, this, v ) > 0 &&
333 // Area( this, _next, v ) > 0 &&
334 // Area( _next, _prev, v ) > 0 );
337 //================================================================================
339 * \brief Triangulate a polygon. Assure correct orientation for concave polygons
341 //================================================================================
343 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
344 const size_t nbNodes)
346 std::vector< PolyVertex >& _pv = _data->_pv;
347 std::vector< size_t >& _nodeIndex = _data->_nodeIndex;
348 PolyVertex::PVSet& _uniqueNodePV = _data->_uniqueNodePV;
350 // connect nodes into a ring
351 _pv.resize( nbNodes );
352 for ( size_t i = 1; i < nbNodes; ++i )
353 _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], /*index=*/i-1 );
354 _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0], nbNodes-1 );
356 // assure correctness of PolyVertex::_index as a node can encounter more than once
357 // within a polygon boundary
358 if ( _optimizer && nbNodes > 4 )
360 _uniqueNodePV.clear();
361 for ( size_t i = 0; i < nbNodes; ++i )
363 PolyVertex::PVSet::iterator pv = _uniqueNodePV.insert( &_pv[i] ).first;
364 _pv[i]._index = (*pv)->_index;
368 // get a polygon normal
369 gp_XYZ normal(0,0,0), p0,v01,v02;
371 v01 = _pv[1]._nxyz - p0;
372 for ( size_t i = 2; i < nbNodes; ++i )
374 v02 = _pv[i]._nxyz - p0;
378 // project nodes to the found plane
381 axes = gp_Ax2( p0, normal, v01 );
383 catch ( Standard_Failure ) {
386 double factor = 1.0, modulus = normal.Modulus();
387 if ( modulus < 1e-2 )
388 factor = 1. / sqrt( modulus );
389 for ( size_t i = 0; i < nbNodes; ++i )
391 gp_XYZ p = _pv[i]._nxyz - p0;
392 _pv[i]._xy.SetX( axes.XDirection().XYZ() * p * factor);
393 _pv[i]._xy.SetY( axes.YDirection().XYZ() * p * factor );
396 // compute minimal triangle area
401 for ( size_t i = 0; i < nbNodes; ++i )
402 sumArea += _pv[i].TriaArea();
403 const double minArea = 1e-6 * sumArea / ( nbNodes - 2 );
405 // in a loop, find triangles with positive area and having no vertices inside
406 int iN = 0, nbTria = nbNodes - 2;
407 nodes.resize( nbTria * 3 );
408 _nodeIndex.resize( nbTria * 3 );
409 PolyVertex* v = &_pv[0], *vi;
410 int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
411 while ( nbBadTria < nbVertices )
413 if (( isGoodTria = v->TriaArea() > minArea ))
415 for ( vi = v->_next->_next;
419 if ( v->IsInsideTria( vi ))
422 isGoodTria = ( vi == v->_prev );
426 v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
429 if ( --nbVertices == 3 )
431 // last triangle remains
432 v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
434 _optimizer->optimize( nodes, _pv, _nodeIndex );
446 // the polygon is invalid; add triangles with positive area
448 while ( nbBadTria < nbVertices )
450 isGoodTria = v->TriaArea() > minArea;
453 v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
456 if ( --nbVertices == 3 )
458 // last triangle remains
459 v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
471 // add all the rest triangles
472 while ( nbVertices >= 3 )
474 v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
484 //================================================================================
488 //================================================================================
490 Triangulate::Triangulate( bool optimize ): _optimizer(0)
494 _optimizer = new Optimizer;
497 //================================================================================
501 //================================================================================
503 Triangulate::~Triangulate()
510 //================================================================================
512 * \brief Return nb triangles in a decomposed mesh face
513 * \retval int - number of triangles
515 //================================================================================
517 int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
519 // WARNING: counting triangles must be coherent with GetTriangles()
520 switch ( face->GetEntityType() )
522 case SMDSEntity_BiQuad_Triangle:
523 case SMDSEntity_BiQuad_Quadrangle:
524 return face->NbNodes() - 1;
525 // case SMDSEntity_Triangle:
526 // case SMDSEntity_Quad_Triangle:
527 // case SMDSEntity_Quadrangle:
528 // case SMDSEntity_Quad_Quadrangle:
529 // case SMDSEntity_Polygon:
530 // case SMDSEntity_Quad_Polygon:
532 return face->NbNodes() - 2;
537 //================================================================================
539 * \brief Decompose a mesh face into triangles
540 * \retval int - number of triangles
542 //================================================================================
544 int Triangulate::GetTriangles( const SMDS_MeshElement* face,
545 std::vector< const SMDS_MeshNode*>& nodes)
547 if ( face->GetType() != SMDSAbs_Face )
550 // WARNING: decomposing into triangles must be coherent with getNbTriangles()
551 int nbTria, i = 0, nbNodes = face->NbNodes();
552 SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
553 nodes.resize( nbNodes * 3 );
554 nodes[ i++ ] = nIt->next();
555 nodes[ i++ ] = nIt->next();
557 const SMDSAbs_EntityType type = face->GetEntityType();
560 case SMDSEntity_BiQuad_Triangle:
561 case SMDSEntity_BiQuad_Quadrangle:
563 nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
564 nodes[ i++ ] = face->GetNode( nbTria );
565 for ( i = 3; i < 3*(nbTria-1); i += 3 )
567 nodes[ i+0 ] = nodes[ i-2 ];
568 nodes[ i+1 ] = nIt->next();
569 nodes[ i+2 ] = nodes[ 2 ];
571 nodes[ i+0 ] = nodes[ i-2 ];
572 nodes[ i+1 ] = nodes[ 0 ];
573 nodes[ i+2 ] = nodes[ 2 ];
576 case SMDSEntity_Triangle:
579 nodes[ i++ ] = nIt->next();
584 // case SMDSEntity_Quad_Triangle:
585 // case SMDSEntity_Quadrangle:
586 // case SMDSEntity_Quad_Quadrangle:
587 // case SMDSEntity_Polygon:
588 // case SMDSEntity_Quad_Polygon:
590 nbTria = nbNodes - 2;
591 while ( nIt->more() )
592 nodes[ i++ ] = nIt->next();
594 if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
596 nIt = face->interlacedNodesIterator();
597 nodes[ 0 ] = nIt->next();
598 nodes[ 1 ] = nIt->next();
599 nodes[ 2 ] = nIt->next();
600 for ( i = 3; i < 3*nbTria; i += 3 )
602 nodes[ i+0 ] = nodes[ 0 ];
603 nodes[ i+1 ] = nodes[ i-1 ];
604 nodes[ i+2 ] = nIt->next();