lib/fpa/Common/SliceBySliceRandomWalker.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5 #ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__
   6 #define __fpa__Common__SliceBySliceRandomWalker__hxx__
   7
   8 #include <vector>
   9 #include <itkBinaryThresholdImageFilter.h>
  10 #include <itkExtractImageFilter.h>
  11 #include <itkImageRegionConstIterator.h>
  12 #include <itkImageRegionIterator.h>
  13 #include <itkJoinSeriesImageFilter.h>
  14 #include <itkMaximumImageFilter.h>
  15 #include <itkMinimumMaximumImageCalculator.h>
  16 #include <fpa/Common/RandomWalker.h>
  17 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
  18
  19 // -------------------------------------------------------------------------
  20 template< class _TImage, class _TLabels, class _TScalarImage >
  21 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  22 SliceBySliceRandomWalker( )
  23   : Superclass( ),
  24     m_Epsilon( TScalar( 1e-5 ) ),
  25     m_Beta( TScalar( 1 ) ),
  26     m_VesselnessThreshold( TScalar( 5 ) )
  27 {
  28   ivqITKInputConfigureMacro( InputLabels, TLabels );
  29   ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
  30 }
  31
  32 // -------------------------------------------------------------------------
  33 template< class _TImage, class _TLabels, class _TScalarImage >
  34 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  35 ~SliceBySliceRandomWalker( )
  36 {
  37 }
  38
  39 // -------------------------------------------------------------------------
  40 template< class _TImage, class _TLabels, class _TScalarImage >
  41 void
  42 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  43 GenerateData( )
  44 {
  45   typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
  46   typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
  47   typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;
  48
  49   // Preare I/O objects
  50   const TImage* input = this->GetInput( );
  51   const TLabels* labels = this->GetInputLabels( );
  52   const TScalarImage* vesselness = this->GetInputVesselness( );
  53   this->AllocateOutputs( );
  54
  55   // Get vesselness range
  56   typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
  57   typename _TMinMax::Pointer vMinMax = _TMinMax::New( );
  58   vMinMax->SetImage( vesselness );
  59   vMinMax->Compute( );
  60   TScalar vMax = vMinMax->GetMaximum( );
  61
  62   // Some values
  63   typename TImage::RegionType region = input->GetRequestedRegion( );
  64   typename TImage::IndexType regionIndex = region.GetIndex( );
  65   typename TImage::SizeType sliceSize = region.GetSize( );
  66   int numSlices = sliceSize[ 2 ];
  67   sliceSize[ 2 ] = 0;
  68   this->m_VesselnessValue  = TScalar( 0.01 );
  69   this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax;
  70
  71   // Composite image
  72   typename TScalarImage::Pointer composite;
  73   this->_Composite( composite, labels, vesselness, vMax );
  74
  75   // Extract slices
  76   std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
  77   std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
  78   std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
  79   for( int i = 0; i < numSlices; ++i )
  80   {
  81     // Prepare extraction region
  82     typename TImage::IndexType sliceIndex = regionIndex;
  83     sliceIndex[ 2 ] += i;
  84     typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );
  85
  86     // Extract regions
  87     typename _TSliceImage::Pointer input_slice;
  88     typename _TSliceLabels::Pointer labels_slice;
  89     typename _TSliceScalarImage::Pointer composite_slice;
  90     this->_Slice( input_slice, input, sliceRegion );
  91     this->_Slice( labels_slice, labels, sliceRegion );
  92     this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion );
  93
  94     // Update vectors with each image
  95     data3D[ i ] = input_slice;
  96     binaryTree[ i ] = labels_slice;
  97     fusion[ i ] = composite_slice;
  98
  99   } // rof
 100
 101   // Random walker slice-by-slice
 102   this->_RandomWalker( binaryTree, data3D, fusion, true );
 103   this->_RandomWalker( binaryTree, data3D, fusion, false );
 104
 105   // Join results
 106   typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin;
 107   typename _TJoin::Pointer join = _TJoin::New( );
 108   for( int i = 0; i < numSlices; ++i )
 109     join->SetInput( i, binaryTree[ i ] );
 110   join->SetOrigin( input->GetOrigin( )[ 2 ] );
 111   join->SetSpacing( input->GetSpacing( )[ 2 ] );
 112   join->Update( );
 113   this->GetOutput( )->Graft( join->GetOutput( ) );
 114 }
 115
 116 // -------------------------------------------------------------------------
 117 template< class _TImage, class _TLabels, class _TScalarImage >
 118 void
 119 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 120 _Composite(
 121   typename TScalarImage::Pointer& composite,
 122   const TLabels* labels,
 123   const TScalarImage* vesselness,
 124   const TScalar& maxVess
 125   )
 126 {
 127   // Some values
 128   typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
 129
 130   // Composite image
 131   composite = TScalarImage::New( );
 132   composite->CopyInformation( vesselness );
 133   composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
 134   composite->Allocate( );
 135   itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
 136   itk::ImageRegionConstIterator< TLabels > l( labels, region );
 137   itk::ImageRegionIterator< TScalarImage > c( composite, region );
 138   v.GoToBegin( );
 139   l.GoToBegin( );
 140   c.GoToBegin( );
 141   for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
 142     c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
 143 }
 144
 145 // -------------------------------------------------------------------------
 146 template< class _TImage, class _TLabels, class _TScalarImage >
 147 template< class _TSlicePtr, class _TInput >
 148 void
 149 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 150 _Slice(
 151   _TSlicePtr& slice,
 152   const _TInput* input,
 153   typename _TInput::RegionType region
 154   )
 155 {
 156   typedef typename _TSlicePtr::ObjectType _TSlice;
 157   typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
 158
 159   typename _TExtract::Pointer extract = _TExtract::New( );
 160   extract->SetInput( input );
 161   extract->SetExtractionRegion( region );
 162   extract->SetDirectionCollapseToIdentity( );
 163   extract->Update( );
 164   slice = extract->GetOutput( );
 165   slice->DisconnectPipeline( );
 166 }
 167
 168 // -------------------------------------------------------------------------
 169 template< class _TImage, class _TLabels, class _TScalarImage >
 170 template< class _TBinaryTree, class _TData3D, class _TFusion >
 171 void
 172 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 173 _RandomWalker(
 174   _TBinaryTree& binaryTree,
 175   const _TData3D& data3D,
 176   const _TFusion& fusion,
 177   bool down
 178   )
 179 {
 180   typedef typename _TBinaryTree::value_type     _TSliceBinaryPtr;
 181   typedef typename _TData3D::value_type         _TSliceData3DPtr;
 182   typedef typename _TFusion::value_type         _TSliceFusionPtr;
 183   typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
 184   typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
 185   typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
 186   typedef typename _TSliceFusion::PixelType     _TSliceScalar;
 187
 188   int numSlices = binaryTree.size( );
 189   int z = ( down )? -1: numSlices;
 190   int o = ( down )? 1: -1;
 191   while( ( down && z < numSlices - 2 ) || ( !down && z > 1 ) )
 192   {
 193     z += o;
 194     _TSliceBinaryPtr tmp = binaryTree[ z ];
 195     _TSliceBinaryPtr next = binaryTree[ z + o ];
 196     _TSliceData3DPtr data = data3D[ z + o ];
 197     _TSliceFusionPtr vess = fusion[ z + o ];
 198     typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
 199
 200     // Create seeds and background
 201     _TSliceBinaryPtr seeds = _TSliceBinary::New( );
 202     seeds->CopyInformation( tmp );
 203     seeds->SetBufferedRegion( region );
 204     seeds->Allocate( );
 205     seeds->FillBuffer( 0 );
 206
 207     itk::ImageRegionConstIterator< _TSliceBinary > tIt( tmp, region );
 208     itk::ImageRegionConstIterator< _TSliceFusion > vIt( vess, region );
 209     itk::ImageRegionIterator< _TSliceBinary > sIt( seeds, region );
 210     tIt.GoToBegin( );
 211     vIt.GoToBegin( );
 212     sIt.GoToBegin( );
 213     int len = 0;
 214     int numO = 0;
 215     int numB = 0;
 216     for( ; !tIt.IsAtEnd( ); ++tIt, ++vIt, ++sIt )
 217     {
 218       if( tIt.Get( ) != 0 )
 219       {
 220         len++;
 221         if( vIt.Get( ) > this->m_VesselnessValue )
 222         {
 223           sIt.Set( 1 );
 224           numO++;
 225
 226         } // fi
 227
 228       } // fi
 229
 230       if( vIt.Get( ) < this->m_VesselnessValue )
 231       {
 232         sIt.Set( 2 );
 233         numB++;
 234
 235       } // fi
 236
 237     } // rof
 238
 239     if( len == 0 )
 240       continue;
 241
 242     // Random walker function
 243     typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
 244     typename _TFunction::Pointer edge = _TFunction::New( );
 245     edge->SetBeta( this->m_Beta );
 246     edge->SetEpsilon( this->m_Epsilon );
 247     edge->TreatAsWeightOff( );
 248
 249     // Real random walker
 250     typedef fpa::Common::RandomWalker< _TSliceData3D, _TSliceBinary, _TSliceScalar > _TRandomWalker;
 251     typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
 252     rw->SetInput( data );
 253     rw->SetInputLabels( seeds );
 254     rw->SetEdgeFunction( edge );
 255
 256     // Extract label
 257     typedef itk::BinaryThresholdImageFilter< _TSliceBinary, _TSliceBinary > _TExtract;
 258     typename _TExtract::Pointer extract = _TExtract::New( );
 259     extract->SetInput( rw->GetOutput( ) );
 260     extract->SetInsideValue( 1 );
 261     extract->SetOutsideValue( 0 );
 262     extract->SetLowerThreshold( 1 );
 263     extract->SetUpperThreshold( 1 );
 264
 265     // Keep maximum
 266     typedef itk::MaximumImageFilter< _TSliceBinary > _TMax;
 267     typename _TMax::Pointer max = _TMax::New();
 268     max->SetInput( 0, extract->GetOutput( ) );
 269     max->SetInput( 1, next );
 270     max->Update( );
 271     binaryTree[ z + o ] = max->GetOutput( );
 272     binaryTree[ z + o ]->DisconnectPipeline( );
 273
 274   } // elihw
 275 }
 276
 277 #endif // __fpa__Common__SliceBySliceRandomWalker__hxx__
 278 // eof - $RCSfile$