[Insight-users] Re: help for GeodesicActiveContourImageFilter
Luis Ibanez
luis.ibanez at kitware.com
Sat Mar 31 11:53:37 EST 2007
Hi d.f.
Please post the min and max values of your speed image.
If the level set is not stopping at the location of
edges in the speed image, the following maybe the
typical suspects:
1) Your propagation weighting parameter is too high.
2) You are running for too many iterations.
Given enough time (iterations) a level set will always
leak to the other sides of contours, unless your speed
image really have zeros all around the contour.
You may want to add a visualization layer to your
program. Add an Observer to the GeodesicActiveContour
filter and every N iterations use the vtkContour filter
and associated classes for displaying the zero set of
the level set. That will let you see how the level set
is progressing as the iterations go by.
For example code on how use VTK for displaying the
contour of an ITK 3D image, please look at:
InsightApplications/Auxiliary/vtk/
itkReadITKImage3DSegmentShowVTK.cxx
Regards,
Luis
-------------
d f wrote:
> 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
> <mailto: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 <mailto:Insight-users at itk.org>
> http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list