[ITK] [ITK-users] Motion under Mean Curvature with LevelSetv4 Framework

Cagatay Bilgin bilgincc at gmail.com
Tue Jul 22 03:58:17 EDT 2014


Dear ITK Community,

I am trying to familiarize myself with the new level set classes. My goal
is to simulate motion under mean curvature using the new design.

I have put together the following, looking at the examples and tests.
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3872740/ is a great source to
understand the design and tests in source directory are very helpful, thank
you for the resources!

The test scenario I have is a L shaped object. I would expect the object to
become somewhat ellipse and then disappear at the end of the evolution.
However I get the following results attached to the email. I don't quite
follow these results. Am I missing something obvious?

Here is the code and CMakeLists.txt. I ran the program with ./Motion 100 0
0 0.05

#include "itkBinaryImageToLevelSetImageAdaptor.h"
#include "itkImageFileReader.h"
#include "itkLevelSetContainer.h"
#include "itkLevelSetEquationPropagationTerm.h"
#include "itkLevelSetEquationAdvectionTerm2.h"
#include "itkLevelSetEquationContainer.h"
#include "itkLevelSetEquationTermContainer.h"
#include "itkLevelSetEvolution.h"
#include "itkLevelSetEvolutionNumberOfIterationsStoppingCriterion.h"
#include "itkLevelSetDenseImage.h"
#include "itkVTKVisualizeImageLevelSetIsoValues.h"
#include "itkSinRegularizedHeavisideStepFunction.h"
#include "itkLevelSetIterationUpdateCommand.h"
#include "itkLevelSetEquationCurvatureTerm.h"
#include "itkCastImageFilter.h"
#include "itkWhitakerSparseLevelSetImage.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"

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

/*
 * L Shape
 */
void CreateSquareImage(FIT::Pointer image)
{
    FIT::RegionType region;
    FIT::IndexType start;
    start[0] = 0;
    start[1] = 0;

    FIT::SizeType size;
    size[0] = 200;
    size[1] = 300;

    region.SetSize(size);
    region.SetIndex(start);

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

    // Set pixels in a square to one value
    for(unsigned int r = 20; r < 160; r++)
    {
        for(unsigned int c = 30; c < 100; c++)
        {
            FIT::IndexType pixelIndex;
            pixelIndex[0] = r;
            pixelIndex[1] = c;

            image->SetPixel(pixelIndex, 255);
        }
    }

    for(unsigned int r = 20; r < 80; r++)
    {
        for(unsigned int c = 100; c < 200; c++)
        {
            FIT::IndexType pixelIndex;
            pixelIndex[0] = r;
            pixelIndex[1] = c;

            image->SetPixel(pixelIndex, 255);
        }
    }
}


/*
 */
