[Insight-users] Segmentation simple circular model with level sets

[silvano]

Hi all,

I would like to use the very usefull ITK classes to segment mammographic
images with level set methods.
In order to learn these classes, I run properly the examples provided with
the library and explained in the software guide.
Therefore, I've tried to apply the filter ShapeDetectionLevelSetImageFilter
to a very simple circular model,
but I cannot obtain any reasonable segmentation.

Hereafter, I report the code used (sorry for the long list of code lines!!
:) ).
The application accepts a list of arguments, precisely :
1_size_model 2_amplitude_model 3_sigma_smoothing 4_propagationscale 5
curvaturescale 6_levelset_initialvalue
I tried a lot of different argument lists, among them you can try the
following one
101 10 1 1 1 10

I would very appreciate to receive any comment, suggestion to solve the
problem.

Thank u in advance,

Silvano


#include <itkIndex.h>
#include <itkImage.h>
#include <itkImageFileWriter.h>
#include <itkCastImageFilter.h>
#include <itkTimeProbesCollectorBase.h>
#include <itkGradientAnisotropicDiffusionImageFilter.h>
#include <itkFastMarchingImageFilter.h>
#include <itkShapeDetectionLevelSetImageFilter.h>
#include <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
#include <itkBoundedReciprocalImageFilter.h>

typedef itk::Image<float,2> TImage;

TImage::Pointer model2D_flat( float sizeM, float amp )
{
    TImage::Pointer model = TImage::New();

    int hsize = (int)((sizeM-1.0)/2.0);

    std::cout << "size model image = " << sizeM;
    std::cout << ", half size model image = " << hsize << std::endl;

    TImage::IndexType start;
    start[0] =   0;
    start[1] =   0;

    TImage::SizeType  size;
    size[0]  = sizeM;
    size[1]  = sizeM;

    TImage::RegionType region;
    region.SetSize( size );
    region.SetIndex( start );

    model->SetRegions( region );
    model->Allocate();

    float radius2 = static_cast<float>(sizeM)/4.0;
    radius2 = radius2 * radius2;
    for ( int y=-hsize; y <= hsize; ++y )
    {
        for ( int x=-hsize; x <= hsize; ++x )
        {
            TImage::IndexType index;
            index[0] = (x+hsize);
            index[1] = (y+hsize);
            if ( (x*x+y*y) < radius2 )
                model->SetPixel( index, amp );
            else
                model->SetPixel( index, amp/2.0 );
        }
    }

    return model;
}


int main(int argc, char **argv)
{
    double sizeM = atof( argv[1] );
    double amplitude = atof( argv[2] );
    double sigma = atof( argv[3] );
    double propscale = atof( argv[4] );
    double curvscale = atof( argv[5] );
    double isovalue = atof( argv[6] );

    std::cout << "size= " << sizeM << std::endl;
    std::cout << "amplitude= " << amplitude << std::endl;
    std::cout << "sigma= " << sigma << std::endl;
    std::cout << "propscale= " << propscale << std::endl;
    std::cout << "curvscale= " << curvscale << std::endl;
    std::cout << "isovalue= " << isovalue << std::endl;

    TImage::Pointer model = model2D_flat( sizeM, amplitude );

    itk::TimeProbesCollectorBase timeCollector;
    timeCollector.Start("Total");

    typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< TImage,
TImage > GradientFilterType;
    GradientFilterType::Pointer gradient = GradientFilterType::New();
    gradient->SetSigma( sigma );
    gradient->SetNormalizeAcrossScale(false);

    typedef itk::BoundedReciprocalImageFilter< TImage, TImage >
BoundFilterType;
    BoundFilterType::Pointer bounder = BoundFilterType::New();

    typedef itk::FastMarchingImageFilter< TImage, TImage >
FastMarchingFilterType;
    FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();

    typedef  itk::ShapeDetectionLevelSetImageFilter< TImage, TImage >
ShapeDetectionFilterType;
    ShapeDetectionFilterType::Pointer shapeDetection =
ShapeDetectionFilterType::New();

    gradient->SetInput( model );
    bounder->SetInput( gradient->GetOutput() );
    shapeDetection->SetInput( fastMarching->GetOutput() );
    shapeDetection->SetFeatureImage( bounder->GetOutput() );

    typedef FastMarchingFilterType::NodeContainer        NodeContainer;
    typedef FastMarchingFilterType::NodeType             NodeType;
    NodeContainer::Pointer seeds = NodeContainer::New();

    NodeType node;
    const double seedValue = -isovalue;


    itk::Index<2> pos;
    int hsize = (int)((sizeM-1.0)/2.0);
    pos[0] = hsize; pos[1] = hsize;
    node.SetValue( seedValue );
    node.SetIndex( pos );

    seeds->Initialize();
    seeds->InsertElement( 0, node );

    fastMarching->SetTrialPoints( seeds );
    fastMarching->SetOutputSize( model->GetLargestPossibleRegion().GetSize()
);
    fastMarching->SetOutputSpacing( model->GetSpacing() );
    fastMarching->SetOutputOrigin ( model->GetOrigin() );
    fastMarching->SetSpeedConstant( 1.0 );

    shapeDetection->SetPropagationScaling( propscale );
    shapeDetection->SetCurvatureScaling( curvscale );
    shapeDetection->SetMaximumRMSError( 0.02 );
    shapeDetection->SetNumberOfIterations( 1000 );

    typedef itk::ImageFileWriter<TImage> WriterType;
    WriterType::Pointer iwriter = WriterType::New();

    timeCollector.Start("process");
    iwriter->SetInput( model );
    iwriter->SetFileName( "model.nii" );
    iwriter->Update();

    iwriter->SetInput( fastMarching->GetOutput() );
    iwriter->SetFileName( "fastMarchOutput.nii" );
    iwriter->Update();

    iwriter->SetInput( shapeDetection->GetOutput() );
    iwriter->SetFileName( "shapeOutput.nii" );
    iwriter->Update();

    iwriter->SetInput( bounder->GetOutput() );
    iwriter->SetFileName( "edge.nii" );
    iwriter->Update();

    timeCollector.Stop("process");

    std::cout << "No. elpased iterations: " <<
shapeDetection->GetElapsedIterations() << std::endl;
    std::cout << "RMS change: " << shapeDetection->GetRMSChange() <<
std::endl;

    timeCollector.Stop("Total");
    std::cout << std::endl;
    std::cout << "Computing Time (sec): " << std::endl;
    timeCollector.Report();
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100929/7b89166e/attachment.htm>