[Insight-users] Help with GeodesicActiveContourImageFilter --> Modification to example code does not seem to work
Kishore Mosaliganti
kishoreraom at gmail.com
Tue Feb 1 16:22:19 EST 2011
Hi Vimal,
The levelset function is supposed to be initialized in such a way that
there is a negative and positive part in the range of the function.
The fast marching filter allows one to specify negative seeds and
therefore it generates a 0 levelset which is a circle.
The output of the Danielsson Distance map is a fully positive function.
To make it work, subtract some value from its output before setting it
as input to the levelset filter.
Kishore
On Tue, Feb 1, 2011 at 4:15 PM, Vimal Singh <vimal at mail.utexas.edu> wrote:
> 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;
> }
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list