...

[FrontAlgorithms.git] / appli / CTArteries / algorithms / RandomWalkSegmentation.hxx

diff --git a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx
index 15fb9aff0135b4571635b700a0e1364f7f9843a9..30738e96e92ba2e6d75e0d6fa6aff3d9b4513031 100644 (file)
--- a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx
+++ b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx
@@ -12,15 +12,7 @@
 #include "DijkstraWithMeanAndVariance.h"
 #include "RandomWalkLabelling.h"
 
-/* TODO
-   #include <itkImageRegionConstIterator.h>
-   #include <itkImageRegionIterator.h>
-
-
-   #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
-
-   #include <itkImageFileWriter.h>
-*/
+#include <itkImageFileWriter.h>
 
 // -------------------------------------------------------------------------
 template< class _TInputImage, class _TOutputImage >
@@ -140,9 +132,12 @@ GenerateData( )
 
   // Smooth input
   typename TOutputImage::Pointer smooth_in;
-  this->_Smooth( this->GetInput( ), smooth_in );
+  std::cout << "smooth" << std::endl;
+  this->_Smooth( this->GetInput( ), smooth_in, 2 );
+  this->_Save( smooth_in, "smooth.mhd" );
 
   // Initial segmentation
+  std::cout << "raw" << std::endl;
   typename TOutputImage::Pointer init_seg;
   typename TPath::Pointer init_axis;
   _TScalar init_mean, init_std;
@@ -153,9 +148,12 @@ GenerateData( )
       this->m_Beta,
       init_seg, init_axis, init_mean, init_std
       );
+  std::cout << "Stat: " << init_mean << " +/- " << init_std << std::endl;
   init_std *= _TScalar( this->m_Sigma );
+  this->_Save( init_seg, "raw.mhd" );
 
   // Extract input ROIs
+  std::cout << "ROI" << std::endl;
   typename TOutputImage::Pointer smooth_in_roi, init_seg_roi;
   roi = this->_ROI( smooth_in.GetPointer( ), roi, 10, smooth_in_roi );
   this->_ROI( init_seg.GetPointer( ), roi, 0, init_seg_roi );
@@ -164,6 +162,7 @@ GenerateData( )
   this->_AxisROI( init_axis.GetPointer( ), roi, init_axis_roi );
 
   // Labelling
+  std::cout << "labelling" << std::endl;
   typename _TLabels::Pointer init_labels;
   _TScalar radius = _TScalar( this->m_Radius );
   this->_Label(
@@ -173,17 +172,20 @@ GenerateData( )
     init_mean, init_std, radius,
     init_labels
     );
+  this->_Save( init_labels, "init_labels.mhd" );
 
   // Random walker
+  std::cout << "random walker " << init_std << " " << this->m_Beta << std::endl;
   typename _TLabels::Pointer rw_seg;
   this->_RandomWalker(
     smooth_in_roi.GetPointer( ),
     init_labels.GetPointer( ),
-    init_std / _TScalar( 2 ),
+    this->m_Beta, // init_std / _TScalar( 2 ),
     rw_seg
     );
 
   // ROI outputs
+  std::cout << "axis" << std::endl;
   typename TOutputImage::Pointer out_dist;
   typename TPath::Pointer out_axis;
   this->_DistanceAndAxis(
@@ -193,37 +195,34 @@ GenerateData( )
     );
 
   // Put everything back to requested region
