[Insight-users] Re: help for GeodesicActiveContourImageFilter

Luis Ibanez luis.ibanez at kitware.com
Fri Mar 23 17:47:34 EST 2007


Hi Oliver,

The first place to look is at the Speed image that you are feeding
into the GeodesicActive contour filter.

You must make sure that this image is as close as possible to a
binary image, with zero values in the edges where you want the
level set to stop.

You may find convenient to split your program in two.

The first half just to compute the speed image as the output
of the sigmoid filter.

The second half reading the speed image from a file and
continuig from there.

The fact that your contours are not stopping at the edges of the
image is an indication that your speed image doesn't have a strong
enough contrast.


    Regards,



       Luis



---------------------
Oliver Gloger wrote:
> Hello Luis,
> 
> I send 2 e-mails to the insight-users mailing list, but got no response for my questions for the GeodesicActiveContourImageFilter. Perhaps you could give an answer for my questions?
> 
> 1.) My segments are growing too large (see attachment)! It seems that this filter does not notice the edge potential. Is it necessary to use the sigmoidfilter? After using my sigmoidfilter my output is a white picture!!! But I don't know why!!! Could that be a problem? Isn't it enough to use the GradientMagnitude directly for the edge potential? So I made the following change in the Code in your Examples: geodesicActiveContour->SetFeatureImage( gradientMagnitude->GetOutput() );
> The OutputImage of the gardientMagnitude is okay!!! But my segments are again too large! Can you notice what is wrong? 
> 
> 2.) I do not know my initialdistance for the contours before! Is it really necessary? I thought this filter orientates at the edge potential?
> 
> I would be very thankful, if you will answer me!
> 
> Best regards,
> Oliver!
> 
> Here is the code:
> 
> #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"
> 
> 
> int main( int argc, char *argv[] )
> {
>   if( argc < 9 )
>     {
>     std::cerr << "Missing Parameters " << std::endl;
>     std::cerr << "Usage: " << argv[0];
>     std::cerr << " inputImageName  boundaryImageName outputImageName";
>     std::cerr << " SigmoidAlpha SigmoidBeta initialdistance";
>     std::cerr << " pairwise seedpoints (x1) (y1) (x2) (y2) ..."  << std::endl;
>     return 1;
>     }
>   
>   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::ImageFileReader< InternalImageType > ReaderType;
>   typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
> 
>   ReaderType::Pointer reader = ReaderType::New();
>   WriterType::Pointer writer = WriterType::New();
> 
>   reader->SetFileName( argv[1] );
>   writer->SetFileName( argv[3] );
> 
>   typedef itk::RescaleIntensityImageFilter< 
>                                InternalImageType, 
>                                OutputImageType >   CastFilterType;
> 
> 
>   
>   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::FastMarchingImageFilter< 
>                               InternalImageType, 
>                               InternalImageType >    FastMarchingFilterType;
> 
>   FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
>   
>   typedef  itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, 
>                 InternalImageType >    GeodesicActiveContourFilterType;
>   GeodesicActiveContourFilterType::Pointer geodesicActiveContour = 
>                                      GeodesicActiveContourFilterType::New();
> 
>   const double propagationScaling = 2.0;
>   geodesicActiveContour->SetPropagationScaling( propagationScaling );
>   geodesicActiveContour->SetCurvatureScaling( 1.0 );
>   geodesicActiveContour->SetAdvectionScaling( 2.0 );
>   geodesicActiveContour->SetMaximumRMSError( 0.02 );
>   geodesicActiveContour->SetNumberOfIterations( 50 );
> 
>   gradientMagnitude->SetInput( reader->GetOutput() );
>   sigmoid->SetInput( gradientMagnitude->GetOutput() );
> 
>   geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
>   geodesicActiveContour->SetFeatureImage( gradientMagnitude->GetOutput() );
> 
>   thresholder->SetInput( geodesicActiveContour->GetOutput() );
>   writer->SetInput( thresholder->GetOutput() );
>   
>   const double sigma = 0.75;
>   gradientMagnitude->SetSigma(  sigma  );
>   sigmoid->SetAlpha( atof(argv[4]) );
>   sigmoid->SetBeta(  atof(argv[5]) );
>   
>   typedef FastMarchingFilterType::NodeContainer  NodeContainer;
>   typedef FastMarchingFilterType::NodeType       NodeType;
> 
>   NodeContainer::Pointer seeds = NodeContainer::New();
> 
>   InternalImageType::IndexType  seedPosition;
> 	NodeType node;
>   const double initialDistance = atof(argv[6]);
>   const double seedValue = - initialDistance;
> 
>   //I use a boundaryImage too visualize my boundaries I use from Matlab for //every segment
>   OutputImageType::Pointer boundaryimage = OutputImageType::New();
>   OutputImageType::IndexType start;
>   start[0] = 0; // first index on X
>   start[1] = 0; // first index on Y
> 
> 	OutputImageType::SizeType size;
> 	size[0] = 256; // size along X
> 	size[1] = 256; // size along Y
> 
>   OutputImageType::RegionType region;
>   region.SetSize( size );
>   region.SetIndex( start );
> boundaryimage->SetRegions( region );
> boundaryimage->Allocate();
> 
> OutputImageType::IndexType pixelIndex;
> //pixelIndex[0] = 27; // x position
> //pixelIndex[1] = 29; // y position
>   //node.SetValue( seedValue );
>   //node.SetIndex( seedPosition );
>       seeds->Initialize();
> 	  int saatpaare=7;
> 	  int nodeindex=0;
> 	  while (saatpaare<argc){
> 		seedPosition[0] = atoi(argv[saatpaare]);
> 		seedPosition[1] = atoi(argv[saatpaare+1]);
> 		pixelIndex[0] = atoi(argv[saatpaare]);
> 		pixelIndex[1] = atoi(argv[saatpaare+1]);
> 		boundaryimage->SetPixel( pixelIndex , 100 );
> 		node.SetValue( seedValue );
> 		node.SetIndex( seedPosition );
>         seeds->InsertElement( nodeindex, node );
> 		saatpaare=saatpaare+2;
> 		nodeindex++;
> 		}
> 
> typedef itk::ImageFileWriter<OutputImageType> BoundaryWriterType;
>   BoundaryWriterType::Pointer boundaryWriter = BoundaryWriterType::New();
>   boundaryWriter->SetFileName(argv[2]);
>   boundaryWriter->SetInput(boundaryimage);
>   try
>     {
>     boundaryWriter->Update();
>     }
>   catch( itk::ExceptionObject & excep )
>     {
>     std::cerr << "Exception caught at boundayrWriter!" << std::endl;
>     std::cerr << excep << std::endl;
>     }
>   
>   fastMarching->SetTrialPoints(  seeds  );
> 
>   fastMarching->SetSpeedConstant( 1.0 );
>   CastFilterType::Pointer caster2 = CastFilterType::New();
>   CastFilterType::Pointer caster3 = CastFilterType::New();
>   CastFilterType::Pointer caster4 = CastFilterType::New();
> 
>   WriterType::Pointer writer2 = WriterType::New();
>   WriterType::Pointer writer3 = WriterType::New();
>   WriterType::Pointer writer4 = WriterType::New();
> 
>   caster2->SetInput( gradientMagnitude->GetOutput() );
>   writer2->SetInput( caster2->GetOutput() );
>   writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");
>   caster2->SetOutputMinimum(   0 );
>   caster2->SetOutputMaximum( 255 );
>   writer2->Update();
> 
> // This Output is always empty!!!!!!! But why????
>   caster3->SetInput( sigmoid->GetOutput() );
>   writer3->SetInput( caster3->GetOutput() );
>   writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");
>   caster3->SetOutputMinimum(   0 );
>   caster3->SetOutputMaximum( 255 );
>   writer3->Update();
> 
>   caster4->SetInput( fastMarching->GetOutput() );
>   writer4->SetInput( caster4->GetOutput() );
>   writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");
>   caster4->SetOutputMinimum(   0 );
>   caster4->SetOutputMaximum( 255 );
>   fastMarching->SetOutputSize( 
>            reader->GetOutput()->GetBufferedRegion().GetSize() );
>    try
>     {
>     writer->Update();
>     }
>   catch( itk::ExceptionObject & excep )
>     {
>     std::cerr << "Exception caught !" << std::endl;
>     std::cerr << excep << std::endl;
>     }
>   
>   writer4->Update();
> 
>   return 0;
> }
> 
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 


More information about the Insight-users mailing list