[ITK] [ITK-users] Different output with dense and sparse level set representations
Cagatay Bilgin
bilgincc at gmail.com
Thu Aug 21 13:33:53 EDT 2014
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140821/a2bfd8a6/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/20140821/a2bfd8a6/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/20140821/a2bfd8a6/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