[Insight-users] Help with GeodesicActiveContourImageFilter --> Modification to example code does not seem to work

[Vimal Singh]

Hi all,

     I'm a new User to ITK and have been learning it for past few days. 
I managed to run the example codes and observe the same results 
available in software guide document.

Currently I'm trying to develop my GeodesicActiveContourImageFilter, 
which is similar to the one provided in the examples/segmentation folder 
but does not use FastMarchingImageFilter for generation of initial Level 
Set. Instead relies on user input for this initial level-set.
I've modified the example code accordingly and observed that the 
DistanceMap obtained using DanielssonDistanceMapImageFilter matches 
exactly the distance-map of initial level set generated by Example code. 
However my final segmentation result is not the same, despite the 
parameters for all other filters in two codes have been kept the same. 
My complete code is  below.

Can anyone give me hints as to what I might be missing. I've removed the 
FastMarchingImageFilter all together. Any kind off help would be highly 
appreciated !!! Or cite me an easier way to debug things ?


#include "itkImage.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
//#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkDanielssonDistanceMapImageFilter.h"  //Include for distance map

int main( int argc, char *argv[] )
{
   if( argc < 10 )
     {
     std::cerr << "Missing Parameters " << std::endl;
     std::cerr << "Usage: " << argv[0];
     std::cerr << " inputImage  initialLevelSet outputImage";
     std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
     std::cerr << " CurvatureScaling AdvectionTerm PropagationScaling" 
<< std::endl;
     return 1;
     }

   std::cout << "Debug Point 1" << std::endl;

   typedef   float           InternalPixelType;
   const     unsigned int    Dimension = 2;
   typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
   typedef unsigned char                            OutputPixelType;
   typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
   typedef itk::BinaryThresholdImageFilter<
                         InternalImageType,
                         OutputImageType >       ThresholdingFilterType;

   ThresholdingFilterType::Pointer thresholder = 
ThresholdingFilterType::New();

   thresholder->SetLowerThreshold( -1000.0 );
   thresholder->SetUpperThreshold(     0.0 );

   thresholder->SetOutsideValue(  0  );
   thresholder->SetInsideValue(  255 );

   typedef itk::DanielssonDistanceMapImageFilter<
                 InternalImageType, InternalImageType > 
DistanceMapFilterType;
   DistanceMapFilterType::Pointer distanceMap = 
DistanceMapFilterType::New();   //Instantiating the DistanceMap filter


   std::cout << "Debug Point 3" << std::endl;
   typedef  itk::ImageFileReader< InternalImageType > ReaderType;
   typedef  itk::ImageFileWriter<  OutputImageType > WriterType;

   ReaderType::Pointer reader = ReaderType::New();  //reads image
   ReaderType::Pointer reader2 = ReaderType::New(); //reads thresholded 
initial level-set   //reading the inital binary level-set image
   WriterType::Pointer writer = WriterType::New();

   reader->SetFileName( argv[1] );
   reader2->SetFileName( argv[2] );
   writer->SetFileName( argv[3] );
   typedef itk::RescaleIntensityImageFilter<
                                InternalImageType,
                                OutputImageType >   CastFilterType;
   typedef   itk::CurvatureAnisotropicDiffusionImageFilter<
                                InternalImageType,
                                InternalImageType >  SmoothingFilterType;

   SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
   typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<
                                InternalImageType,
                                InternalImageType >  GradientFilterType;
   typedef   itk::SigmoidImageFilter<
                                InternalImageType,
                                InternalImageType >  SigmoidFilterType;

   GradientFilterType::Pointer  gradientMagnitude = 
GradientFilterType::New();

   SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
   sigmoid->SetOutputMinimum(  0.0  );
   sigmoid->SetOutputMaximum(  1.0  );

   typedef  itk::GeodesicActiveContourLevelSetImageFilter< 
InternalImageType,
                 InternalImageType >    GeodesicActiveContourFilterType;
   GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
                                      
GeodesicActiveContourFilterType::New();

   const double propagationScaling = atof( argv[9] );
   const double curvatureScaling = atof( argv[7] );
   const double advectionScaling = atof( argv[8] );
   geodesicActiveContour->SetPropagationScaling( propagationScaling );
   geodesicActiveContour->SetCurvatureScaling( curvatureScaling );
   geodesicActiveContour->SetAdvectionScaling( advectionScaling );
   std::cout << "Debug Point 7" << std::endl;
   geodesicActiveContour->SetMaximumRMSError( 0.02 );
   geodesicActiveContour->SetNumberOfIterations(800);

   smoothing->SetInput( reader->GetOutput() );
   gradientMagnitude->SetInput( smoothing->GetOutput() );
   sigmoid->SetInput( gradientMagnitude->GetOutput() );

   distanceMap->SetInput( reader2->GetOutput() );
   geodesicActiveContour->SetInput(  distanceMap->GetOutput() );        
//our distance Mapped initial level-set
   geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );

   thresholder->SetInput( geodesicActiveContour->GetOutput() );
   writer->SetInput( thresholder->GetOutput() );

   smoothing->SetTimeStep( 0.125 );
   smoothing->SetNumberOfIterations(  5 );
   smoothing->SetConductanceParameter( 9.0 );
   std::cout << "Debug Point 9" << std::endl;
   const double sigma = atof( argv[4] );
   gradientMagnitude->SetSigma(  sigma  );
   const double alpha =  atof( argv[5] );
   const double beta  =  atof( argv[6] );
   sigmoid->SetAlpha( alpha );
   sigmoid->SetBeta(  beta  );
   std::cout << "Debug Point 9.1" << std::endl;

   CastFilterType::Pointer caster1 = CastFilterType::New();
   CastFilterType::Pointer caster2 = CastFilterType::New();
   CastFilterType::Pointer caster3 = CastFilterType::New();
   CastFilterType::Pointer caster4 = CastFilterType::New();

   WriterType::Pointer writer1 = WriterType::New();
   WriterType::Pointer writer2 = WriterType::New();
   WriterType::Pointer writer3 = WriterType::New();
   WriterType::Pointer writer4 = WriterType::New();
   std::cout << "Debug Point 9.15" << std::endl;

   caster1->SetInput( smoothing->GetOutput() );
   writer1->SetInput( caster1->GetOutput() );
   writer1->SetFileName("myGeodesicActiveContourImageFilterOutput1.png");
   caster1->SetOutputMinimum(   0 );
   caster1->SetOutputMaximum( 255 );
   writer1->Update();
   std::cout << "Debug Point 9.2" << std::endl;

   caster2->SetInput( gradientMagnitude->GetOutput() );
   writer2->SetInput( caster2->GetOutput() );
   writer2->SetFileName("myGeodesicActiveContourImageFilterOutput2.png");
   caster2->SetOutputMinimum(   0 );
   caster2->SetOutputMaximum( 255 );
   writer2->Update();
   std::cout << "Debug Point 9.3" << std::endl;

   caster3->SetInput( sigmoid->GetOutput() );
   writer3->SetInput( caster3->GetOutput() );
   writer3->SetFileName("myGeodesicActiveContourImageFilterOutput3.png");
   caster3->SetOutputMinimum(   0 );
   caster3->SetOutputMaximum( 255 );
   writer3->Update();

   //caster4->SetInput( distanceMap->GetOutput() );
   caster4->SetInput( geodesicActiveContour->GetInput());
   writer4->SetInput( caster4->GetOutput() );
   writer4->SetFileName("myGeodesicActiveContourImageFilterOutput4.png");
   caster4->SetOutputMinimum(   0 );
   caster4->SetOutputMaximum( 255 );
   //writer4->Update();

   std::cout << "Debug Point 10" << std::endl;

   try
     {
     writer->Update();
     }
   catch( itk::ExceptionObject & excep )
     {
     std::cerr << "Exception caught !" << std::endl;
     std::cerr << excep << std::endl;
     }


   // Print out some useful information
   std::cout << std::endl;
   std::cout << "Max. no. iterations: " << 
geodesicActiveContour->GetNumberOfIterations() << std::endl;
   std::cout << "Max. RMS error: " << 
geodesicActiveContour->GetMaximumRMSError() << std::endl;
   std::cout << std::endl;
   std::cout << "No. elpased iterations: " << 
geodesicActiveContour->GetElapsedIterations() << std::endl;
   std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() 
<< std::endl;

   writer4->Update();
   return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110201/b0d44a59/attachment.htm>