[Insight-users] Problems with geodesic active contours

[Luis Ibanez]

Hi Anders,


You may want to follow more closely the model illustrated
in the software guide

   http://www . itk . org/ItkSoftwareGuide . pdf



Your pipeline now looks like:

A- CurvatureAnisotropicDiffusion
B- GradientMagnitude(RecursiveGaussian)
C- Threshold Gradient
D- Sigmoid
E- GeodesicActiveContor


The threshold gradient is unecessary, you are sacrifycing the
smoothness of the gradient, and end up with a binary image.
Please remove this filter from the pipeline and connect the
output of the gradient magnitude directly to tha input of the
sigmoid. Like in figure 9.12, pdf-page 365.

The value of beta that you are setting for the sigmoid
seems to be too low.

You will find guidelines for this value on the SoftwareGuide.
Please look at pdf-page 369, paper-page 395. The paragraph
talking about intensity values K1 and K2.

Look a the images in figure 9.13, pdf-page 372.  When you
tune the parameters of the Sigmoid, the output should look
like the image in the lower right of 9.13.  It should be
smoothed but almost binarized, with dark values in the edges
of the structure you want to segment (the aorta in your case).


---

You may want to check your pipeline up to the output of the
FastMarching filter and make sure that upto that level, the
image is reasonable (e.g. a coarse segmentation).  Once you
do this, you can move to refine this segmentation using
the GeodesicActiveContour filter.


-----


I couldn't figure out what you meant by:

> The rescaling intensity is used when I write the result to an image..
> I get good results for my fast marching algorithm, but when I use and with
> the sigmoid function.. but when I put these into GeodesicActiveContour I
> get a result image with maximal pixel values (255 for me..). When I check
> each slice I get 3 of 128 images with a result of an black region on a
> white background(probably inverted result..)


Maybe you  could rephrase this question....


Regards,



   Luis



