[ITK-users] Motion under Mean Curvature with LevelSetv4 Framework
Cagatay Bilgin
bilgincc at gmail.com
Mon Sep 15 11:50:51 EDT 2014
Hi Arnaud,
Sorry for the late response, I somehow missed this email. I have not done
the same experiment in the old framework to compare against. I can giving
it a shot but the old framework is more mysterious to me.
Best,
Cagatay Bilgin
On Fri, Aug 8, 2014 at 1:44 AM, Arnaud Gelas <arnaudgelas at gmail.com> wrote:
> 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
>>
>>
>
--
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/20140915/18b53990/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/20140915/18b53990/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/20140915/18b53990/attachment-0001.png>
More information about the Insight-users
mailing list