]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESHUtils/SMESH_Triangulate.cxx
Salome HOME
Merge branch 'master' into gni/adaptation
[modules/smesh.git] / src / SMESHUtils / SMESH_Triangulate.cxx
1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22 // File      : SMESH_Triangulate.cxx
23 // Created   : Thu Jan 18 18:00:13 2018
24 // Author    : Edward AGAPOV (eap)
25
26 // Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx
27
28 #include "SMESH_MeshAlgos.hxx"
29
30 #include <Standard_ErrorHandler.hxx>
31 #include <Standard_Failure.hxx>
32 #include <gp_Ax2.hxx>
33
34 #include <boost/container/flat_set.hpp>
35
36 using namespace SMESH_MeshAlgos;
37
38 namespace
39 {
40   struct Node // node of a triangle
41   {
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]
44
45     //! return node index within the node array
46     size_t Index() const { return  _triaIndex + _nodeIndex; }
47
48     //! return local 3-d index [0-2]
49     static size_t ThirdIndex( size_t i1, size_t i2 )
50     {
51       size_t i3 = ( i2 + 1 ) % 3;
52       if ( i3 == i1 )
53         i3 = ( i2 + 2 ) % 3;
54       return i3;
55     }
56     //! return 3-d node index within the node array
57     static size_t ThirdIndex( const Node& n1, const Node& n2 )
58     {
59       return n1._triaIndex + ThirdIndex( n1._nodeIndex, n2._nodeIndex );
60     }
61     bool operator<(const Node& other) const { return _triaIndex < other._triaIndex; }
62   };
63   typedef boost::container::flat_set< Node >                 TriaNodeSet;
64
65 }
66 /*!
67  * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
68  */
69 struct Triangulate::PolyVertex
70 {
71   SMESH_NodeXYZ _nxyz;
72   size_t        _index;
73   gp_XY         _xy;
74   PolyVertex*   _prev;
75   PolyVertex*   _next;
76
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 );
81   PolyVertex* Delete();
82
83   // compare PolyVertex'es by node
84   bool operator()(const PolyVertex* a, const PolyVertex* b) const
85   {
86     return ( a->_nxyz.Node() <  b->_nxyz.Node() );
87   }
88   // set of PolyVertex sorted by mesh node
89   typedef boost::container::flat_set< PolyVertex*, PolyVertex > PVSet;
90 };
91
92 struct Triangulate::Data
93 {
94   std::vector< PolyVertex > _pv;
95   std::vector< size_t >     _nodeIndex;
96   PolyVertex::PVSet         _uniqueNodePV;
97 };
98
99 struct Triangulate::Optimizer
100 {
101   std::vector< TriaNodeSet > _nodeUsage; // inclusions of a node in triangles
102
103   //================================================================================
104   /*!
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
109    */
110   //================================================================================
111
112   void optimize( std::vector< const SMDS_MeshNode*>& nodes,
113                  std::vector< PolyVertex > &         points,
114                  std::vector< size_t > &             nodeIndices)
115   {
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
119     {
120       _nodeUsage[ i ].clear();
121     }
122     for ( size_t i = 0, iTria = 0; i < nodeIndices.size(); ++iTria )
123     {
124       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 0 });
125       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 1 });
126       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 2 });
127     }
128
129     // optimization
130     for ( size_t iTria = 0; iTria < nodeIndices.size(); iTria += 3 )
131     {
132       double badness1 = computeBadness( nodeIndices[ iTria + 0 ],
133                                         nodeIndices[ iTria + 1 ],
134                                         nodeIndices[ iTria + 2 ],
135                                         points );
136       for ( size_t i = 0; i < 3; ++i ) // loop on triangle edges to find a neighbor triangle
137       {
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 ||
145              usage2.size() < 2 )
146           continue;
147
148         // look for another triangle using two nodes
149         TriaNodeSet::iterator usIt1 = usage1.begin();
150         for ( ; usIt1 != usage1.end(); ++usIt1 )
151         {
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
157
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 ];
162
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 );
166
167           if ( Max( badness1, badness2 ) < Max( badness3, badness4 ))
168             continue;
169
170           // swap edge by modifying nodeIndices
171
172           nodeIndices[ i2 ] = ind4;
173           _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria });
174           _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
175
176           i1 = usIt1->Index();
177           nodeIndices[ i1 ] = ind3;
178           _nodeUsage[ ind3 ].insert( *usIt1 );
179           _nodeUsage[ ind1 ].erase ( *usIt1 );
180
181           --i; // to re-check a current edge
182           badness1 = badness3;
183           break;
184         }
185       }
186     }
187
188     // update nodes by updated nodeIndices
189     for ( size_t i = 0; i < nodeIndices.size(); ++i )
190       nodes[ i ] = points[ nodeIndices[ i ]]._nxyz.Node();
191
192     return;
193   }
194
195   //================================================================================
196   /*!
197    * \brief Return 1./area. Initially: max cos^2 of triangle angles
198    */
199   //================================================================================
200
201   double computeBadness( size_t i1, size_t i2, size_t i3,
202                          std::vector< PolyVertex > & points,
203                          bool                        checkArea = false )
204   {
205     if ( checkArea )
206     {
207       points[ i2 ]._prev = & points[ i1 ];
208       points[ i2 ]._next = & points[ i3 ];
209       double a = points[ i2 ].TriaArea();
210       // if ( a < 0 )
211       //   return std::numeric_limits<double>::max();
212       // return 1. / a;
213
214       if ( a < 0 )
215         return 2;
216     }
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,
221                      p3 - p2,
222                      p1 - p3 };
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() )
229       return 2;
230
231     double maxCos2 = 0;
232     for ( int i = 0; i < 3; ++i )
233     {
234       int i2 = ( i+1 ) % 3;
235       double dot = -vec[ i ] * vec[ i2 ];
236       if ( dot > 0 )
237         maxCos2 = Max( maxCos2, dot * dot / len[ i ] / len[ i2 ] );
238     }
239     return maxCos2;
240   }
241 };
242
243 //================================================================================
244 /*!
245  * \brief Initialization
246  */
247 //================================================================================
248
249 void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
250                                               PolyVertex&          v,
251                                               size_t               index )
252 {
253   _nxyz.Set( n );
254   _next = &v;
255   v._prev = this;
256   _index = index;
257 }
258 //================================================================================
259 /*!
260  * \brief Remove self from a polygon
261  */
262 //================================================================================
263
264 Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
265 {
266   _prev->_next = _next;
267   _next->_prev = _prev;
268   return _next;
269 }
270
271 //================================================================================
272 /*!
273  * \brief Return nodes of a triangle
274  */
275 //================================================================================
276
277   void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes,
278                                               size_t*               nodeIndices) const
279 {
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;
286 }
287
288 //================================================================================
289 /*!
290  * \brief Compute triangle area
291  */
292 //================================================================================
293
294 inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
295 {
296   gp_XY vPrev = xy0 - xy1;
297   gp_XY vNext = xy2 - xy1;
298   return vNext ^ vPrev;
299 }
300
301 //================================================================================
302 /*!
303  * \brief Compute triangle area
304  */
305 //================================================================================
306
307 double Triangulate::PolyVertex::TriaArea() const
308 {
309   return Area( _prev->_xy, this->_xy, _next->_xy );
310 }
311
312 //================================================================================
313 /*!
314  * \brief Check if a vertex is inside a triangle
315  */
316 //================================================================================
317
318 bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
319 {
320   if ( this ->_nxyz == v->_nxyz ||
321        _prev->_nxyz == v->_nxyz ||
322        _next->_nxyz == v->_nxyz )
323     return false;
324
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 &&
330           ( t ^ n ) >= tol &&
331           ( n ^ p ) >= tol );
332   // return ( Area( _prev, this, v ) > 0 &&
333   //          Area( this, _next, v ) > 0 &&
334   //          Area( _next, _prev, v ) > 0 );
335 }
336
337 //================================================================================
338 /*!
339  * \brief Triangulate a polygon. Assure correct orientation for concave polygons
340  */
341 //================================================================================
342
343 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
344                                const size_t                        nbNodes)
345 {
346   std::vector< PolyVertex >&    _pv = _data->_pv;
347   std::vector< size_t >& _nodeIndex = _data->_nodeIndex;
348   PolyVertex::PVSet&  _uniqueNodePV = _data->_uniqueNodePV;
349
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 );
355
356   // assure correctness of PolyVertex::_index as a node can encounter more than once
357   // within a polygon boundary
358   if ( _optimizer && nbNodes > 4 )
359   {
360     _uniqueNodePV.clear();
361     for ( size_t i = 0; i < nbNodes; ++i )
362     {
363       PolyVertex::PVSet::iterator pv = _uniqueNodePV.insert( &_pv[i] ).first;
364       _pv[i]._index = (*pv)->_index;
365     }
366   }
367
368   // get a polygon normal
369   gp_XYZ normal(0,0,0), p0,v01,v02;
370   p0  = _pv[0]._nxyz;
371   v01 = _pv[1]._nxyz - p0;
372   for ( size_t i = 2; i < nbNodes; ++i )
373   {
374     v02 = _pv[i]._nxyz - p0;
375     normal += v01 ^ v02;
376     v01 = v02;
377   }
378   // project nodes to the found plane
379   gp_Ax2 axes;
380   try {
381     axes = gp_Ax2( p0, normal, v01 );
382   }
383   catch ( Standard_Failure ) {
384     return false;
385   }
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 )
390   {
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 );
394   }
395
396   // compute minimal triangle area
397   double sumArea = 0;
398   if ( factor == 1.0 )
399     sumArea = modulus;
400   else
401     for ( size_t i = 0; i < nbNodes; ++i )
402       sumArea += _pv[i].TriaArea();
403   const double minArea = 1e-6 * sumArea / ( nbNodes - 2 );
404
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 )
412   {
413     if (( isGoodTria = v->TriaArea() > minArea ))
414     {
415       for ( vi = v->_next->_next;
416             vi != v->_prev;
417             vi = vi->_next )
418       {
419         if ( v->IsInsideTria( vi ))
420           break;
421       }
422       isGoodTria = ( vi == v->_prev );
423     }
424     if ( isGoodTria )
425     {
426       v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
427       iN += 3;
428       v = v->Delete();
429       if ( --nbVertices == 3 )
430       {
431         // last triangle remains
432         v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
433         if ( _optimizer )
434           _optimizer->optimize( nodes, _pv, _nodeIndex );
435         return true;
436       }
437       nbBadTria = 0;
438     }
439     else
440     {
441       v = v->_next;
442       ++nbBadTria;
443     }
444   }
445
446   // the polygon is invalid; add triangles with positive area
447   nbBadTria = 0;
448   while ( nbBadTria < nbVertices )
449   {
450     isGoodTria = v->TriaArea() > minArea;
451     if ( isGoodTria )
452     {
453       v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
454       iN += 3;
455       v = v->Delete();
456       if ( --nbVertices == 3 )
457       {
458         // last triangle remains
459         v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
460         return true;
461       }
462       nbBadTria = 0;
463     }
464     else
465     {
466       v = v->_next;
467       ++nbBadTria;
468     }
469   }
470
471   // add all the rest triangles
472   while ( nbVertices >= 3 )
473   {
474     v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
475     iN += 3;
476     v = v->Delete();
477     --nbVertices;
478   }
479
480   return true;
481
482 } // triangulate()
483
484 //================================================================================
485 /*!
486  * \brief Constructor
487  */
488 //================================================================================
489
490 Triangulate::Triangulate( bool optimize ): _optimizer(0)
491 {
492   _data = new Data;
493   if ( optimize )
494     _optimizer = new Optimizer;
495 }
496
497 //================================================================================
498 /*!
499  * \brief Destructor
500  */
501 //================================================================================
502
503 Triangulate::~Triangulate()
504 {
505   delete _data;
506   delete _optimizer;
507   _optimizer = 0;
508 }
509
510 //================================================================================
511 /*!
512  * \brief Return nb triangles in a decomposed mesh face
513  *  \retval int - number of triangles
514  */
515 //================================================================================
516
517 int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
518 {
519   // WARNING: counting triangles must be coherent with GetTriangles()
520   switch ( face->GetEntityType() )
521   {
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:
531   default:
532     return face->NbNodes() - 2;
533   }
534   return 0;
535 }
536
537 //================================================================================
538 /*!
539  * \brief Decompose a mesh face into triangles
540  *  \retval int - number of triangles
541  */
542 //================================================================================
543
544 int Triangulate::GetTriangles( const SMDS_MeshElement*             face,
545                                std::vector< const SMDS_MeshNode*>& nodes)
546 {
547   if ( face->GetType() != SMDSAbs_Face )
548     return 0;
549
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();
556
557   const SMDSAbs_EntityType type = face->GetEntityType();
558   switch ( type )
559   {
560   case SMDSEntity_BiQuad_Triangle:
561   case SMDSEntity_BiQuad_Quadrangle:
562
563     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
564     nodes[ i++ ] = face->GetNode( nbTria );
565     for ( i = 3; i < 3*(nbTria-1); i += 3 )
566     {
567       nodes[ i+0 ] = nodes[ i-2 ];
568       nodes[ i+1 ] = nIt->next();
569       nodes[ i+2 ] = nodes[ 2 ];
570     }
571     nodes[ i+0 ] = nodes[ i-2 ];
572     nodes[ i+1 ] = nodes[ 0 ];
573     nodes[ i+2 ] = nodes[ 2 ];
574     break;
575
576   case SMDSEntity_Triangle:
577
578     nbTria = 1;
579     nodes[ i++ ] = nIt->next();
580     break;
581
582   default:
583
584     // case SMDSEntity_Quad_Triangle:
585     // case SMDSEntity_Quadrangle:
586     // case SMDSEntity_Quad_Quadrangle:
587     // case SMDSEntity_Polygon:
588     // case SMDSEntity_Quad_Polygon:
589
590     nbTria = nbNodes - 2;
591     while ( nIt->more() )
592       nodes[ i++ ] = nIt->next();
593
594     if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
595     {
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 )
601       {
602         nodes[ i+0 ] = nodes[ 0 ];
603         nodes[ i+1 ] = nodes[ i-1 ];
604         nodes[ i+2 ] = nIt->next();
605       }
606     }
607   }
608
609   return nbTria;
610 }