[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/insight-users/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/insight-users/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/insight-users/attachments/20140915/fd3524b9/attachment-0003.png>
More information about the Insight-users
mailing list