-  /* TODO
-     std::cout << "6" << std::endl;
-     { // begin
-     TOutputImage* output = this->GetOutput( );
-     output->SetBufferedRegion( output->GetRequestedRegion( ) );
-     output->Allocate( );
-     output->FillBuffer( -std::numeric_limits< _TScalar >::max( ) );
-
-     itk::ImageRegionConstIterator< TOutputImage > rIt(
-     output_roi, output_roi->GetRequestedRegion( )
-     );
-     itk::ImageRegionIterator< TOutputImage > oIt( output, roi );
-     rIt.GoToBegin( );
-     oIt.GoToBegin( );
-     for( ; !rIt.IsAtEnd( ); ++rIt, ++oIt )
-     oIt.Set( rIt.Get( ) );
-
-     TPath* output_axis = this->GetOutputAxis( );
-     output_axis->SetReferenceImage( output );
-     for( unsigned long i = 0; i < output_axis_roi->GetSize( ); ++i )
-     {
-     TIndex v = output_axis_roi->GetVertex( i );
-     for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
-     v[ d ] += roi.GetIndex( )[ d ];
-     output_axis->AddVertex( v );
-
-     } // rof
-
-     } // end
-     std::cout << "7" << std::endl;
-  */
+  std::cout << "output" << std::endl;
+  { // begin
+    TOutputImage* output = this->GetOutput( );
+    output->SetBufferedRegion( output->GetRequestedRegion( ) );
+    output->Allocate( );
+    output->FillBuffer( -std::numeric_limits< _TScalar >::max( ) );
+
+    itk::ImageRegionConstIterator< TOutputImage > rIt(
+      out_dist, out_dist->GetRequestedRegion( )
+      );
+    itk::ImageRegionIterator< TOutputImage > oIt( output, roi );
+    rIt.GoToBegin( );
+    oIt.GoToBegin( );
+    for( ; !rIt.IsAtEnd( ); ++rIt, ++oIt )
+      oIt.Set( rIt.Get( ) );
+
+    TPath* output_axis = this->GetOutputAxis( );
+    output_axis->SetReferenceImage( output );
+    for( unsigned long i = 0; i < out_axis->GetSize( ); ++i )
+    {
+      TIndex v = out_axis->GetVertex( i );
+      for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
+        v[ d ] += roi.GetIndex( )[ d ];
+      output_axis->AddVertex( v );
+
+    } // rof
+
+  } // end
 }
 
 // -------------------------------------------------------------------------
@@ -243,7 +242,7 @@ _SynchSeed( const _TIn* in, _TSeed& seed )
 template< class _TInputImage, class _TOutputImage >
 template< class _TIn, class _TOutPtr >
 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
-_Smooth( const _TIn* in, _TOutPtr& out )
+_Smooth( const _TIn* in, _TOutPtr& out, double s )
 {
   typedef typename _TOutPtr::ObjectType _TOut;
   typedef itk::SmoothingRecursiveGaussianImageFilter< _TIn, _TOut > _TSmooth;
@@ -251,7 +250,7 @@ _Smooth( const _TIn* in, _TOutPtr& out )
   typename _TSmooth::Pointer smooth = _TSmooth::New( );
   smooth->SetInput( in );
   smooth->SetNormalizeAcrossScale( true );
-  smooth->SetSigmaArray( in->GetSpacing( ) * double( 2 ) );
+  smooth->SetSigmaArray( in->GetSpacing( ) * s );
   smooth->Update( );
   out = smooth->GetOutput( );
   out->DisconnectPipeline( );
@@ -380,6 +379,10 @@ _Label(
   label->Update( );
   out = label->GetOutputLabels( );
   out->DisconnectPipeline( );
+
+  std::cout << label->GetLowerThreshold( ) << std::endl;
+  std::cout << label->GetUpperThreshold( ) << std::endl;
+
 }
 
 // -------------------------------------------------------------------------
@@ -447,7 +450,20 @@ _DistanceAndAxis(
   extract->Update( );
   out_axis = TPath::New( );
   out_axis->Graft( extract->GetOutput( ) );
-  this->_Smooth( extract->GetCenterness( )->GetOutput( ), out_dist );
+  this->_Smooth( extract->GetCenterness( )->GetOutput( ), out_dist, 1 );
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+template< class _TInPtr >
+void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
+_Save( const _TInPtr& in, const std::string& fname )
+{
+  typedef itk::ImageFileWriter< typename _TInPtr::ObjectType > _TWriter;
+  typename _TWriter::Pointer w = _TWriter::New( );
+  w->SetInput( in );
+  w->SetFileName( fname );
+  w->Update( );
 }
 
 #endif // __RandomWalkSegmentation__hxx__