src/SMESHUtils/SMESH_Triangulate.cxx

   1 // Copyright (C) 2007-2022  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 }