[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