[Insight-users] Re: help for GeodesicActiveContourImageFilter

d f drf8888 at gmail.com
Thu Mar 29 20:51:47 EST 2007


Hello,

I am having a similar problem with the geodesicactivecontourfilter.

I am using ITK 3.0.0 with cygwin.  I modified the program to handle 3D
images, but the rest is the same. and I can monitor the outputs of each
stage (smoothing, gradient, sigmoid, fastmarching).  They are all correct.

However, the speed image seems to have no effect on the segmentation- the
surface propagates right through the regions of low values. It is as if the
level set function is ignoring the speed image completely- has anyone else
run into this? I am able to correctly return the speed image back from the
geodesicActiveContour with GetfeatureImage(), so I know its at least getting
there correctly.

Thanks for any help!

On 3/23/07, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
> 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;
> > }
> >
> >
> >
> >
> >
> > ------------------------------------------------------------------------
> >
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070329/0151c077/attachment.html


More information about the Insight-users mailing list