[ITK] [ITK-users] Different output with dense and sparse level set representations

Cagatay Bilgin bilgincc at gmail.com
Mon Sep 15 11:53:06 EDT 2014


Arnaud, should I also try this with the old framework? This really looks
like a bug to me regardless of the framework used.

Best,
Cagatay Bilgin

On Thu, Aug 21, 2014 at 10:33 AM, Cagatay Bilgin <bilgincc at gmail.com> wrote:

> Hello everybody,
>
> I continue with my exploration of the new level set framework. I have
> modified one of the examples and have an advection field in the positive x
> direction. I would expect the initial contour to move in that direction and
> that's the case with the sparse representation. Using the dense
> representation however, only half of the contour moves to my surprise. Is
> this a bug or I am missing something here. Here is the code and the cmake
> file and snapshots of the movements.
>
> [image: Inline image 2][image: Inline image 1]
>
>
> #include "itkBinaryImageToLevelSetImageAdaptor.h"
> #include "itkLevelSetContainer.h"
> #include "itkLevelSetEquationPropagationTerm.h"
> #include "itkLevelSetEquationAdvectionTerm.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 "itkWhitakerSparseLevelSetImage.h"
> #include "itkLevelSetSparseImage.h"
>
> typedef itk::Image< float, 2 >  FIT;
>
> void CreateSquareImage(FIT::Pointer image, FIT::Pointer prop, FIT::Pointer
> advec)
> {
>     FIT::RegionType region;
>     FIT::IndexType start;
>     start[0] = 0;
>     start[1] = 0;
>
>     FIT::SizeType size;
>     size[0] = 100;
>     size[1] = 100;
>
>     region.SetSize(size);
>     region.SetIndex(start);
>
>     image->SetRegions(region);
>     image->Allocate();
>     image->FillBuffer(0);
>
>     //constant grow in all directions
>     prop->SetRegions(image->GetLargestPossibleRegion());
>     prop->Allocate();
>     prop->FillBuffer(1);
>
>     //advec will increase in positive x direction
>     advec->SetRegions(image->GetLargestPossibleRegion());
>     advec->Allocate();
>     advec->FillBuffer(0);
>
>     // Set pixels in a square to one value
>     for(unsigned int x = 35; x < 65; x++){
>         for(unsigned int y = 35; y < 65; y++){
>             FIT::IndexType pixelIndex;
>             pixelIndex[0] = x;
>             pixelIndex[1] = y;
>
>             image->SetPixel(pixelIndex, 255);
>         }
>     }
>
>     //advection in positive x direction
>     for(unsigned int x = 0; x < 100; x++){
>         for(unsigned int y = 0; y < 100; y++){
>             FIT::IndexType pixelIndex;
>             pixelIndex[0] = x;
>             pixelIndex[1] = y;
>
>             advec->SetPixel(pixelIndex, x);
>         }
>     }
> }
>
>
> /*
>  */
> int main(int argc, char* argv[] )
> {
>     using namespace itk;
>     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;
>     }
>
>     FIT::Pointer input = FIT::New();
>     FIT::Pointer prop  = FIT::New();
>     FIT::Pointer advec  = FIT::New();
>
>     CreateSquareImage(input, prop, advec);
>
>     int numberOfIterations = atoi(argv[1]);
>
>     typedef float LSPT;
>     typedef Image< LSPT, 2 >    LSIT;
>     typedef LevelSetDenseImage< LSIT >  LST;
>     //typedef WhitakerSparseLevelSetImage< LSPT, 2 > LST;
>
>     typedef LST::OutputRealType  LSOutputT;
>
>     // convert a binary mask to a level-set function
>     typedef BinaryImageToLevelSetImageAdaptor<FIT, LST >  BI2LST;
>     BI2LST::Pointer adaptor = BI2LST::New();
>     adaptor->SetInputImage(input );
>     adaptor->Initialize();
>     LST::Pointer levelSet = adaptor->GetLevelSet();
>
>     // The Heaviside function
>     typedef SinRegularizedHeavisideStepFunction<LSOutputT, LSOutputT>
> HeavisideFunctionType;
>     HeavisideFunctionType::Pointer heaviside =
> HeavisideFunctionType::New();
>     heaviside->SetEpsilon(1 );
>
>     // Create the level set container
>     typedef LevelSetContainer< IdentifierType, LST > LSContainerT;
>     LSContainerT::Pointer levelSetContainer = LSContainerT::New();
>     levelSetContainer->SetHeaviside(heaviside );
>     levelSetContainer->AddLevelSet(0, levelSet );
>
>     // Create the terms.
>     typedef LevelSetEquationPropagationTerm<FIT, LSContainerT>  PropT;
>     PropT::Pointer propagationTerm = PropT::New();
>     propagationTerm->SetInput(prop);
>     propagationTerm->SetCoefficient(atof(argv[2]));
>
>     typedef LevelSetEquationAdvectionTerm<FIT, LSContainerT>  AdvecT;
>     AdvecT::Pointer advectionTerm = AdvecT::New();
>     advectionTerm->SetInput(advec);
>     advectionTerm->SetCoefficient(atof(argv[3]));
>
>     typedef LevelSetEquationCurvatureTerm<FIT, LSContainerT>  CurvT;
>     CurvT::Pointer curvatureTerm = CurvT::New();
>     curvatureTerm->SetCoefficient(atof(argv[4]));
>
>     // Create term container (equation rhs)
>     typedef LevelSetEquationTermContainer< FIT, LSContainerT >
> TermContainerT;
>     TermContainerT::Pointer termContainer = TermContainerT::New();
>     termContainer->SetLevelSetContainer(levelSetContainer );
>     termContainer->SetInput(input );
>     termContainer->AddTerm(0, curvatureTerm );
>     termContainer->AddTerm(1, advectionTerm );
>     termContainer->AddTerm(2, propagationTerm );
>
>     // Create equation container
>     typedef LevelSetEquationContainer< TermContainerT >
> EquationContainerType;
>     EquationContainerType::Pointer equationContainer =
> EquationContainerType::New();
>     equationContainer->SetLevelSetContainer(levelSetContainer );
>     equationContainer->AddEquation(0, termContainer );
>
>     // Create stopping criteria
>     typedef LevelSetEvolutionNumberOfIterationsStoppingCriterion<
> LSContainerT > StoppingCriterionType;
>     StoppingCriterionType::Pointer criterion =
> StoppingCriterionType::New();
>     criterion->SetNumberOfIterations(numberOfIterations );
>
>     // Create the visualizer
>     typedef VTKVisualizeImageLevelSetIsoValues< FIT, LST > VisT;
>     VisT::Pointer visualizer = VisT::New();
>     visualizer->SetInputImage(input );
>     visualizer->SetLevelSet(levelSet );
>     visualizer->SetScreenCapture(true );
>
>     // Create evolution class
>     typedef LevelSetEvolution< EquationContainerType, LST > LSEvolT;
>     LSEvolT::Pointer evolution = LSEvolT::New();
>     evolution->SetEquationContainer(equationContainer );
>     evolution->SetStoppingCriterion(criterion );
>     evolution->SetLevelSetContainer(levelSetContainer );
>
>     typedef LevelSetIterationUpdateCommand< LSEvolT, VisT >
> IterationUpdateCommandType;
>     IterationUpdateCommandType::Pointer iterationUpdateCommand =
> IterationUpdateCommandType::New();
>     iterationUpdateCommand->SetFilterToUpdate(visualizer );
>     iterationUpdateCommand->SetUpdatePeriod(1 );
>
>     evolution->AddObserver(IterationEvent(), iterationUpdateCommand );
>     evolution->Update();
>
>     return EXIT_SUCCESS;
> }
>
>
> cmake_minimum_required(VERSION 2.8)
> project(Motion)
>
> 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})
>
> Params used for the experiment:
> Ran sparse representation with
> ./Motion 15 0 1 1
>
> Ran dense representation with
> ./Motion 15 0 1 0
>
>
> Cemal Cagatay Bilgin, PhD
> Life Sciences Division
> Lawrence Berkeley National Lab
> MS977, One Cyclotron Road
> Berkeley, CA 94720, USA
> Email: ccbilgin at lbl.gov
>



-- 
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/20140915/fd3524b9/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sparse_00015.png
Type: image/png
Size: 2848 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/community/attachments/20140915/fd3524b9/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dense_00010.png
Type: image/png
Size: 2737 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/community/attachments/20140915/fd3524b9/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