src/SMESHUtils/SMESH_Triangulate.cxx

   1 // Copyright (C) 2007-2016  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 > TNodeSet;
  64
  65 }
  66
  67 struct Triangulate::Optimizer
  68 {
  69   std::vector< TNodeSet > _nodeUsage; // inclusions of a node in triangles
  70
  71   //================================================================================
  72   /*!
  73    * \brief Optimize triangles by edge swapping
  74    *  \param [inout] nodes - polygon triangulation, i.e. connectivity of all triangles to optimize
  75    *  \param [in] points - coordinates of nodes of the input polygon
  76    *  \param [in] nodeIndices - indices of triangulation nodes within the input polygon
  77    */
  78   //================================================================================
  79
  80   void optimize( std::vector< const SMDS_MeshNode*>& nodes,
  81                  std::vector< PolyVertex > &         points,
  82                  std::vector< size_t > &             nodeIndices)
  83   {
  84     // for each node of the polygon, remember triangles using it
  85     _nodeUsage.resize( points.size() );
  86     for ( size_t i = 0; i < points.size(); ++i ) // clear old data
  87     {
  88       _nodeUsage[ i ].clear();
  89     }
  90     for ( size_t i = 0, iTria = 0; i < nodeIndices.size(); ++iTria )
  91     {
  92       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 0 });
  93       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 1 });
  94       _nodeUsage[ nodeIndices[ i++ ]].insert({ iTria * 3, 2 });
  95     }
  96
  97     // optimization
  98     for ( size_t iTria = 0; iTria < nodeIndices.size(); iTria += 3 )
  99     {
 100       double badness1 = computeBadness( nodeIndices[ iTria + 0 ],
 101                                         nodeIndices[ iTria + 1 ],
 102                                         nodeIndices[ iTria + 2 ],
 103                                         points );
 104       for ( size_t i = 0; i < 3; ++i ) // loop on triangle edges to find a neighbor triangle
 105       {
 106         size_t i1 = iTria + i; // node index in nodeIndices
 107         size_t i2 = iTria + ( i + 1 ) % 3;
 108         size_t ind1 = nodeIndices[ i1 ]; // node index in points
 109         size_t ind2 = nodeIndices[ i2 ];
 110         TNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
 111         TNodeSet & usage2 = _nodeUsage[ ind2 ];
 112         if ( usage1.size() < 2 ||
 113              usage2.size() < 2 )
 114           continue;
 115
 116         // look for another triangle using two nodes
 117         TNodeSet::iterator usIt1 = usage1.begin();
 118         for ( ; usIt1 != usage1.end(); ++usIt1 )
 119         {
 120           if ( usIt1->_triaIndex == iTria )
 121             continue; // current triangle
 122           TNodeSet::iterator usIt2 = usage2.find( *usIt1 );
 123           if ( usIt2 == usage2.end() )
 124             continue; // no common _triaIndex in two usages
 125
 126           size_t i3 = iTria + ( i + 2 ) % 3;
 127           size_t i4 = Node::ThirdIndex( *usIt1, *usIt2 ); // 4th node of quadrangle
 128           size_t ind3 = nodeIndices[ i3 ];
 129           size_t ind4 = nodeIndices[ i4 ];
 130
 131           double badness2 = computeBadness( ind2, ind1, ind4, points );
 132           double badness3 = computeBadness( ind1, ind4, ind3, points, /*checkArea=*/true );
 133           double badness4 = computeBadness( ind2, ind3, ind4, points, /*checkArea=*/true );
 134
 135           if ( Max( badness1, badness2 ) < Max( badness3, badness4 ))
 136             continue;
 137
 138           // swap edge by modifying nodeIndices
 139
 140           nodeIndices[ i2 ] = ind4;
 141           _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
 142           _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria });
 143
 144           i1 = usIt1->Index();
 145           nodeIndices[ i1 ] = ind3;
 146           _nodeUsage[ ind1 ].erase ( *usIt1 );
 147           _nodeUsage[ ind3 ].insert( *usIt1 );
 148
 149           --i; // to re-check a current edge
 150           badness1 = badness3;
 151           break;
 152         }
 153       }
 154     }
 155
 156     // update nodes by updated nodeIndices
 157     for ( size_t i = 0; i < nodeIndices.size(); ++i )
 158       nodes[ i ] = points[ nodeIndices[ i ]]._nxyz.Node();
 159
 160     return;
 161   }
 162
 163   //================================================================================
 164   /*!
 165    * \brief Return 1./area. Initially: max cos^2 of triangle angles
 166    */
 167   //================================================================================
 168
 169   double computeBadness( size_t i1, size_t i2, size_t i3,
 170                          std::vector< PolyVertex > & points,
 171                          bool                        checkArea = false )
 172   {
 173     //if ( checkArea )
 174     {
 175       points[ i2 ]._prev = & points[ i1 ];
 176       points[ i2 ]._next = & points[ i3 ];
 177       double a = points[ i2 ].TriaArea();
 178       if ( a < 0 )
 179         return std::numeric_limits<double>::max();
 180       return 1. / a;
 181
 182       if ( points[ i2 ].TriaArea() < 0 )
 183         return 2;
 184     }
 185     const gp_XY & p1 = points[ i1 ]._xy;
 186     const gp_XY & p2 = points[ i2 ]._xy;
 187     const gp_XY & p3 = points[ i3 ]._xy;
 188     gp_XY vec[3] = { p2 - p1,
 189                      p3 - p2,
 190                      p1 - p3 };
 191     double len[3] = { vec[0].SquareModulus(),
 192                       vec[1].SquareModulus(),
 193                       vec[2].SquareModulus() };
 194     if ( len[0] < gp::Resolution() ||
 195          len[1] < gp::Resolution() ||
 196          len[2] < gp::Resolution() )
 197       return 2;
 198
 199     double maxCos2 = 0;
 200     for ( int i = 0; i < 3; ++i )
 201     {
 202       int i2 = ( i+1 ) % 3;
 203       double dot = -vec[ i ] * vec[ i2 ];
 204       if ( dot > 0 )
 205         maxCos2 = Max( maxCos2, dot * dot / len[ i ] / len[ i2 ] );
 206     }
 207     return maxCos2;
 208   }
 209 };
 210
 211 //================================================================================
 212 /*!
 213  * \brief Initialization
 214  */
 215 //================================================================================
 216
 217 void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
 218                                               PolyVertex&          v,
 219                                               size_t               index )
 220 {
 221   _nxyz.Set( n );
 222   _next = &v;
 223   v._prev = this;
 224   _index = index;
 225 }
 226 //================================================================================
 227 /*!
 228  * \brief Remove self from a polygon
 229  */
 230 //================================================================================
 231
 232 Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
 233 {
 234   _prev->_next = _next;
 235   _next->_prev = _prev;
 236   return _next;
 237 }
 238
 239 //================================================================================
 240 /*!
 241  * \brief Return nodes of a triangle
 242  */
 243 //================================================================================
 244
 245   void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes,
 246                                               size_t*               nodeIndices) const
 247 {
 248   nodes[0] = _prev->_nxyz._node;
 249   nodes[1] =  this->_nxyz._node;
 250   nodes[2] = _next->_nxyz._node;
 251   nodeIndices[0] = _prev->_index;
 252   nodeIndices[1] =  this->_index;
 253   nodeIndices[2] = _next->_index;
 254 }
 255
 256 //================================================================================
 257 /*!
 258  * \brief Compute triangle area
 259  */
 260 //================================================================================
 261
 262 inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
 263 {
 264   gp_XY vPrev = xy0 - xy1;
 265   gp_XY vNext = xy2 - xy1;
 266   return vNext ^ vPrev;
 267 }
 268
 269 //================================================================================
 270 /*!
 271  * \brief Compute triangle area
 272  */
 273 //================================================================================
 274
 275 double Triangulate::PolyVertex::TriaArea() const
 276 {
 277   return Area( _prev->_xy, this->_xy, _next->_xy );
 278 }
 279
 280 //================================================================================
 281 /*!
 282  * \brief Check if a vertex is inside a triangle
 283  */
 284 //================================================================================
 285
 286 bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
 287 {
 288   if ( this ->_nxyz == v->_nxyz ||
 289        _prev->_nxyz == v->_nxyz ||
 290        _next->_nxyz == v->_nxyz )
 291     return false;
 292
 293   gp_XY p = _prev->_xy - v->_xy;
 294   gp_XY t =  this->_xy - v->_xy;
 295   gp_XY n = _next->_xy - v->_xy;
 296   const double tol = -1e-7;
 297   return (( p ^ t ) >= tol &&
 298           ( t ^ n ) >= tol &&
 299           ( n ^ p ) >= tol );
 300   // return ( Area( _prev, this, v ) > 0 &&
 301   //          Area( this, _next, v ) > 0 &&
 302   //          Area( _next, _prev, v ) > 0 );
 303 }
 304
 305 //================================================================================
 306 /*!
 307  * \brief Triangulate a polygon. Assure correct orientation for concave polygons
 308  */
 309 //================================================================================
 310
 311 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
 312                                const size_t                        nbNodes)
 313 {
 314   // connect nodes into a ring
 315   _pv.resize( nbNodes );
 316   for ( size_t i = 1; i < nbNodes; ++i )
 317     _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], i-1 );
 318   _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0], nbNodes-1 );
 319
 320   // get a polygon normal
 321   gp_XYZ normal(0,0,0), p0,v01,v02;
 322   p0  = _pv[0]._nxyz;
 323   v01 = _pv[1]._nxyz - p0;
 324   for ( size_t i = 2; i < nbNodes; ++i )
 325   {
 326     v02 = _pv[i]._nxyz - p0;
 327     normal += v01 ^ v02;
 328     v01 = v02;
 329   }
 330   // project nodes to the found plane
 331   gp_Ax2 axes;
 332   try {
 333     axes = gp_Ax2( p0, normal, v01 );
 334   }
 335   catch ( Standard_Failure ) {
 336     return false;
 337   }
 338   for ( size_t i = 0; i < nbNodes; ++i )
 339   {
 340     gp_XYZ p = _pv[i]._nxyz - p0;
 341     _pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
 342     _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
 343   }
 344
 345   // in a loop, find triangles with positive area and having no vertices inside
 346   int iN = 0, nbTria = nbNodes - 2;
 347   nodes.resize( nbTria * 3 );
 348   _nodeIndex.resize( nbTria * 3 );
 349   const double minArea = 1e-6;
 350   PolyVertex* v = &_pv[0], *vi;
 351   int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
 352   while ( nbBadTria < nbVertices )
 353   {
 354     if (( isGoodTria = v->TriaArea() > minArea ))
 355     {
 356       for ( vi = v->_next->_next;
 357             vi != v->_prev;
 358             vi = vi->_next )
 359       {
 360         if ( v->IsInsideTria( vi ))
 361           break;
 362       }
 363       isGoodTria = ( vi == v->_prev );
 364     }
 365     if ( isGoodTria )
 366     {
 367       v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
 368       iN += 3;
 369       v = v->Delete();
 370       if ( --nbVertices == 3 )
 371       {
 372         // last triangle remains
 373         v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
 374         if ( _optimizer )
 375           _optimizer->optimize( nodes, _pv, _nodeIndex );
 376         return true;
 377       }
 378       nbBadTria = 0;
 379     }
 380     else
 381     {
 382       v = v->_next;
 383       ++nbBadTria;
 384     }
 385   }
 386
 387   // the polygon is invalid; add triangles with positive area
 388   nbBadTria = 0;
 389   while ( nbBadTria < nbVertices )
 390   {
 391     isGoodTria = v->TriaArea() > minArea;
 392     if ( isGoodTria )
 393     {
 394       v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
 395       iN += 3;
 396       v = v->Delete();
 397       if ( --nbVertices == 3 )
 398       {
 399         // last triangle remains
 400         v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
 401         return true;
 402       }
 403       nbBadTria = 0;
 404     }
 405     else
 406     {
 407       v = v->_next;
 408       ++nbBadTria;
 409     }
 410   }
 411
 412   // add all the rest triangles
 413   while ( nbVertices >= 3 )
 414   {
 415     v->GetTriaNodes( &nodes[ iN ], &_nodeIndex[ iN ] );
 416     iN += 3;
 417     v = v->Delete();
 418     --nbVertices;
 419   }
 420
 421   return true;
 422
 423 } // triangulate()
 424
 425 //================================================================================
 426 /*!
 427  * \brief Constructor
 428  */
 429 //================================================================================
 430
 431 Triangulate::Triangulate( bool optimize ): _optimizer(0)
 432 {
 433   if ( optimize )
 434     _optimizer = new Optimizer;
 435 }
 436
 437 //================================================================================
 438 /*!
 439  * \brief Destructor
 440  */
 441 //================================================================================
 442
 443 Triangulate::~Triangulate()
 444 {
 445   delete _optimizer;
 446   _optimizer = 0;
 447 }
 448
 449 //================================================================================
 450 /*!
 451  * \brief Return nb triangles in a decomposed mesh face
 452  *  \retval int - number of triangles
 453  */
 454 //================================================================================
 455
 456 int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
 457 {
 458   // WARNING: counting triangles must be coherent with GetTriangles()
 459   switch ( face->GetEntityType() )
 460   {
 461   case SMDSEntity_BiQuad_Triangle:
 462   case SMDSEntity_BiQuad_Quadrangle:
 463     return face->NbNodes() - 1;
 464     // case SMDSEntity_Triangle:
 465     // case SMDSEntity_Quad_Triangle:
 466     // case SMDSEntity_Quadrangle:
 467     // case SMDSEntity_Quad_Quadrangle:
 468     // case SMDSEntity_Polygon:
 469     // case SMDSEntity_Quad_Polygon:
 470   default:
 471     return face->NbNodes() - 2;
 472   }
 473   return 0;
 474 }
 475
 476 //================================================================================
 477 /*!
 478  * \brief Decompose a mesh face into triangles
 479  *  \retval int - number of triangles
 480  */
 481 //================================================================================
 482
 483 int Triangulate::GetTriangles( const SMDS_MeshElement*             face,
 484                                std::vector< const SMDS_MeshNode*>& nodes)
 485 {
 486   if ( face->GetType() != SMDSAbs_Face )
 487     return 0;
 488
 489   // WARNING: decomposing into triangles must be coherent with getNbTriangles()
 490   int nbTria, i = 0, nbNodes = face->NbNodes();
 491   SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
 492   nodes.resize( nbNodes * 3 );
 493   nodes[ i++ ] = nIt->next();
 494   nodes[ i++ ] = nIt->next();
 495
 496   const SMDSAbs_EntityType type = face->GetEntityType();
 497   switch ( type )
 498   {
 499   case SMDSEntity_BiQuad_Triangle:
 500   case SMDSEntity_BiQuad_Quadrangle:
 501
 502     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
 503     nodes[ i++ ] = face->GetNode( nbTria );
 504     for ( i = 3; i < 3*(nbTria-1); i += 3 )
 505     {
 506       nodes[ i+0 ] = nodes[ i-2 ];
 507       nodes[ i+1 ] = nIt->next();
 508       nodes[ i+2 ] = nodes[ 2 ];
 509     }
 510     nodes[ i+0 ] = nodes[ i-2 ];
 511     nodes[ i+1 ] = nodes[ 0 ];
 512     nodes[ i+2 ] = nodes[ 2 ];
 513     break;
 514
 515   case SMDSEntity_Triangle:
 516
 517     nbTria = 1;
 518     nodes[ i++ ] = nIt->next();
 519     break;
 520
 521   default:
 522
 523     // case SMDSEntity_Quad_Triangle:
 524     // case SMDSEntity_Quadrangle:
 525     // case SMDSEntity_Quad_Quadrangle:
 526     // case SMDSEntity_Polygon:
 527     // case SMDSEntity_Quad_Polygon:
 528
 529     nbTria = nbNodes - 2;
 530     while ( nIt->more() )
 531       nodes[ i++ ] = nIt->next();
 532
 533     if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
 534     {
 535       nIt = face->interlacedNodesIterator();
 536       nodes[ 0 ] = nIt->next();
 537       nodes[ 1 ] = nIt->next();
 538       nodes[ 2 ] = nIt->next();
 539       for ( i = 3; i < 3*nbTria; i += 3 )
 540       {
 541         nodes[ i+0 ] = nodes[ 0 ];
 542         nodes[ i+1 ] = nodes[ i-1 ];
 543         nodes[ i+2 ] = nIt->next();
 544       }
 545     }
 546   }
 547
 548   return nbTria;
 549 }