lib/cpm/DataStructures/SimplexMesh.hxx

   1 #ifndef __CPM__DATASTRUCTURES__SIMPLEXMESH__HXX__
   2 #define __CPM__DATASTRUCTURES__SIMPLEXMESH__HXX__
   3
   4 #include <cmath>
   5 #include <vnl/vnl_math.h>
   6 #include <itkMatrix.h>
   7
   8 // -------------------------------------------------------------------------
   9 template< typename P, unsigned int K, unsigned int D, typename T >
  10 bool cpm::DataStructures::SimplexMesh< P, K, D, T >::
  11 ComputeSimplexParameters(
  12   const PointIdentifier& pId,
  13   VectorType& A,
  14   VectorType& B,
  15   VectorType& C,
  16   VectorType& N,
  17   BarycenterVector& barycenter,
  18   ValueType& r,
  19   ValueType& d,
  20   ValueType& dp,
  21   ValueType& phi
  22   ) const
  23 {
  24   typedef itk::Matrix< ValueType, 4, 4 > TTetraMatrix;
  25
  26   static const ValueType _0 = ValueType( 0 );
  27   static const ValueType _1 = ValueType( 1 );
  28   static const ValueType _2 = ValueType( 2 );
  29   static const ValueType _3 = ValueType( 3 );
  30   static const ValueType _4 = ValueType( 4 );
  31
  32   // Get topological neighborhood
  33   const TPrimalEdge* eA = this->FindEdge( pId );
  34   if( eA == NULL )
  35     return( false );
  36   const TPrimalEdge* eB = eA->GetOnext( );
  37   if( eB == NULL )
  38     return( false );
  39   const TPrimalEdge* eC = eB->GetOnext( );
  40   if( eC == NULL )
  41     return( false );
  42
  43   // Get geometrical neighborhood
  44   VectorType Q;
  45   Q = this->GetPoint( eA->GetOrigin( ) ).GetVectorFromOrigin( );
  46   A = this->GetPoint( eA->GetDestination( ) ).GetVectorFromOrigin( );
  47   B = this->GetPoint( eB->GetDestination( ) ).GetVectorFromOrigin( );
  48   C = this->GetPoint( eC->GetDestination( ) ).GetVectorFromOrigin( );
  49
  50   // Auxiliary values
  51   VectorType BA = B - A;
  52   VectorType CA = C - A;
  53   VectorType BAxCA = itk::CrossProduct( BA, CA );
  54   ValueType nBAxCA = BAxCA.GetNorm( );
  55   if( !( _0 < nBAxCA ) )
  56     return( false );
  57
  58   // More auxiliary values
  59   VectorType QA = Q - A;
  60
  61   // Compute barycenter
  62   ValueType nBAxCA2 = nBAxCA * nBAxCA;
  63   barycenter[ 2 ] = ( itk::CrossProduct( BA, QA ) * BAxCA ) / nBAxCA2;
  64   barycenter[ 1 ] = ( itk::CrossProduct( QA, CA ) * BAxCA ) / nBAxCA2;
  65   barycenter[ 0 ] = _1 - barycenter[ 2 ] - barycenter[ 1 ];
  66
  67   // Compute Q projection onto neighboring plane
  68   VectorType Qp =
  69     ( A * barycenter[ 0 ] ) +
  70     ( B * barycenter[ 1 ] ) +
  71     ( C * barycenter[ 2 ] );
  72
  73   // Compute normal
  74   N = BAxCA / nBAxCA;
  75
  76   // Compute circumcircle center and radius (r)
  77   VectorType CC = BA * CA.GetSquaredNorm( );
  78   CC -= CA * BA.GetSquaredNorm( );
  79   CC = itk::CrossProduct( BAxCA, CC );
  80   CC /= _2 * nBAxCA2;
  81   r = CC.GetNorm( );
  82   CC += A;
  83
  84   // Barycenter distance
  85   d = ( Qp - CC ).GetNorm( );
  86
  87   // Desired barycenter distance
  88   dp = ( ( ( A + B + C ) / _3 ) - CC ).GetNorm( );
  89
  90   // Compute tetrahedron's circumsphere
  91   TTetraMatrix mA;
  92   for( unsigned int dId = 0; dId < 3; dId++ )
  93   {
  94     mA[ 0 ][ dId ] = A[ dId ];
  95     mA[ 1 ][ dId ] = B[ dId ];
  96     mA[ 2 ][ dId ] = C[ dId ];
  97     mA[ 3 ][ dId ] = Q[ dId ];
  98
  99   } // rof
 100   mA[ 0 ][ 3 ] = _1;
 101   mA[ 1 ][ 3 ] = _1;
 102   mA[ 2 ][ 3 ] = _1;
 103   mA[ 3 ][ 3 ] = _1;
 104
 105   ValueType dA = vnl_determinant( mA.GetVnlMatrix( ) );
 106   if( dA > _0 ) // Tetrahedron IS NOT planar
 107   {
 108     TTetraMatrix mD;
 109     mD[ 0 ][ 3 ] = _1;
 110     mD[ 1 ][ 3 ] = _1;
 111     mD[ 2 ][ 3 ] = _1;
 112     mD[ 3 ][ 3 ] = _1;
 113     mD[ 0 ][ 0 ] = A.GetSquaredNorm( );
 114     mD[ 1 ][ 0 ] = B.GetSquaredNorm( );
 115     mD[ 2 ][ 0 ] = C.GetSquaredNorm( );
 116     mD[ 3 ][ 0 ] = Q.GetSquaredNorm( );
 117
 118     mD[ 0 ][ 1 ] = A[ 1 ];
 119     mD[ 1 ][ 1 ] = B[ 1 ];
 120     mD[ 2 ][ 1 ] = C[ 1 ];
 121     mD[ 3 ][ 1 ] = Q[ 1 ];
 122     mD[ 0 ][ 2 ] = A[ 2 ];
 123     mD[ 1 ][ 2 ] = B[ 2 ];
 124     mD[ 2 ][ 2 ] = C[ 2 ];
 125     mD[ 3 ][ 2 ] = Q[ 2 ];
 126     ValueType Dx = vnl_determinant( mD.GetVnlMatrix( ) );
 127
 128     mD[ 0 ][ 1 ] = A[ 0 ];
 129     mD[ 1 ][ 1 ] = B[ 0 ];
 130     mD[ 2 ][ 1 ] = C[ 0 ];
 131     mD[ 3 ][ 1 ] = Q[ 0 ];
 132     mD[ 0 ][ 2 ] = A[ 2 ];
 133     mD[ 1 ][ 2 ] = B[ 2 ];
 134     mD[ 2 ][ 2 ] = C[ 2 ];
 135     mD[ 3 ][ 2 ] = Q[ 2 ];
 136     ValueType Dy = vnl_determinant( mD.GetVnlMatrix( ) );
 137
 138     mD[ 0 ][ 1 ] = A[ 0 ];
 139     mD[ 1 ][ 1 ] = B[ 0 ];
 140     mD[ 2 ][ 1 ] = C[ 0 ];
 141     mD[ 3 ][ 1 ] = Q[ 0 ];
 142     mD[ 0 ][ 2 ] = A[ 1 ];
 143     mD[ 1 ][ 2 ] = B[ 1 ];
 144     mD[ 2 ][ 2 ] = C[ 1 ];
 145     mD[ 3 ][ 2 ] = Q[ 1 ];
 146     ValueType Dz = vnl_determinant( mD.GetVnlMatrix( ) );
 147
 148     mD[ 0 ][ 2 ] = A[ 2 ];
 149     mD[ 1 ][ 2 ] = B[ 2 ];
 150     mD[ 2 ][ 2 ] = C[ 2 ];
 151     mD[ 3 ][ 2 ] = Q[ 2 ];
 152     ValueType dC = vnl_determinant( mD.GetVnlMatrix( ) );
 153
 154     // Circumcenter
 155     VectorType CS;
 156     CS[ 0 ] = Dx;
 157     CS[ 1 ] = Dy * -_1;
 158     CS[ 2 ] = Dz;
 159     CS /= _2 * dA;
 160
 161     // Circumradius
 162     /*
 163      * WARNING: this is not necessary since sin(phi) and cos(phi)
 164      * depend on R as divisor:
 165      *
 166      * ValueType R = ( Dx * Dx ) + ( Dy * Dy ) + ( Dz * Dz );
 167      * R -= _4 * dA * dC;
 168      * R = ValueType( std::sqrt( double( R ) ) );
 169      * R /= _2 * ValueType( std::fabs( double( dA ) ) );
 170      */
 171
 172     // Phi sin and cos
 173     ValueType sphi = r;
 174     VectorType CCCS = CC - CS;
 175     ValueType cphi = CCCS.GetNorm( );
 176     if( ( QA * N ) < _0 )
 177       sphi *= -_1;
 178     if( ( CCCS * N ) < _0 )
 179       cphi *= -_1;
 180     phi = ValueType( std::atan2( double( sphi ), double( cphi ) ) );
 181   }
 182   else // Tetrahedron IS planar
 183     phi = _0;
 184
 185   return( true );
 186 }
 187
 188 // -------------------------------------------------------------------------
 189 template< typename P, unsigned int K, unsigned int D, typename T >
 190 bool cpm::DataStructures::SimplexMesh< P, K, D, T >::
 191 VerifySimplexOrder( ) const
 192 {
 193   if( this->m_OnextRings.size( ) == 0 )
 194     return( true );
 195
 196   typename Superclass::CntOnextRings::const_iterator rIt =
 197     this->m_OnextRings.begin( );
 198   unsigned int minOrder = ( *rIt )->GetOnextRingSize( );
 199   unsigned int maxOrder = minOrder;
 200   ++rIt;
 201   for( ; rIt != this->m_OnextRings.end( ); ++rIt )
 202   {
 203     unsigned int order = ( *rIt )->GetOnextRingSize( );
 204     minOrder = ( order < minOrder )? order: minOrder;
 205     maxOrder = ( maxOrder < order )? order: minOrder;
 206
 207   } // rof
 208
 209   if( K == 1 )
 210     return(
 211       ( minOrder == 1 || minOrder == 2 ) &&
 212       ( maxOrder == 1 || maxOrder == 2 )
 213       );
 214   else if( K == 2 )
 215     return( ( minOrder == maxOrder ) && ( minOrder == K + 1 ) );
 216   else
 217     return( false );
 218 }
 219
 220 // -------------------------------------------------------------------------
 221 template< typename P, unsigned int K, unsigned int D, typename T >
 222 bool cpm::DataStructures::SimplexMesh< P, K, D, T >::
 223 OperatorDeleteEdge( const PointIdentifier& a, const PointIdentifier& b )
 224 {
 225   TPrimalEdge* e0 = this->FindEdge( a, b );
 226   if( e0 == NULL )
 227     return( false );
 228
 229   // Name edges and faces
 230   TPrimalEdge* e1 = e0->GetOnext( );
 231   TPrimalEdge* e2 = e1->GetOnext( );
 232   TPrimalEdge* e3 = e0->GetLnext( );
 233   TPrimalEdge* e4 = e0->GetRprev( );
 234   TPrimalEdge* e5 = e2->GetSym( )->GetOprev( );
 235   TPrimalEdge* e6 = e4->GetLnext( );
 236
 237   CellIdentifier f0 = e0->GetLeft( );
 238   CellIdentifier f1 = e0->GetRight( );
 239   CellIdentifier f2 = e1->GetLeft( );
 240   CellIdentifier f3 = e3->GetRight( );
 241
 242   // Disconnect edges e0, e2, e4 and origins of e1 and e3
 243   TPrimalEdge::Splice( e2, e0 );
 244   TPrimalEdge::Splice( e2, e1 );
 245   TPrimalEdge::Splice( e5, e2->GetSym( ) );
 246   TPrimalEdge::Splice( e3, e0->GetSym( ) );
 247   TPrimalEdge::Splice( e3, e4 );
 248   TPrimalEdge::Splice( e6, e4->GetSym( ) );
 249
 250   // Connect e1 and e3 origins
 251   TPrimalEdge::Splice( e5, e1 );
 252   TPrimalEdge::Splice( e6, e3 );
 253
 254   // Change e1 and e3 origins
 255   e1->SetOrigin( e5->GetOrigin( ) );
 256   e3->SetOrigin( e6->GetOrigin( ) );
 257
 258   // Assign face id (even if f1 is a "not face")
 259   typename TPrimalEdge::Iterator eIt = e3->BeginLnext( );
 260   while( eIt != e3->EndLnext( ) && *eIt != e1->GetSym( ) )
 261   {
 262     ( *eIt )->SetLeft( f1 );
 263     eIt++;
 264
 265   } // elihw
 266
 267   // Update points
 268   this->m_OnextRings[ e5->GetOrigin( ) ] = e5;
 269   this->m_OnextRings[ e6->GetOrigin( ) ] = e6;
 270
 271   // Update cells entry points
 272   CellAutoPointer c1, c2, c3;
 273   this->GetCell( f1, c1 );
 274   this->GetCell( f2, c2 );
 275   this->GetCell( f3, c3 );
 276
 277   TQuadEdgeCell* ec1 = dynamic_cast< TQuadEdgeCell* >( c1.GetPointer( ) );
 278   TQuadEdgeCell* ec2 = dynamic_cast< TQuadEdgeCell* >( c2.GetPointer( ) );
 279   TQuadEdgeCell* ec3 = dynamic_cast< TQuadEdgeCell* >( c3.GetPointer( ) );
 280
 281   ec1->SetEntryPrimalEdge( e3 );
 282   ec2->SetEntryPrimalEdge( e1 );
 283   ec3->SetEntryPrimalEdge( e3->GetSym( ) );
 284
 285   // Delete face f0
 286   this->_DeleteFace( f0 );
 287
 288   // Delete e0 origin and destination
 289   this->_DeletePoint( e0->GetOrigin( ) );
 290   this->_DeletePoint( e0->GetDestination( ) );
 291
 292   // Delete edges e0, e2 and e4
 293   this->_DeleteEdge( e0 );
 294   this->_DeleteEdge( e2 );
 295   this->_DeleteEdge( e4 );
 296
 297   return( true );
 298 }
 299
 300 // -------------------------------------------------------------------------
 301 template< typename P, unsigned int K, unsigned int D, typename T >
 302 cpm::DataStructures::SimplexMesh< P, K, D, T >::
 303 SimplexMesh( )
 304   : Superclass( )
 305 {
 306 }
 307
 308 // -------------------------------------------------------------------------
 309 template< typename P, unsigned int K, unsigned int D, typename T >
 310 cpm::DataStructures::SimplexMesh< P, K, D, T >::
 311 ~SimplexMesh( )
 312 {
 313 }
 314
 315 // -------------------------------------------------------------------------
 316 template< typename P, unsigned int K, unsigned int D, typename T >
 317 typename cpm::DataStructures::SimplexMesh< P, K, D, T >::
 318 VectorType cpm::DataStructures::SimplexMesh< P, K, D, T >::
 319 _ComputePointNormal( const TPrimalEdge* e ) const
 320 {
 321   VectorType N( ValueType( 0 ) );
 322   if( K == 1 )
 323   {
 324     /* TODO
 325        PointIdentifier iQp0 = e->GetOrigin( );
 326        PointIdentifier iQp1 = e->GetDestination( );
 327        PointIdentifier iQm1 = e->GetOnext( )->GetDestination( );
 328        PointIdentifier iQp2 = e->GetSym( )->GetOnext( )->GetDestination( );
 329        PointIdentifier iQm2 =
 330        e->GetOnext( )->GetSym( )->GetOnext( )->GetDestination( );
 331
 332        VectorType Qp0 = this->GetPoint( iQp0 ).GetVectorFromOrigin( );
 333        VectorType Qp1 = this->GetPoint( iQp1 ).GetVectorFromOrigin( );
 334        VectorType Qm1 = this->GetPoint( iQm1 ).GetVectorFromOrigin( );
 335        VectorType Qp2 = this->GetPoint( iQp2 ).GetVectorFromOrigin( );
 336        VectorType Qm2 = this->GetPoint( iQm2 ).GetVectorFromOrigin( );
 337
 338        // Compute tangent
 339        VectorType T( ValueType( 0 ) );
 340        T = Qp1 - Qm1;
 341        ValueType nT = T.GetNorm( );
 342        if( nT > ValueType( 0 ) )
 343        T /= nT;
 344
 345        // Compute auxiliary vector
 346        VectorType R( ValueType( 0 ) );
 347        R = itk::CrossProduct( T, Qm1 - Qm2 );
 348        R = itk::CrossProduct( R, Qp2 - Qp1 );
 349        ValueType nR = R.GetNorm( );
 350        if( nR > ValueType( 0 ) )
 351        R /= nR;
 352     */
 353     // TODO : psi??? see documentation
 354   }
 355   else if( K == 2 )
 356   {
 357     PointIdentifier iA = e->GetDestination( );
 358     PointIdentifier iB = e->GetOnext( )->GetDestination( );
 359     PointIdentifier iC = e->GetOnext( )->GetOnext( )->GetDestination( );
 360
 361     VectorType A = this->GetPoint( iA ).GetVectorFromOrigin( );
 362     VectorType B = this->GetPoint( iB ).GetVectorFromOrigin( );
 363     VectorType C = this->GetPoint( iC ).GetVectorFromOrigin( );
 364
 365     N  = itk::CrossProduct( B - A, C - A );
 366     ValueType nN = N.GetNorm( );
 367     if( nN > ValueType( 0 ) )
 368       N /= nN;
 369
 370   } // fi
 371   return( N );
 372 }
 373
 374 #endif // __CPM__DATASTRUCTURES__SIMQLEXMESH__HXX__
 375
 376 // eof - $RCSfile$