[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/insight-users/attachments/20140722/87e3963c/attachment.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/insight-users/attachments/20140722/87e3963c/attachment.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/insight-users/attachments/20140722/87e3963c/attachment-0001.png>
More information about the Insight-users
mailing list