[ITK] [ITK-users] Motion under Mean Curvature with LevelSetv4 Framework
Arnaud Gelas
arnaudgelas at gmail.com
Fri Aug 8 04:44:26 EDT 2014
Hi Cagatay,
I just came back from vacation... Have you solved this problem?
It would be great if we could make the equivalent code in the old framework
to compare and fix what could be wrong. Have you done something like this?
I'll start looking at it...
Best,
Arnaud
On Tue, Jul 22, 2014 at 9:58 AM, Cagatay Bilgin <bilgincc at gmail.com> wrote:
> 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
>
> _____________________________________
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140808/d9e4d640/attachment-0001.html>
-------------- 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/20140808/d9e4d640/attachment-0002.png>
-------------- 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/20140808/d9e4d640/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