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

[Cagatay Bilgin]

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/insight-users/attachments/20140821/a2bfd8a6/attachment.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/20140821/a2bfd8a6/attachment.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/20140821/a2bfd8a6/attachment-0001.png>