Hello,<br><br>I am having a similar problem with the geodesicactivecontourfilter.<br><br>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.
<br><br>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.
<br><br>Thanks for any help!<br><br><div><span class="gmail_quote">On 3/23/07, <b class="gmail_sendername">Luis Ibanez</b> <<a href="mailto:luis.ibanez@kitware.com">luis.ibanez@kitware.com</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi Oliver,<br><br>The first place to look is at the Speed image that you are feeding<br>into the GeodesicActive contour filter.<br><br>You must make sure that this image is as close as possible to a<br>binary image, with zero values in the edges where you want the
<br>level set to stop.<br><br>You may find convenient to split your program in two.<br><br>The first half just to compute the speed image as the output<br>of the sigmoid filter.<br><br>The second half reading the speed image from a file and
<br>continuig from there.<br><br>The fact that your contours are not stopping at the edges of the<br>image is an indication that your speed image doesn't have a strong<br>enough contrast.<br><br><br> Regards,<br><br>
<br><br> Luis<br><br><br><br>---------------------<br>Oliver Gloger wrote:<br>> Hello Luis,<br>><br>> 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?
<br>><br>> 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() );
<br>> The OutputImage of the gardientMagnitude is okay!!! But my segments are again too large! Can you notice what is wrong?<br>><br>> 2.) I do not know my initialdistance for the contours before! Is it really necessary? I thought this filter orientates at the edge potential?
<br>><br>> I would be very thankful, if you will answer me!<br>><br>> Best regards,<br>> Oliver!<br>><br>> Here is the code:<br>><br>> #include "itkImage.h"<br>> #include "itkGeodesicActiveContourLevelSetImageFilter.h
"<br>> #include "itkCurvatureAnisotropicDiffusionImageFilter.h"<br>> #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"<br>> #include "itkSigmoidImageFilter.h"<br>> #include "
itkFastMarchingImageFilter.h"<br>> #include "itkRescaleIntensityImageFilter.h"<br>> #include "itkBinaryThresholdImageFilter.h"<br>> #include "itkImageFileReader.h"<br>> #include "
itkImageFileWriter.h"<br>><br>><br>> int main( int argc, char *argv[] )<br>> {<br>> if( argc < 9 )<br>> {<br>> std::cerr << "Missing Parameters " << std::endl;<br>
> std::cerr << "Usage: " << argv[0];<br>> std::cerr << " inputImageName boundaryImageName outputImageName";<br>> std::cerr << " SigmoidAlpha SigmoidBeta initialdistance";
<br>> std::cerr << " pairwise seedpoints (x1) (y1) (x2) (y2) ..." << std::endl;<br>> return 1;<br>> }<br>><br>> typedef float InternalPixelType;<br>> const unsigned int Dimension = 2;
<br>> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;<br>><br>> typedef unsigned char OutputPixelType;<br>> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
<br>> typedef itk::BinaryThresholdImageFilter<<br>> InternalImageType,<br>> OutputImageType > ThresholdingFilterType;<br>><br>> ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
<br>><br>> thresholder->SetLowerThreshold( -1000.0 );<br>> thresholder->SetUpperThreshold( 0.0 );<br>><br>> thresholder->SetOutsideValue( 0 );<br>> thresholder->SetInsideValue( 255 );
<br>><br>> typedef itk::ImageFileReader< InternalImageType > ReaderType;<br>> typedef itk::ImageFileWriter< OutputImageType > WriterType;<br>><br>> ReaderType::Pointer reader = ReaderType::New();
<br>> WriterType::Pointer writer = WriterType::New();<br>><br>> reader->SetFileName( argv[1] );<br>> writer->SetFileName( argv[3] );<br>><br>> typedef itk::RescaleIntensityImageFilter<<br>
> InternalImageType,<br>> OutputImageType > CastFilterType;<br>><br>><br>><br>> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
<br>> InternalImageType,<br>> InternalImageType > GradientFilterType;<br>><br>> typedef itk::SigmoidImageFilter<<br>> InternalImageType,
<br>> InternalImageType > SigmoidFilterType;<br>><br>> GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();<br>><br>> SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
<br>><br>> sigmoid->SetOutputMinimum( 0.0 );<br>> sigmoid->SetOutputMaximum( 1.0 );<br>> typedef itk::FastMarchingImageFilter<<br>> InternalImageType,<br>> InternalImageType > FastMarchingFilterType;
<br>><br>> FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();<br>><br>> typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,<br>> InternalImageType > GeodesicActiveContourFilterType;
<br>> GeodesicActiveContourFilterType::Pointer geodesicActiveContour =<br>> GeodesicActiveContourFilterType::New();<br>><br>> const double propagationScaling = 2.0;<br>
> geodesicActiveContour->SetPropagationScaling( propagationScaling );<br>> geodesicActiveContour->SetCurvatureScaling( 1.0 );<br>> geodesicActiveContour->SetAdvectionScaling( 2.0 );<br>> geodesicActiveContour->SetMaximumRMSError(
0.02 );<br>> geodesicActiveContour->SetNumberOfIterations( 50 );<br>><br>> gradientMagnitude->SetInput( reader->GetOutput() );<br>> sigmoid->SetInput( gradientMagnitude->GetOutput() );<br>
><br>> geodesicActiveContour->SetInput( fastMarching->GetOutput() );<br>> geodesicActiveContour->SetFeatureImage( gradientMagnitude->GetOutput() );<br>><br>> thresholder->SetInput( geodesicActiveContour->GetOutput() );
<br>> writer->SetInput( thresholder->GetOutput() );<br>><br>> const double sigma = 0.75;<br>> gradientMagnitude->SetSigma( sigma );<br>> sigmoid->SetAlpha( atof(argv[4]) );<br>> sigmoid->SetBeta( atof(argv[5]) );
<br>><br>> typedef FastMarchingFilterType::NodeContainer NodeContainer;<br>> typedef FastMarchingFilterType::NodeType NodeType;<br>><br>> NodeContainer::Pointer seeds = NodeContainer::New();<br>
><br>> InternalImageType::IndexType seedPosition;<br>> NodeType node;<br>> const double initialDistance = atof(argv[6]);<br>> const double seedValue = - initialDistance;<br>><br>> //I use a boundaryImage too visualize my boundaries I use from Matlab for //every segment
<br>> OutputImageType::Pointer boundaryimage = OutputImageType::New();<br>> OutputImageType::IndexType start;<br>> start[0] = 0; // first index on X<br>> start[1] = 0; // first index on Y<br>><br>> OutputImageType::SizeType size;
<br>> size[0] = 256; // size along X<br>> size[1] = 256; // size along Y<br>><br>> OutputImageType::RegionType region;<br>> region.SetSize( size );<br>> region.SetIndex( start );<br>> boundaryimage->SetRegions( region );
<br>> boundaryimage->Allocate();<br>><br>> OutputImageType::IndexType pixelIndex;<br>> //pixelIndex[0] = 27; // x position<br>> //pixelIndex[1] = 29; // y position<br>> //node.SetValue( seedValue );
<br>> //node.SetIndex( seedPosition );<br>> seeds->Initialize();<br>> int saatpaare=7;<br>> int nodeindex=0;<br>> while (saatpaare<argc){<br>> seedPosition[0] = atoi(argv[saatpaare]);
<br>> seedPosition[1] = atoi(argv[saatpaare+1]);<br>> pixelIndex[0] = atoi(argv[saatpaare]);<br>> pixelIndex[1] = atoi(argv[saatpaare+1]);<br>> boundaryimage->SetPixel( pixelIndex , 100 );
<br>> node.SetValue( seedValue );<br>> node.SetIndex( seedPosition );<br>> seeds->InsertElement( nodeindex, node );<br>> saatpaare=saatpaare+2;<br>> nodeindex++;
<br>> }<br>><br>> typedef itk::ImageFileWriter<OutputImageType> BoundaryWriterType;<br>> BoundaryWriterType::Pointer boundaryWriter = BoundaryWriterType::New();<br>> boundaryWriter->SetFileName(argv[2]);
<br>> boundaryWriter->SetInput(boundaryimage);<br>> try<br>> {<br>> boundaryWriter->Update();<br>> }<br>> catch( itk::ExceptionObject & excep )<br>> {<br>> std::cerr << "Exception caught at boundayrWriter!" << std::endl;
<br>> std::cerr << excep << std::endl;<br>> }<br>><br>> fastMarching->SetTrialPoints( seeds );<br>><br>> fastMarching->SetSpeedConstant( 1.0 );<br>> CastFilterType::Pointer caster2 = CastFilterType::New();
<br>> CastFilterType::Pointer caster3 = CastFilterType::New();<br>> CastFilterType::Pointer caster4 = CastFilterType::New();<br>><br>> WriterType::Pointer writer2 = WriterType::New();<br>> WriterType::Pointer writer3 = WriterType::New();
<br>> WriterType::Pointer writer4 = WriterType::New();<br>><br>> caster2->SetInput( gradientMagnitude->GetOutput() );<br>> writer2->SetInput( caster2->GetOutput() );<br>> writer2->SetFileName("
GeodesicActiveContourImageFilterOutput2.png");<br>> caster2->SetOutputMinimum( 0 );<br>> caster2->SetOutputMaximum( 255 );<br>> writer2->Update();<br>><br>> // This Output is always empty!!!!!!! But why????
<br>> caster3->SetInput( sigmoid->GetOutput() );<br>> writer3->SetInput( caster3->GetOutput() );<br>> writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");<br>> caster3->SetOutputMinimum( 0 );
<br>> caster3->SetOutputMaximum( 255 );<br>> writer3->Update();<br>><br>> caster4->SetInput( fastMarching->GetOutput() );<br>> writer4->SetInput( caster4->GetOutput() );<br>> writer4->SetFileName("
GeodesicActiveContourImageFilterOutput4.png");<br>> caster4->SetOutputMinimum( 0 );<br>> caster4->SetOutputMaximum( 255 );<br>> fastMarching->SetOutputSize(<br>> reader->GetOutput()->GetBufferedRegion().GetSize() );
<br>> try<br>> {<br>> writer->Update();<br>> }<br>> catch( itk::ExceptionObject & excep )<br>> {<br>> std::cerr << "Exception caught !" << std::endl;
<br>> std::cerr << excep << std::endl;<br>> }<br>><br>> writer4->Update();<br>><br>> return 0;<br>> }<br>><br>><br>><br>><br>><br>> ------------------------------------------------------------------------
<br>><br>_______________________________________________<br>Insight-users mailing list<br><a href="mailto:Insight-users@itk.org">Insight-users@itk.org</a><br><a href="http://www.itk.org/mailman/listinfo/insight-users">
http://www.itk.org/mailman/listinfo/insight-users</a><br></blockquote></div><br>