]> Creatis software - FrontAlgorithms.git/blobdiff - lib/fpa/Common/SliceBySliceRandomWalker.hxx
...
[FrontAlgorithms.git] / lib / fpa / Common / SliceBySliceRandomWalker.hxx
index d1fbd9cfdf495e39a5036db6fdbc1e2ba96c9c12..3b86d023b9b56f65de65e71121648c6fafce3025 100644 (file)
@@ -6,10 +6,12 @@
 #define __fpa__Common__SliceBySliceRandomWalker__hxx__
 
 #include <vector>
+#include <itkBinaryThresholdImageFilter.h>
 #include <itkExtractImageFilter.h>
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionIterator.h>
 #include <itkJoinSeriesImageFilter.h>
+#include <itkMaximumImageFilter.h>
 #include <itkMinimumMaximumImageCalculator.h>
 #include <fpa/Common/RandomWalker.h>
 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
 template< class _TImage, class _TLabels, class _TScalarImage >
 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 SliceBySliceRandomWalker( )
-  : Superclass( )
+  : Superclass( ),
+    m_Epsilon( TScalar( 1e-5 ) ),
+    m_Beta( TScalar( 1 ) ),
+    m_VesselnessThreshold( TScalar( 5 ) )
 {
   ivqITKInputConfigureMacro( InputLabels, TLabels );
   ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
@@ -47,21 +52,25 @@ GenerateData( )
   const TScalarImage* vesselness = this->GetInputVesselness( );
   this->AllocateOutputs( );
 
+  // Get vesselness range
+  typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+  typename _TMinMax::Pointer vMinMax = _TMinMax::New( );
+  vMinMax->SetImage( vesselness );
+  vMinMax->Compute( );
+  TScalar vMax = vMinMax->GetMaximum( );
+
   // Some values
   typename TImage::RegionType region = input->GetRequestedRegion( );
-  typename TImage::SizeType regionIndex = region.GetIndex( );
+  typename TImage::IndexType regionIndex = region.GetIndex( );
   typename TImage::SizeType sliceSize = region.GetSize( );
   int numSlices = sliceSize[ 2 ];
   sliceSize[ 2 ] = 0;
-  this->m_VesselnessValue =
-    this->m_VesselnessThreshold * this->m_VesselnessMax / TScalar( 100 );
+  this->m_VesselnessValue  = TScalar( 0.01 );
+  this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax;
 
   // Composite image
-  typename TScalarImage::Pointer composite = TScalarImage::New( );
-  composite->CopyInformation( vesselness );
-  composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
-  composite->Allocate( );
-  this->_Composite( composite, labels, vesselness );
+  typename TScalarImage::Pointer composite;
+  this->_Composite( composite, labels, vesselness, vMax );
 
   // Extract slices
   std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
@@ -80,7 +89,7 @@ GenerateData( )
     typename _TSliceScalarImage::Pointer composite_slice;
     this->_Slice( input_slice, input, sliceRegion );
     this->_Slice( labels_slice, labels, sliceRegion );
-    this->_Slice( composite_slice, composite, sliceRegion );
+    this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion );
 
     // Update vectors with each image
     data3D[ i ] = input_slice;
@@ -90,14 +99,16 @@ GenerateData( )
   } // rof
 
   // Random walker slice-by-slice
-  this->_GoDown( binaryTree, data3D, fusion, numSlices );
-  this->_GoUp( binaryTree, data3D, fusion, numSlices );
+  this->_RandomWalker( binaryTree, data3D, fusion, true );
+  this->_RandomWalker( binaryTree, data3D, fusion, false );
 
   // Join results
-  typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
+  typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin;
   typename _TJoin::Pointer join = _TJoin::New( );
   for( int i = 0; i < numSlices; ++i )
     join->SetInput( i, binaryTree[ i ] );
+  join->SetOrigin( input->GetOrigin( )[ 2 ] );
+  join->SetSpacing( input->GetSpacing( )[ 2 ] );
   join->Update( );
   this->GetOutput( )->Graft( join->GetOutput( ) );
 }
