[Insight-users] Active Contour not Propagating

[Andrey Sklyar]

I am trying to use GeodesicActiveContourLevelSetImageFilter to extract
the boundary of the lungs from a series of sagittal MRI slices through
a patient.  For testing purposes, I'm only using 6 slices.

After extracting an initial rough boundary of the lungs, I generate a
signed distance map using SignedDanielssonDistanceMapImageFilter to
use as the initial level set to the active contour filter.  The
feature image is generated using the same pipeline as in Software
Guide, pg 537, and is fully defined in the code excerpt below.

The problem that I am running into is that the boundary does not
expand even though the propagation force is positive, and the
speed/feature image has values close to 1 around the boundary.  The
curvature force has a visible effect though - the boundary is smoothed
out.

I can think of a couple of possible causes of these problems, but I am
not intimately familiar with the implementation of these filters, so I
would like to get the advice of others who are more experienced in
their use than I am.  Here are my hypotheses:

- The image is too skinny in the z-dimension - the resolution in each
xy slice is high (~128x128), but I'm only using 6 images in the z
direction.
- The voxels have a different size in the z-dimension (the slice
thickness is large compared to the in-slice pixels' size) - does the
active contour filter expect the voxels to be the same size in all
dimensions?
- The contour is touching the image boundary in the z-dimesion and is
stopping because of this - I ran some experiments in 2D to check this
case, but the boundary still keeps evolving even after it hits the
border, so this is unlikely to be the problem

I would also not be surprised if it were a completely different
problem all together.  Any insight into why the propagation scaling
would not be having any effect on the boundary evolution would be of
great help!

- Andrey

// ============================== Code Excerpt ================================
// ===== You can copy and paste this into a text editor for easier viewing ====
// ============================================================================
// Seed Region Extraction #include statements
// ...
#include <itkBinaryThresholdImageFilter.h>

// Active Contour Include Statements
#include <itkSignedDanielssonDistanceMapImageFilter.h>
#include <itkCurvatureAnisotropicDiffusionImageFilter.h>
#include <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
#include <itkSigmoidImageFilter.h>
#include <itkGeodesicActiveContourLevelSetImageFilter.h>

class Segmenter
{
public:
    typedef signed short PixelType;
    typedef float FloatPixelType;

    // Setter Methods for Filter Parameters
    // ...

    // Getter Methods for Intermediate Images in the Pipeline
    // ...

private:
    // Seed Region Extraction Types
    // ...
    typedef itk::Image<PixelType, num_dims>            ImageType;

    typedef itk::BinaryThresholdImageFilter<
        ImageType, ImageType>                        Extractor;

    // Active Contour Types
    typedef itk::Image<
        FloatPixelType, num_dims>                     FloatImageType;

    typedef itk::SignedDanielssonDistanceMapImageFilter<
        ImageType, FloatImageType>                    ContourInitializer;

    typedef itk::CurvatureAnisotropicDiffusionImageFilter<
        ImageType, FloatImageType>                    Smoother;

    typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
        FloatImageType, FloatImageType>                Gradient;

    typedef itk::SigmoidImageFilter<
        FloatImageType, FloatImageType>              Sigmoid;

    typedef itk::GeodesicActiveContourLevelSetImageFilter<
        FloatImageType, FloatImageType>                ActiveContourFilter;

    typedef itk::BinaryThresholdImageFilter<
        FloatImageType, ImageType>                    ContourExtractor;

    // Seed Extraction Variables
    // ...
    Extractor::Pointer extractor;

    // Active Contour Variables
    ContourInitializer::Pointer initializer;
    Smoother::Pointer smoother;
    Gradient::Pointer gradient;
    Sigmoid::Pointer sigmoid;
    ActiveContourFilter::Pointer activeContourFilter;
    ContourExtractor::Pointer contourExtractor;
};

Segmenter::Segmenter()
{
    // Initializers
    initializer = ContourInitializer::New();
    smoother = Smoother::New();
    gradient = Gradient::New();
    sigmoid = Sigmoid::New();
    activeContourFilter = ActiveContourFilter::New();
    contourExtractor = ContourExtractor::New();

    // Parameters
    smoother->SetTimeStep( 0.0625 );
    smoother->SetConductanceParameter( 9.0 );

    gradient->SetSigma(1.0);

    sigmoid->SetOutputMinimum(0.00001);
    sigmoid->SetOutputMaximum(1.0);

    sigmoid->SetAlpha(-5.0);
    sigmoid->SetBeta(100.0);

    activeContourFilter->SetPropagationScaling(1.0);
    activeContourFilter->SetAdvectionScaling(1.0);
    activeContourFilter->SetCurvatureScaling(0.1);
    activeContourFilter->SetMaximumRMSError(0.05);
    activeContourFilter->SetNumberOfIterations(NUM_ITERS);

    contourExtractor->SetLowerThreshold( -1000.0 );
    contourExtractor->SetUpperThreshold( 0.0 );
    contourExtractor->SetOutsideValue(0);
    contourExtractor->SetInsideValue(thresholdValue);

    // Pipeline
    initializer->SetInput(extractor->GetOutput());
    smoother->SetInput(reader->GetOutput());
    gradient->SetInput(smoother->GetOutput());
    sigmoid->SetInput(gradient->GetOutput());
    activeContourFilter->SetInput(initializer->GetOutput());
    activeContourFilter->SetFeatureImage(sigmoid->GetOutput());
    contourExtractor->SetInput(activeContourFilter->GetOutput());
}

// Segmenter Setter Function Definitions
// ...

// Segmenter Image Getter Function Definitions
// ...