-------------------------------
Anders Almås wrote:
> Hi!
> I have a problem in implementing geodesic active contours. My problem is
> to segment aorta. Here is my file with typedef's and some other parameters
> that might give you a feeling of what I am doing and what might be wrong:
> 
> // Number of dimensions in ImageType
>   enum{ ImageDimension = 3 };
>  // PixelType is the raw image pixel type
>   typedef float PixelType;
>   // ImageType is the raw image type
>   typedef itk::Image<PixelType, ImageDimension> ImageType;
>   // Derivative filter type, */
>   typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType>
> DerivativeFilterType;
>  //Rescale intensity for increase intnsityscale in images for viewing...
>   typedef itk::RescaleIntensityImageFilter<ImageType, ImageType>
> RescaleRangeType;
>   // SigmoidFilterType for doing a sigmoid filtration
>   typedef itk::SigmoidImageFilter<ImageType,ImageType> SigmoidFilterType;
>   // FastMarchingFilterType for doing fastmarching filtration
>   typedef itk::FastMarchingImageFilter<ImageType,ImageType>
>                                                    FastMarchingFilterType;
>   // GeodesicActiveContourFilterType for doing the geodesicActiveContour
> filtration
>   typedef itk::GeodesicActiveContourLevelSetImageFilter<ImageType,
>                                 ImageType,PixelType>
> GeodesicActiveContourFilterType;
>   // ThresholdType for doing thresholding
>   typedef itk::ThresholdImageFilter<ImageType> ThresholdType;
> 
> std::cout<< "Starts preprocessing"<<std::endl;
>   // Discrete CurvatureAnistropicDiffusionImageFilter
> 
>   typedef itk::CurvatureAnisotropicDiffusionImageFilter<ImageType,
> ImageType> SmoothingFilterType;
>   SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
>   smoothing->SetTimeStep(0.0625); // Default time step for 3D images with
> itk...
>   smoothing->SetIterations(1);    // Testable parameter...
>   smoothing->SetConductanceParameter(1);
>   smoothing->SetInput(input);
>   smoothing->Update();
> 
>   RescaleRangeType::Pointer rangeSmoothing = RescaleRangeType::New();
>   rangeSmoothing->SetOutputMinimum(min_intensity);
>   rangeSmoothing->SetOutputMaximum(max_intensity);
>   rangeSmoothing->SetInput(smoothing->GetOutput());
>   rangeSmoothing->Update();
> 
>   //Setting up a gradientMagnitudeRecursiveGaussianImageFilter
>   DerivativeFilterType::Pointer gradientMagnitude =
> DerivativeFilterType::New();
> 
>   gradientMagnitude->SetSigma(1);
>   gradientMagnitude->SetInput(smoothing->GetOutput());
>   gradientMagnitude->Update();
> 
>   ThresholdType::Pointer thresholdGradient = ThresholdType::New();
>   thresholdGradient->ThresholdBelow(50);
>   thresholdGradient->SetOutsideValue(255);
>   thresholdGradient->SetInput(rangeGradient->GetOutput());
>   thresholdGradient->Update();
> 
>   //Setting up a sigmoid filter
> 
>   SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
> 
>   sigmoid->SetAlpha(-0.6);
>   sigmoid->SetBeta(-0.05);
>   sigmoid->SetInput(thresholdGradient->GetOutput());
>   sigmoid->SetOutputMaximum(1.0);
>   sigmoid->SetOutputMinimum(0.0);
> 
>   //Setting up a fast marching Filter
> 
>   FastMarchingFilterType::Pointer fastMarching =
> FastMarchingFilterType::New();
>   typedef FastMarchingFilterType::NodeType NodeType;
>   typedef FastMarchingFilterType::NodeContainer NodeContainer;
> 
>   NodeType node;
>   NodeType node2;
>   NodeType node3;
> 
>   NodeContainer::Pointer trialpoints = NodeContainer::New();
> 
>   node.SetValue(0.0);
>   node2.SetValue(0.0);
>   node3.SetValue(0.0);
> 
>   ImageType::IndexType seed = {{lumenseed[0], lumenseed[1],
> lumenseed[2]}};
>   ImageType::IndexType seed2 = {{lumenseed2[0], lumenseed2[1],
> lumenseed2[2]}};
>   ImageType::IndexType seed3 = {{lumenseed3[0], lumenseed3[1],
> lumenseed3[2]}};
> 
>   node.SetIndex(seed);
>   node2.SetIndex(seed2);
>   node3.SetIndex(seed3);
> 
>   trialpoints->InsertElement(0, node);
>   trialpoints->InsertElement(1, node2);
>   trialpoints->InsertElement(2, node3);
> 
>   fastMarching->SetTrialPoints(trialpoints);
>   fastMarching->SetSpeedConstant(1.0);
>   fastMarching->SetOutputSize(input->GetBufferedRegion().GetSize());
>   //fastMarching->SetStoppingValue(10.0);
>   fastMarching->Update();
> 
>    // Set GeodesicActiveContourImageFilter
>   GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
> GeodesicActiveContourFilterType::New();
> 
>   geodesicActiveContour->SetPropagationScaling(5);
>   geodesicActiveContour->SetCurvatureScaling(0.1);
>   geodesicActiveContour->SetAdvectionScaling(1.0);
>   geodesicActiveContour->SetMaximumRMSError( 0.02 );
>   geodesicActiveContour->SetMaximumIterations( 100);
>   geodesicActiveContour->SetInput(fastMarching->GetOutput());
>   //geodesicActiveContour->SetFeatureImage(sigmoid->GetOutput());
>   geodesicActiveContour->SetFeatureImage(thresholdGradient->GetOutput());
>   geodesicActiveContour->Update();
> 
>   //Done GeodesicActiveContourImageFilter...
> 
> The rescaling intensity is used when I write the result to an image..
> I get good results for my fast marching algorithm, but when I use and with
> the sigmoid function.. but when I put these into GeodesicActiveContour I
> get a result image with maximal pixel values (255 for me..). When I check
> each slice I get 3 of 128 images with a result of an black region on a
> white background(probably inverted result..)
> What is then wrong here? Is it that my parameters are wrong or the code it
> self?
> 
> 
> 
> -------------------------------------------------------------------------------
> Navn:   Anders Almås
> E-post: andersal at stud . ntnu . no
> Tlf:    90550374
> -------------------------------------------------------------------------------
> 
>