lib/cpm/DataStructures/QuadEdgeMesh.hxx

   1 #ifndef __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__
   2 #define __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__
   3
   4 // -------------------------------------------------------------------------
   5 template< typename P, unsigned int D, typename T >
   6 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
   7 TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
   8 FindEdge( const PointIdentifier& a ) const
   9 {
  10   if( a < this->m_OnextRings.size( ) )
  11     return( this->m_OnextRings[ a ] );
  12   else
  13     return( NULL );
  14 }
  15
  16 // -------------------------------------------------------------------------
  17 template< typename P, unsigned int D, typename T >
  18 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  19 TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  20 FindEdge( const PointIdentifier& a, const PointIdentifier& b ) const
  21 {
  22   TPrimalEdge* found = NULL;
  23   TPrimalEdge* e = this->FindEdge( a );
  24   if( e != NULL )
  25   {
  26     typename TPrimalEdge::Iterator eIt = e->BeginOnext( );
  27     for( ; eIt != e->EndOnext( ) && found == NULL; eIt++ )
  28       if( ( *eIt )->GetDestination( ) == b )
  29         found = *eIt;
  30
  31   } // fi
  32   return( found );
  33 }
  34
  35 // -------------------------------------------------------------------------
  36 template< typename P, unsigned int D, typename T >
  37 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  38 TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  39 FindEntryEdge( const CellIdentifier& i ) const
  40 {
  41   if( i < this->GetNumberOfCells( ) )
  42   {
  43     CellAutoPointer c;
  44     this->GetCell( i, c );
  45
  46     TQuadEdgeCell* ec = dynamic_cast< TQuadEdgeCell* >( c.GetPointer( ) );
  47     if( ec != NULL )
  48       return( ec->GetEntryPrimalEdge( ) );
  49     else
  50       return( NULL );
  51   }
  52   else
  53     return( NULL );
  54 }
  55
  56 // -------------------------------------------------------------------------
  57 template< typename P, unsigned int D, typename T >
  58 const typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  59 CntNormals& cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  60 GetPointNormalsContainer( ) const
  61 {
  62   return( this->m_PointNormals );
  63 }
  64
  65 // -------------------------------------------------------------------------
  66 template< typename P, unsigned int D, typename T >
  67 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  68 NormalsIterator cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  69 BeginPointNormals( ) const
  70 {
  71   return( this->m_PointNormals.begin( ) );
  72 }
  73
  74 // -------------------------------------------------------------------------
  75 template< typename P, unsigned int D, typename T >
  76 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  77 NormalsIterator cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  78 EndPointNormals( ) const
  79 {
  80   return( this->m_PointNormals.end( ) );
  81 }
  82
  83 // -------------------------------------------------------------------------
  84 template< typename P, unsigned int D, typename T >
  85 const typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  86 VectorType& cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  87 GetPointNormal( const PointIdentifier& id ) const
  88 {
  89   static const VectorType zero( TScalar( 0 ) );
  90   if( id < this->m_PointNormals.size( ) )
  91     return( this->m_PointNormals[ id ] );
  92   else
  93     return( zero );
  94 }
  95
  96 // -------------------------------------------------------------------------
  97 template< typename P, unsigned int D, typename T >
  98 bool cpm::DataStructures::QuadEdgeMesh< P, D, T >::
  99 RequestedRegionIsOutsideOfTheBufferedRegion( )
 100 {
 101   // Mesh regions don't have meaning with QuadEdges, this is important
 102   // in order to guarantee correct pipeline execution
 103   return( false );
 104 }
 105
 106 // -------------------------------------------------------------------------
 107 template< typename P, unsigned int D, typename T >
 108 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 109 SetPoints( PointsContainer* points )
 110 {
 111   this->_ReleaseQuadEdgeObjects( );
 112   this->m_OnextRings.resize( points->Size( ), NULL );
 113   this->m_PointNormals.resize(
 114     points->Size( ), VectorType( TScalar( 0 ) )
 115     );
 116   this->Superclass::SetPoints( points );
 117 }
 118
 119 // -------------------------------------------------------------------------
 120 template< typename P, unsigned int D, typename T >
 121 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 122 SetPoint( PointIdentifier id , PointType point )
 123 {
 124   // If the point is brand new, add some non-existent edges and normals
 125   unsigned long N = this->GetNumberOfPoints( );
 126   N = ( N < id + 1 )? id + 1: N;
 127   if( this->m_OnextRings.size( ) < N )
 128   {
 129     this->m_OnextRings.resize( N, NULL );
 130     this->m_PointNormals.resize( N, VectorType( TScalar( 0 ) ) );
 131
 132   } // fi
 133
 134   // Ok pass it to itk
 135   this->Superclass::SetPoint( id, point );
 136 }
 137
 138 // -------------------------------------------------------------------------
 139 template< typename P, unsigned int D, typename T >
 140 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 141 Initialize( )
 142 {
 143   this->Superclass::Initialize( );
 144   this->_ReleaseQuadEdgeObjects( );
 145 }
 146
 147 // -------------------------------------------------------------------------
 148 template< typename P, unsigned int D, typename T >
 149 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 150 Graft( const itk::DataObject* data )
 151 {
 152   this->Superclass::Graft( data );
 153   const Self* mesh = dynamic_cast< const Self* >( data );
 154   if( mesh != NULL )
 155   {
 156     this->_ReleaseQuadEdgeObjects( );
 157     this->m_PrimalEdges    = mesh->m_PrimalEdges;
 158     this->m_DualEdges      = mesh->m_DualEdges;
 159     this->m_OnextRings     = mesh->m_OnextRings;
 160     this->m_PointNormals   = mesh->m_PointNormals;
 161
 162   } // fi
 163 }
 164
 165 // -------------------------------------------------------------------------
 166 template< typename P, unsigned int D, typename T >
 167 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 168 SetCells( CellsContainer* cells )
 169 {
 170   // Clear cells
 171   if( this->GetCells( ) != NULL )
 172     this->SetCells( CellsContainer::New( ) );
 173   this->_ReleaseQuadEdgeObjects( );
 174   this->m_OnextRings.resize( this->GetNumberOfPoints( ), NULL );
 175
 176   // Copy cells
 177   if( cells != NULL )
 178   {
 179     for( unsigned long cId = 0; cId < cells->Size( ); ++cId )
 180     {
 181       CellAutoPointer cell( cells->GetElement( cId ), false );
 182       this->SetCell( cId, cell );
 183
 184     } // rof
 185
 186   } // fi
 187 }
 188
 189 // -------------------------------------------------------------------------
 190 template< typename P, unsigned int D, typename T >
 191 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 192 SetCellData( CellDataContainer* data )
 193 {
 194   /* TODO
 195      No need for CellData in initial application:
 196      how should this be overridden?
 197   */
 198 }
 199
 200 // -------------------------------------------------------------------------
 201 template< typename P, unsigned int D, typename T >
 202 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 203 SetCell( CellIdentifier id, CellAutoPointer& ptr )
 204 {
 205   // Overwrite an existing cell?
 206   if( id < this->GetNumberOfCells( ) )
 207   {
 208     itkExceptionMacro(
 209       << "TODO: What do I do here? Trying to overwrite a cell!"
 210       );
 211
 212   } // fi
 213
 214   // Do not allow to add vertices or edges as cells for the sake of
 215   // correct mesh construction
 216   if( ptr->GetNumberOfPoints( ) < 3 )
 217   {
 218     itkExceptionMacro(
 219       << "Face have less than 3 points ("
 220       << ptr->GetNumberOfPoints( )
 221       << ")"
 222       );
 223
 224   } // fi
 225
 226   // Fill up geometrical prerrequisites
 227   if( !( this->_CheckPoints( ptr ) ) )
 228   {
 229     itkExceptionMacro( << "Points missing in the mesh" );
 230
 231   } // fi
 232
 233   // Get (or create) individual edges
 234   _TEdges edges;
 235   this->_ConstructEdges( edges, ptr );
 236
 237   // Configure new face's geometry
 238   edges[ 0 ]->SetLnextRingGeometry( id );
 239
 240   // Update normals
 241   typename _TEdges::const_iterator eIt = edges.begin( );
 242   for( ; eIt != edges.end( ); ++eIt )
 243     this->m_PointNormals[ ( *eIt )->GetOrigin( ) ] =
 244       this->_ComputePointNormal( *eIt );
 245
 246   // Save cell
 247   TQuadEdgeCell* qcell = new TQuadEdgeCell( edges[ 0 ] );
 248   CellAutoPointer cell;
 249   cell.TakeOwnership( qcell );
 250   this->Superclass::SetCell( id, cell );
 251 }
 252
 253 // -------------------------------------------------------------------------
 254 template< typename P, unsigned int D, typename T >
 255 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 256 SetCellData( CellIdentifier id, CellPixelType data )
 257 {
 258   /* TODO
 259      No need for CellData in initial application:
 260      how should this be overridden?
 261   */
 262 }
 263
 264 // -------------------------------------------------------------------------
 265 template< typename P, unsigned int D, typename T >
 266 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 267 BuildCellLinks( ) const
 268 {
 269   // This is no longer required when using QuadEdges
 270 }
 271
 272 // -------------------------------------------------------------------------
 273 template< typename P, unsigned int D, typename T >
 274 cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 275 QuadEdgeMesh( )
 276   : Superclass( )
 277 {
 278   this->_ReleaseQuadEdgeObjects( );
 279 }
 280
 281 // -------------------------------------------------------------------------
 282 template< typename P, unsigned int D, typename T >
 283 cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 284 ~QuadEdgeMesh( )
 285 {
 286   this->_ReleaseQuadEdgeObjects( );
 287 }
 288
 289 // -------------------------------------------------------------------------
 290 template< typename P, unsigned int D, typename T >
 291 bool cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 292 _CheckPoints( const CellAutoPointer& ptr ) const
 293 {
 294   const typename Superclass::PointsContainer* pnts = this->GetPoints( );
 295
 296   bool all_exist = true;
 297   typename CellType::PointIdConstIterator pIt = ptr->PointIdsBegin( );
 298   for( ; pIt != ptr->PointIdsEnd( ) && all_exist; ++pIt )
 299     all_exist &= pnts->IndexExists( *pIt );
 300   return( all_exist );
 301 }
 302
 303 // -------------------------------------------------------------------------
 304 template< typename P, unsigned int D, typename T >
 305 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 306 _DeletePoint( const PointIdentifier& pId )
 307 {
 308   PointIdentifier lastId = this->GetNumberOfPoints( ) - 1;
 309   if( lastId < pId )
 310     return;
 311
 312   // Modify point container
 313   PointType to_delete = this->GetPoint( pId );
 314   PointType last = this->GetPoint( lastId );
 315   this->Superclass::SetPoint( pId, last );
 316
 317   // Modify Onext ring and normals containers
 318   this->m_OnextRings[ pId ] = this->m_OnextRings[ lastId ];
 319   this->m_PointNormals[ pId ] = this->m_PointNormals[ lastId ];
 320   this->m_OnextRings.pop_back( );
 321   this->m_PointNormals.pop_back( );
 322
 323   // Update origins
 324   typename TPrimalEdge::Iterator oIt =
 325     this->m_OnextRings[ pId ]->BeginOnext( );
 326   for( ; oIt != this->m_OnextRings[ pId ]->EndOnext( ); ++oIt )
 327     ( *oIt )->SetOrigin( pId );
 328
 329   // Erase the point
 330   this->GetPoints( )->CastToSTLContainer( ).pop_back( );
 331 }
 332
 333 // -------------------------------------------------------------------------
 334 template< typename P, unsigned int D, typename T >
 335 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 336 TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 337 _CreateQuadEdge( const PointIdentifier& a, const PointIdentifier& b )
 338 {
 339   // Create brand new object
 340   typename TPrimalEdge::Pointer edge = TPrimalEdge::New( );
 341   edge->CreateRings( );
 342
 343   // Configure geometry
 344   edge->SetOrigin( a );
 345   edge->SetDestination( b );
 346
 347   // Keep trace of all reserved memory
 348   this->m_PrimalEdges.insert( edge );
 349   this->m_PrimalEdges.insert( edge->GetSym( ) );
 350   this->m_DualEdges.insert( edge->GetRot( ) );
 351   this->m_DualEdges.insert( edge->GetInvRot( ) );
 352
 353   return( edge );
 354 }
 355
 356 // -------------------------------------------------------------------------
 357 template< typename P, unsigned int D, typename T >
 358 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 359 _DeleteEdge( TPrimalEdge* edge )
 360 {
 361   this->m_DualEdges.erase( edge->GetInvRot( ) );
 362   this->m_PrimalEdges.erase( edge->GetSym( ) );
 363   this->m_DualEdges.erase( edge->GetRot( ) );
 364   this->m_PrimalEdges.erase( edge );
 365 }
 366
 367 // -------------------------------------------------------------------------
 368 template< typename P, unsigned int D, typename T >
 369 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 370 _DeleteFace( const CellIdentifier& f )
 371 {
 372   unsigned long cellId = this->GetNumberOfCells( ) - 1;
 373   if( cellId < f )
 374     return;
 375
 376   // Exchange desired face with last one
 377   CellAutoPointer to_delete, last;
 378   this->GetCell( f, to_delete );
 379   this->GetCell( cellId, last );
 380   this->Superclass::SetCell( f, last );
 381   this->Superclass::SetCell( cellId, to_delete );
 382
 383   // This will delete cell's memory when to_delete is no longer used
 384   to_delete.TakeOwnership( );
 385
 386   // Squeeze container
 387   this->GetCells( )->CastToSTLContainer( ).pop_back( );
 388 }
 389
 390 // -------------------------------------------------------------------------
 391 template< typename P, unsigned int D, typename T >
 392 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 393 _ConstructEdges( _TEdges& edges, const CellAutoPointer& ptr )
 394 {
 395   edges.clear( );
 396   typename CellType::PointIdConstIterator pIt = ptr->PointIdsBegin( );
 397   for( pIt = ptr->PointIdsBegin( ); pIt != ptr->PointIdsEnd( ); ++pIt )
 398   {
 399     typename CellType::PointIdConstIterator nIt = pIt;
 400     ++nIt;
 401     if( nIt == ptr->PointIdsEnd( ) )
 402       nIt = ptr->PointIdsBegin( );
 403
 404     TPrimalEdge* edge = this->FindEdge( *pIt, *nIt );
 405     if( edge == NULL )
 406     {
 407       edge = this->_CreateQuadEdge( *pIt, *nIt );
 408       if( this->m_OnextRings[ *pIt ] != NULL )
 409       {
 410         if( this->m_OnextRings[ *pIt ]->IsOriginInternal( ) )
 411         {
 412           itkExceptionMacro(
 413             << "Edge ["
 414             << *pIt << ", " << *nIt << "] is already internal."
 415             );
 416
 417         } // fi
 418
 419         // Put edge in Onext ring
 420         typename TPrimalEdge::Iterator eIt =
 421           this->m_OnextRings[ *pIt ]->BeginOnext( );
 422         while( eIt != this->m_OnextRings[ *pIt ]->EndOnext( ) )
 423         {
 424           if( !( ( *eIt )->IsLeftSet( ) ) )
 425           {
 426             TPrimalEdge::Splice( *eIt, edge );
 427             eIt = this->m_OnextRings[ *pIt ]->EndOnext( );
 428           }
 429           else
 430             ++eIt;
 431
 432         } // elihw
 433       }
 434       else
 435         this->m_OnextRings[ *pIt ] = edge;
 436     }
 437     else
 438     {
 439       if( edge->IsLeftSet( ) )
 440       {
 441         itkExceptionMacro(
 442           << "Edge ["
 443           << edge->GetOrigin( ) << ", "
 444           << edge->GetDestination( ) << "] already have a left face ("
 445           << edge->GetLeft( )
 446           << "). Face already exits (at least in part)."
 447           );
 448
 449       } // fi
 450
 451     } // fi
 452     edges.push_back( edge );
 453
 454   } // rof
 455
 456   // Reorder Onext rings
 457   typename _TEdges::iterator ecIt;
 458   for( ecIt = edges.begin( ); ecIt != edges.end( ); ++ecIt )
 459   {
 460     TPrimalEdge* e = *ecIt;
 461     TPrimalEdge* eOprev;
 462     if( ecIt != edges.begin( ) )
 463     {
 464       typename _TEdges::iterator eOprevIt = ecIt;
 465       --eOprevIt;
 466       eOprev = ( *eOprevIt )->GetSym( );
 467     }
 468     else
 469       eOprev = edges.back( )->GetSym( );
 470     TPrimalEdge::ReorderOnextRing( e, eOprev );
 471
 472   } // rof
 473 }
 474
 475 // -------------------------------------------------------------------------
 476 template< typename P, unsigned int D, typename T >
 477 typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 478 VectorType cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 479 _ComputePointNormal( const TPrimalEdge* e ) const
 480 {
 481   static const TScalar zero = TScalar( 0 );
 482
 483   VectorType n( zero );
 484   if( Superclass::PointDimension == 3 )
 485   {
 486     PointType p0 = this->GetPoint( e->GetOrigin( ) );
 487     typename TPrimalEdge::ConstIterator eIt = e->BeginOnext( );
 488     unsigned int count = 0;
 489     for( ; eIt != e->EndOnext( ); eIt++ )
 490     {
 491       if( ( *eIt )->IsLeftSet( ) )
 492       {
 493         typename TPrimalEdge::ConstIterator nIt = eIt;
 494         nIt++;
 495         if( nIt == e->EndOnext( ) )
 496           nIt = e->BeginOnext( );
 497
 498         n += itk::CrossProduct(
 499           this->GetPoint( ( *eIt )->GetDestination( ) ) - p0,
 500           this->GetPoint( ( *nIt )->GetDestination( ) ) - p0
 501           );
 502         count++;
 503
 504       } // fi
 505
 506     } // rof
 507     TScalar nn = n.GetNorm( );
 508     if( nn > zero && count > 0 )
 509       n /= nn * TScalar( count );
 510
 511   } // fi
 512   return( n );
 513 }
 514
 515 // -------------------------------------------------------------------------
 516 template< typename P, unsigned int D, typename T >
 517 void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
 518 _ReleaseQuadEdgeObjects( )
 519 {
 520   this->m_PrimalEdges.clear( );
 521   this->m_OnextRings.clear( );
 522   this->m_DualEdges.clear( );
 523   this->m_PointNormals.clear( );
 524 }
 525
 526 #endif // __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__
 527
 528 // eof - $RCSfile$