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