int main( int argc, char* argv[] )
{
    if( argc != 5)
    {
        std::cerr << "Missing Arguments" << std::endl;
        std::cerr << argv[0] << std::endl;
        std::cerr << "1- Number of Iterations" << std::endl;
        std::cerr << "2- Propagation Term" << std::endl;
        std::cerr << "3- Advection Term" << std::endl;
        std::cerr << "4- Curvature Term" << std::endl;
        return EXIT_FAILURE;
    }

    // Image Dimension
    const unsigned int Dimension = 2;

    typedef unsigned char                            InputPixelType;
    typedef itk::Image< InputPixelType, Dimension >  IIT;
    typedef itk::Image< float, 2 >  FIT;

    FIT::Pointer input = FIT::New();
    CreateSquareImage(input);

    int numberOfIterations = atoi( argv[1]);

    typedef float                                         LevelSetPixelType;
    typedef itk::Image< LevelSetPixelType, Dimension >    LSIT;
    typedef itk::LevelSetDenseImage< LSIT >  LST;
    //typedef itk::WhitakerSparseLevelSetImage< LevelSetPixelType, 2 > LST;

    typedef LST::OutputType                      LevelSetOutputType;
    typedef LST::OutputRealType                  LevelSetRealType;

    // convert a binary mask to a level-set function
    typedef itk::BinaryImageToLevelSetImageAdaptor<FIT, LST >  BI2LST;
    BI2LST::Pointer adaptor = BI2LST::New();
    adaptor->SetInputImage( input );
    adaptor->Initialize();
    LST::Pointer levelSet = adaptor->GetLevelSet();

    // The Heaviside function
    typedef itk::SinRegularizedHeavisideStepFunction< LevelSetRealType,
LevelSetRealType > HeavisideFunctionType;
    HeavisideFunctionType::Pointer heaviside = HeavisideFunctionType::New();
    heaviside->SetEpsilon( 1 );

    // Create the level set container
    typedef itk::LevelSetContainer< itk::IdentifierType, LST > LSContainerT;
    LSContainerT::Pointer levelSetContainer = LSContainerT::New();
    levelSetContainer->SetHeaviside( heaviside );
    levelSetContainer->AddLevelSet( 0, levelSet );

    // Create the terms.
    typedef itk::LevelSetEquationPropagationTerm<FIT, LSContainerT>
 PropagationTermType;
    PropagationTermType::Pointer propagationTerm =
PropagationTermType::New();
    propagationTerm->SetInput(input);
    propagationTerm->SetCoefficient(atof(argv[2]));

    typedef itk::LevelSetEquationAdvectionTerm2<FIT, LSContainerT>
 AdvectionTermType;
    AdvectionTermType::Pointer advectionTerm = AdvectionTermType::New();
    advectionTerm->SetInput(input);
    advectionTerm->SetCoefficient(atof(argv[3]));

    typedef itk::LevelSetEquationCurvatureTerm<FIT, LSContainerT>
 CurvatureTermType;
    CurvatureTermType::Pointer curvatureTerm = CurvatureTermType::New();
    //curvatureTerm->SetInput(input);
    curvatureTerm->SetCoefficient(atof(argv[4]));


    // Create term container (equation rhs)
    typedef itk::LevelSetEquationTermContainer< FIT, LSContainerT >
TermContainerType;
    TermContainerType::Pointer termContainer = TermContainerType::New();
    termContainer->SetLevelSetContainer( levelSetContainer );
    termContainer->SetInput( input );
    //termContainer->AddTerm( 0, propagationTerm );
    //termContainer->AddTerm( 1, advectionTerm );
    termContainer->AddTerm( 0, curvatureTerm );

    // Create equation container
    typedef itk::LevelSetEquationContainer< TermContainerType >
EquationContainerType;
    EquationContainerType::Pointer equationContainer =
EquationContainerType::New();
    equationContainer->SetLevelSetContainer( levelSetContainer );
    equationContainer->AddEquation( 0, termContainer );

    // Create stopping criteria
    typedef itk::LevelSetEvolutionNumberOfIterationsStoppingCriterion<
LSContainerT > StoppingCriterionType;
    StoppingCriterionType::Pointer criterion = StoppingCriterionType::New();
    criterion->SetNumberOfIterations( numberOfIterations );

    // Create the visualizer
    typedef itk::VTKVisualizeImageLevelSetIsoValues< FIT, LST >
VisualizationType;
    VisualizationType::Pointer visualizer = VisualizationType::New();
    visualizer->SetInputImage( input );
    visualizer->SetLevelSet( levelSet );
    visualizer->SetScreenCapture( true );

    // Create evolution class
    typedef itk::LevelSetEvolution< EquationContainerType, LST >
LevelSetEvolutionType;
    LevelSetEvolutionType::Pointer evolution = LevelSetEvolutionType::New();
    evolution->SetEquationContainer( equationContainer );
    evolution->SetStoppingCriterion( criterion );
    evolution->SetLevelSetContainer( levelSetContainer );

    typedef itk::LevelSetIterationUpdateCommand< LevelSetEvolutionType,
VisualizationType > IterationUpdateCommandType;
    IterationUpdateCommandType::Pointer iterationUpdateCommand =
IterationUpdateCommandType::New();
    iterationUpdateCommand->SetFilterToUpdate( visualizer );
    iterationUpdateCommand->SetUpdatePeriod( 1 );

    evolution->AddObserver( itk::IterationEvent(), iterationUpdateCommand );

    evolution->Update();

    return EXIT_SUCCESS;
}



cmake_minimum_required(VERSION 2.8)

project(CurvatureMotion)

find_package(ITK REQUIRED)
include(${ITK_USE_FILE})

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

add_executable(Motion main.cpp)
target_link_libraries(Motion ${ITK_LIBRARIES}  ${VTK_LIBRARIES})

[image: Inline image 4][image: Inline image 3]
-- 
Cemal Cagatay Bilgin, PhD
Life Sciences Division
Lawrence Berkeley National Lab
MS977, One Cyclotron Road
Berkeley, CA 94720, USA
Email: ccbilgin at lbl.gov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140722/87e3963c/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: levelset_000.png
Type: image/png
Size: 2652 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/community/attachments/20140722/87e3963c/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: levelset_040.png
Type: image/png
Size: 3421 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/community/attachments/20140722/87e3963c/attachment-0003.png>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list