lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx

   1 #ifndef __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   2 #define __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   3
   4 #include <cmath>
   5 #include <itkImageRegionIteratorWithIndex.h>
   6
   7 // -------------------------------------------------------------------------
   8 template< class _TImage, class _TMST >
   9 const typename
  10 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  11 TImage*
  12 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  13 GetCostsImage( )
  14 {
  15   return(
  16     dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 1 ) )
  17     );
  18 }
  19
  20 // -------------------------------------------------------------------------
  21 template< class _TImage, class _TMST >
  22 const typename
  23 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  24 TImage*
  25 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  26 GetDistanceMap( )
  27 {
  28   return(
  29     dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 2 ) )
  30     );
  31 }
  32
  33 // -------------------------------------------------------------------------
  34 template< class _TImage, class _TMST >
  35 typename
  36 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  37 TSkeleton*
  38 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  39 GetSkeleton( )
  40 {
  41   return(
  42     dynamic_cast< TSkeleton* >(
  43       this->itk::ProcessObject::GetOutput(
  44         this->GetNumberOfRequiredOutputs( ) - 1
  45         )
  46       )
  47     );
  48 }
  49
  50 // -------------------------------------------------------------------------
  51 template< class _TImage, class _TMST >
  52 void
  53 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  54 SetCostsImage( TImage* image )
  55 {
  56   this->itk::ProcessObject::SetNthInput( 1, const_cast< _TImage* >( image ) );
  57 }
  58
  59 // -------------------------------------------------------------------------
  60 template< class _TImage, class _TMST >
  61 void
  62 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  63 SetDistanceMap( TImage* image )
  64 {
  65   this->itk::ProcessObject::SetNthInput( 2, const_cast< _TImage* >( image ) );
  66 }
  67
  68 // -------------------------------------------------------------------------
  69 template< class _TImage, class _TMST >
  70 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  71 ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  72   : Superclass( ),
  73     m_SquaredDistanceMap( false )
  74 {
  75   this->SetNumberOfRequiredInputs( 3 );
  76   unsigned int nOuts = this->GetNumberOfRequiredOutputs( );
  77   this->SetNumberOfRequiredOutputs( nOuts + 1 );
  78   typename TSkeleton::Pointer sk = TSkeleton::New( );
  79   this->itk::ProcessObject::SetNthOutput( nOuts, sk.GetPointer( ) );
  80 }
  81
  82 // -------------------------------------------------------------------------
  83 template< class _TImage, class _TMST >
  84 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  85 ~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  86 {
  87 }
  88
  89 // -------------------------------------------------------------------------
  90 template< class _TImage, class _TMST >
  91 void
  92 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  93 GenerateData( )
  94 {
  95   // Create auxiliary objects
  96   auto image = this->GetCostsImage( );
  97   this->m_MarkImage = TMarkImage::New( );
  98   this->m_MarkImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
  99   this->m_MarkImage->SetRequestedRegion( image->GetRequestedRegion( ) );
 100   this->m_MarkImage->SetBufferedRegion( image->GetBufferedRegion( ) );
 101   this->m_MarkImage->SetDirection( image->GetDirection( ) );
 102   this->m_MarkImage->SetOrigin( image->GetOrigin( ) );
 103   this->m_MarkImage->SetSpacing( image->GetSpacing( ) );
 104   this->m_MarkImage->Allocate( );
 105   this->m_MarkImage->FillBuffer( 0 );
 106   this->m_SkeletonImage = TMarkImage::New( );
 107   this->m_SkeletonImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
 108   this->m_SkeletonImage->SetRequestedRegion( image->GetRequestedRegion( ) );
 109   this->m_SkeletonImage->SetBufferedRegion( image->GetBufferedRegion( ) );
 110   this->m_SkeletonImage->SetDirection( image->GetDirection( ) );
 111   this->m_SkeletonImage->SetOrigin( image->GetOrigin( ) );
 112   this->m_SkeletonImage->SetSpacing( image->GetSpacing( ) );
 113   this->m_SkeletonImage->Allocate( );
 114   this->m_SkeletonImage->FillBuffer( 0 );
 115
 116   // Real execution
 117   this->Superclass::GenerateData( );
 118
 119   // Build skeleton
 120   auto sk = this->GetSkeleton( );
 121   sk->SetMinimumSpanningTree( this->GetMinimumSpanningTree( ) );
 122   auto bIt = this->m_Branches.begin( );
 123   for( ; bIt != this->m_Branches.end( ); ++bIt )
 124     sk->AddBranch( bIt->first, bIt->second );
 125 }
 126
 127 // -------------------------------------------------------------------------
 128 template< class _TImage, class _TMST >
 129 void
 130 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 131 _MarkSkeleton( const TVertex& v, const unsigned long& l )
 132 {
 133   if( this->m_SkeletonImage->GetRequestedRegion( ).IsInside( v ) )
 134     this->m_SkeletonImage->SetPixel( v, l );
 135 }
 136
 137 // -------------------------------------------------------------------------
 138 template< class _TImage, class _TMST >
 139 void
 140 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 141 _MarkSphere(
 142   const TVertex& v,
 143   const double& r,
 144   const unsigned long& l
 145   )
 146 {
 147   double rr = r * double( 1.5 );
 148
 149   // Get marking region
 150   auto rreg = this->m_MarkImage->GetRequestedRegion( );
 151   auto spac = this->m_MarkImage->GetSpacing( );
 152   TVertex ireg = rreg.GetIndex( );
 153   TVertex jreg = ireg + rreg.GetSize( );
 154   TVertex idx, jdx;
 155   typename _TImage::SizeType size;
 156   bool all = false;
 157   for( unsigned int d = 0; d < _TImage::ImageDimension; ++d )
 158   {
 159     unsigned long s = std::ceil( rr / double( spac[ d ] ) );
 160     if( s < 3 )
 161     {
 162       s = 3;
 163       all = true;
 164
 165     } // fi
 166     s += 3;
 167     idx[ d ] = v[ d ] - s;
 168     jdx[ d ] = v[ d ] + s;
 169
 170     jreg[ d ]--;
 171     if( idx[ d ] < ireg[ d ] ) idx[ d ] = ireg[ d ];
 172     if( idx[ d ] > jreg[ d ] ) idx[ d ] = jreg[ d ];
 173     if( jdx[ d ] < ireg[ d ] ) jdx[ d ] = ireg[ d ];
 174     if( jdx[ d ] > jreg[ d ] ) jdx[ d ] = jreg[ d ];
 175     size[ d ] = jdx[ d ] - idx[ d ] + 1;
 176
 177   } // rof
 178   typename _TImage::RegionType region;
 179   region.SetIndex( idx );
 180   region.SetSize( size );
 181
 182   // Mark region
 183   typename _TImage::PointType c;
 184   this->m_MarkImage->TransformIndexToPhysicalPoint( v, c );
 185   itk::ImageRegionIteratorWithIndex< TMarkImage >
 186     spIt( this->m_MarkImage, region );
 187   for( spIt.GoToBegin( ); !spIt.IsAtEnd( ); ++spIt )
 188   {
 189     typename _TImage::PointType pnt;
 190     this->m_MarkImage->TransformIndexToPhysicalPoint( spIt.GetIndex( ), pnt );
 191     if( double( pnt.EuclideanDistanceTo( c ) ) <= rr || all )
 192       if( spIt.Get( ) == 0 )
 193         spIt.Set( l );
 194
 195   } // rof
 196 }
 197
 198 // -------------------------------------------------------------------------
 199 template< class _TImage, class _TMST >
 200 unsigned long
 201 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 202 _Mark( const TVertex& v )
 203 {
 204   return( this->m_MarkImage->GetPixel( v ) );
 205 }
 206
 207 // -------------------------------------------------------------------------
 208 template< class _TImage, class _TMST >
 209 unsigned long fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 210 _SkeletonMark( const TVertex& v )
 211 {
 212   return( this->m_SkeletonImage->GetPixel( v ) );
 213 }
 214
 215 // -------------------------------------------------------------------------
 216 template< class _TImage, class _TMST >
 217 double
 218 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 219 _Radius( const TVertex& v )
 220 {
 221   double d = double( this->GetDistanceMap( )->GetPixel( v ) );
 222   double s = std::fabs( d );
 223   if( s > double( 0 ) )
 224     s /= d;
 225   if( this->m_SquaredDistanceMap )
 226     return( s * std::sqrt( std::fabs( d ) ) );
 227   else
 228     return( d );
 229 }
 230
 231 #endif // __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
 232
 233 // eof - $RCSfile$