@@ -109,19 +120,15 @@ fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 _Composite(
   typename TScalarImage::Pointer& composite,
   const TLabels* labels,
-  const TScalarImage* vesselness
+  const TScalarImage* vesselness,
+  const TScalar& maxVess
   )
 {
-  // Min-Max vesselness values
-  typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
-  typename _TMinMax::Pointer minMax = _TMinMax::New( );
-  minMax->SetImage( vesselness );
-  minMax->Compute( );
-  TScalar maxVess = minMax->GetMaximum( );
+  // Some values
   typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
 
   // Composite image
-  typename TScalarImage::Pointer composite = TScalarImage::New( );
+  composite = TScalarImage::New( );
   composite->CopyInformation( vesselness );
   composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
   composite->Allocate( );
@@ -146,7 +153,7 @@ _Slice(
   typename _TInput::RegionType region
   )
 {
-  typedef typename _TSlicePTr::ObjectType _TSlice;
+  typedef typename _TSlicePtr::ObjectType _TSlice;
   typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
 
   typename _TExtract::Pointer extract = _TExtract::New( );
@@ -163,11 +170,11 @@ template< class _TImage, class _TLabels, class _TScalarImage >
 template< class _TBinaryTree, class _TData3D, class _TFusion >
 void
 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
-_GoDown(
+_RandomWalker(
   _TBinaryTree& binaryTree,
   const _TData3D& data3D,
   const _TFusion& fusion,
-  int numSlices
+  bool down
   )
 {
   typedef typename _TBinaryTree::value_type     _TSliceBinaryPtr;
@@ -178,20 +185,22 @@ _GoDown(
   typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
   typedef typename _TSliceFusion::PixelType     _TSliceScalar;
 
-  int z = -1;
-  while( z < numSlices - 2 )
+  int numSlices = binaryTree.size( );
+  int z = ( down )? -1: numSlices;
+  int o = ( down )? 1: -1;
+  while( ( down && z < numSlices - 2 ) || ( !down && z > 1 ) )
   {
-    z++;
+    z += o;
     _TSliceBinaryPtr tmp = binaryTree[ z ];
-    _TSliceBinaryPtr next = binaryTree[ z  + 1 ];
-    _TSliceData3DPtr data = data3D[ z + 1 ];
-    _TSliceFusionPtr vess = fusion[ z + 1 ];
+    _TSliceBinaryPtr next = binaryTree[ z + o ];
+    _TSliceData3DPtr data = data3D[ z + o ];
+    _TSliceFusionPtr vess = fusion[ z + o ];
     typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
 
     // Create seeds and background
     _TSliceBinaryPtr seeds = _TSliceBinary::New( );
     seeds->CopyInformation( tmp );
-    seeds->SetRequestedRegion( region );
+    seeds->SetBufferedRegion( region );
     seeds->Allocate( );
     seeds->FillBuffer( 0 );
 
@@ -232,9 +241,9 @@ _GoDown(
 
     // Random walker function
     typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
-    typename TFunction::Pointer edge = TFunction::New( );
+    typename _TFunction::Pointer edge = _TFunction::New( );
     edge->SetBeta( this->m_Beta );
-    edge->SetEpsilon( this->m_Eps );
+    edge->SetEpsilon( this->m_Epsilon );
     edge->TreatAsWeightOff( );
 
     // Real random walker
@@ -244,32 +253,25 @@ _GoDown(
     rw->SetInputLabels( seeds );
     rw->SetEdgeFunction( edge );
 
-    // Keep maximum
-    /* TODO
-       typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
-       MaxFilterType::Pointer maxFilter = MaxFilterType::New();
-       maxFilter->SetInput(0, img2);
-       maxFilter->SetInput(1, next);
-       maxFilter->Update();
-
-       SliceType::Pointer sth = maxFilter->GetOutput();
-       ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());
-       delete[] probs;
-    */
-  }
-}
+    // Extract label
+    typedef itk::BinaryThresholdImageFilter< _TSliceBinary, _TSliceBinary > _TExtract;
+    typename _TExtract::Pointer extract = _TExtract::New( );
+    extract->SetInput( rw->GetOutput( ) );
+    extract->SetInsideValue( 1 );
+    extract->SetOutsideValue( 0 );
+    extract->SetLowerThreshold( 1 );
+    extract->SetUpperThreshold( 1 );
 
-// -------------------------------------------------------------------------
-template< class _TImage, class _TLabels, class _TScalarImage >
-template< class _TBinaryTree, class _TData3D, class _TFusion >
-void
-_GoUp(
-  _TBinaryTree& binaryTree,
-  const _TData3D& data3D,
-  const _TFusion& fusion,
-  int numSlices
-  )
-{
+    // Keep maximum
+    typedef itk::MaximumImageFilter< _TSliceBinary > _TMax;
+    typename _TMax::Pointer max = _TMax::New();
+    max->SetInput( 0, extract->GetOutput( ) );
+    max->SetInput( 1, next );
+    max->Update( );
+    binaryTree[ z + o ] = max->GetOutput( );
+    binaryTree[ z + o ]->DisconnectPipeline( );
+
+  } // elihw
 }
 
 #endif // __fpa__Common__SliceBySliceRandomWalker__hxx__