[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