<div dir="ltr"><div><div>Hello everybody, </div><div><br></div><div>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. </div>
<div><br></div><div><img src="cid:ii_147f9895b639cfbd" alt="Inline image 2" width="300" height="300"><img src="cid:ii_147f9890c7a614f6" alt="Inline image 1" width="300" height="300"></div></div><div><br></div><div><br></div>
<div>#include "itkBinaryImageToLevelSetImageAdaptor.h"</div><div>#include "itkLevelSetContainer.h"</div><div>#include "itkLevelSetEquationPropagationTerm.h"</div><div>#include "itkLevelSetEquationAdvectionTerm.h"</div>
<div>#include "itkLevelSetEquationContainer.h"</div><div>#include "itkLevelSetEquationTermContainer.h"</div><div>#include "itkLevelSetEvolution.h"</div><div>#include "itkLevelSetEvolutionNumberOfIterationsStoppingCriterion.h"</div>
<div>#include "itkLevelSetDenseImage.h"</div><div>#include "itkVTKVisualizeImageLevelSetIsoValues.h"</div><div>#include "itkSinRegularizedHeavisideStepFunction.h"</div><div>#include "itkLevelSetIterationUpdateCommand.h"</div>
<div>#include "itkLevelSetEquationCurvatureTerm.h"</div><div>#include "itkWhitakerSparseLevelSetImage.h"</div><div>#include "itkLevelSetSparseImage.h"</div><div><br></div><div>typedef itk::Image< float, 2 >  FIT;</div>
<div><br></div><div>void CreateSquareImage(FIT::Pointer image, FIT::Pointer prop, FIT::Pointer advec)</div><div>{</div><div>    FIT::RegionType region;</div><div>    FIT::IndexType start;</div><div>    start[0] = 0;</div>
<div>    start[1] = 0;</div><div> </div><div>    FIT::SizeType size;</div><div>    size[0] = 100;</div><div>    size[1] = 100;</div><div> </div><div>    region.SetSize(size);</div><div>    region.SetIndex(start);</div><div>
 </div><div>    image->SetRegions(region);</div><div>    image->Allocate();</div><div>    image->FillBuffer(0);</div><div>    </div><div>    //constant grow in all directions</div><div>    prop->SetRegions(image->GetLargestPossibleRegion());</div>
<div>    prop->Allocate();</div><div>    prop->FillBuffer(1);</div><div><br></div><div>    //advec will increase in positive x direction</div><div>    advec->SetRegions(image->GetLargestPossibleRegion());</div>
<div>    advec->Allocate();</div><div>    advec->FillBuffer(0);</div><div><br></div><div>    // Set pixels in a square to one value</div><div>    for(unsigned int x = 35; x < 65; x++){</div><div>        for(unsigned int y = 35; y < 65; y++){</div>
<div>            FIT::IndexType pixelIndex;</div><div>            pixelIndex[0] = x;</div><div>            pixelIndex[1] = y;</div><div> </div><div>            image->SetPixel(pixelIndex, 255);</div><div>        }</div>
<div>    }</div><div>    </div><div>    //advection in positive x direction</div><div>    for(unsigned int x = 0; x < 100; x++){</div><div>        for(unsigned int y = 0; y < 100; y++){</div><div>            FIT::IndexType pixelIndex;</div>
<div>            pixelIndex[0] = x;</div><div>            pixelIndex[1] = y;</div><div>  </div><div>            advec->SetPixel(pixelIndex, x);</div><div>        }</div><div>    }</div><div>}</div><div><br></div><div><br>
</div><div>/*</div><div> */</div><div>int main(int argc, char* argv[] )</div><div>{</div><div>    using namespace itk;</div><div>    if(argc != 5)</div><div>    {</div><div>        std::cerr << "Missing Arguments" << std::endl;</div>
<div>        std::cerr << argv[0] << std::endl;</div><div>        std::cerr << "1- Number of Iterations" << std::endl;</div><div>        std::cerr << "2- Propagation Term" << std::endl;</div>
<div>        std::cerr << "3- Advection Term" << std::endl;</div><div>        std::cerr << "4- Curvature Term" << std::endl;</div><div>        return EXIT_FAILURE;</div><div>    }</div>
<div><br></div><div>    FIT::Pointer input = FIT::New();</div><div>    FIT::Pointer prop  = FIT::New();</div><div>    FIT::Pointer advec  = FIT::New();</div><div>    </div><div>    CreateSquareImage(input, prop, advec);</div>
<div>    </div><div>    int numberOfIterations = atoi(argv[1]);</div><div><br></div><div>    typedef float LSPT;</div><div>    typedef Image< LSPT, 2 >    LSIT;</div><div>    typedef LevelSetDenseImage< LSIT >  LST;</div>
<div>    //typedef WhitakerSparseLevelSetImage< LSPT, 2 > LST;</div><div>        </div><div>    typedef LST::OutputRealType  LSOutputT;</div><div>    </div><div>    // convert a binary mask to a level-set function</div>
<div>    typedef BinaryImageToLevelSetImageAdaptor<FIT, LST >  BI2LST;</div><div>    BI2LST::Pointer adaptor = BI2LST::New();</div><div>    adaptor->SetInputImage(input );</div><div>    adaptor->Initialize();</div>
<div>    LST::Pointer levelSet = adaptor->GetLevelSet();</div><div><br></div><div>    // The Heaviside function</div><div>    typedef SinRegularizedHeavisideStepFunction<LSOutputT, LSOutputT> HeavisideFunctionType;</div>
<div>    HeavisideFunctionType::Pointer heaviside = HeavisideFunctionType::New();</div><div>    heaviside->SetEpsilon(1 );</div><div><br></div><div>    // Create the level set container</div><div>    typedef LevelSetContainer< IdentifierType, LST > LSContainerT;</div>
<div>    LSContainerT::Pointer levelSetContainer = LSContainerT::New();</div><div>    levelSetContainer->SetHeaviside(heaviside );</div><div>    levelSetContainer->AddLevelSet(0, levelSet );</div><div><br></div><div>
    // Create the terms.</div><div>    typedef LevelSetEquationPropagationTerm<FIT, LSContainerT>  PropT;</div><div>    PropT::Pointer propagationTerm = PropT::New();</div><div>    propagationTerm->SetInput(prop);</div>
<div>    propagationTerm->SetCoefficient(atof(argv[2]));</div><div>  </div><div>    typedef LevelSetEquationAdvectionTerm<FIT, LSContainerT>  AdvecT;</div><div>    AdvecT::Pointer advectionTerm = AdvecT::New();</div>
<div>    advectionTerm->SetInput(advec);</div><div>    advectionTerm->SetCoefficient(atof(argv[3]));</div><div>  </div><div>    typedef LevelSetEquationCurvatureTerm<FIT, LSContainerT>  CurvT;</div><div>    CurvT::Pointer curvatureTerm = CurvT::New();</div>
<div>    curvatureTerm->SetCoefficient(atof(argv[4]));</div><div>    </div><div>    // Create term container (equation rhs)</div><div>    typedef LevelSetEquationTermContainer< FIT, LSContainerT > TermContainerT;</div>
<div>    TermContainerT::Pointer termContainer = TermContainerT::New();</div><div>    termContainer->SetLevelSetContainer(levelSetContainer );</div><div>    termContainer->SetInput(input );</div><div>    termContainer->AddTerm(0, curvatureTerm );</div>
<div>    termContainer->AddTerm(1, advectionTerm );</div><div>    termContainer->AddTerm(2, propagationTerm );</div><div>    </div><div>    // Create equation container</div><div>    typedef LevelSetEquationContainer< TermContainerT > EquationContainerType;</div>
<div>    EquationContainerType::Pointer equationContainer = EquationContainerType::New();</div><div>    equationContainer->SetLevelSetContainer(levelSetContainer );</div><div>    equationContainer->AddEquation(0, termContainer );</div>
<div><br></div><div>    // Create stopping criteria</div><div>    typedef LevelSetEvolutionNumberOfIterationsStoppingCriterion< LSContainerT > StoppingCriterionType;</div><div>    StoppingCriterionType::Pointer criterion = StoppingCriterionType::New();</div>
<div>    criterion->SetNumberOfIterations(numberOfIterations );</div><div><br></div><div>    // Create the visualizer</div><div>    typedef VTKVisualizeImageLevelSetIsoValues< FIT, LST > VisT;</div><div>    VisT::Pointer visualizer = VisT::New();</div>
<div>    visualizer->SetInputImage(input );</div><div>    visualizer->SetLevelSet(levelSet );</div><div>    visualizer->SetScreenCapture(true );</div><div><br></div><div>    // Create evolution class</div><div>    typedef LevelSetEvolution< EquationContainerType, LST > LSEvolT;</div>
<div>    LSEvolT::Pointer evolution = LSEvolT::New();</div><div>    evolution->SetEquationContainer(equationContainer );</div><div>    evolution->SetStoppingCriterion(criterion );</div><div>    evolution->SetLevelSetContainer(levelSetContainer );</div>
<div><br></div><div>    typedef LevelSetIterationUpdateCommand< LSEvolT, VisT > IterationUpdateCommandType;</div><div>    IterationUpdateCommandType::Pointer iterationUpdateCommand = IterationUpdateCommandType::New();</div>
<div>    iterationUpdateCommand->SetFilterToUpdate(visualizer );</div><div>    iterationUpdateCommand->SetUpdatePeriod(1 );</div><div><br></div><div>    evolution->AddObserver(IterationEvent(), iterationUpdateCommand );</div>
<div>    evolution->Update();</div><div><br></div><div>    return EXIT_SUCCESS;</div><div>}</div><div><br></div><div><br></div><div><div>cmake_minimum_required(VERSION 2.8)</div><div>project(Motion)</div><div><br></div>
<div>find_package(ITK REQUIRED)</div><div>include(${ITK_USE_FILE})</div><div><br></div><div>find_package(VTK REQUIRED)</div><div>include(${VTK_USE_FILE})</div><div><br></div><div>add_executable(Motion main.cpp)</div><div>
target_link_libraries(Motion ${ITK_LIBRARIES}  ${VTK_LIBRARIES})</div></div><div><br></div><div>Params used for the experiment:</div><div>Ran sparse representation with </div><div><div>./Motion 15 0 1 1</div></div><div><br>
</div><div><div>Ran dense representation with </div><div>./Motion 15 0 1 0</div></div><div><br></div><div> <br></div><span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">Cemal Cagatay Bilgin</span>, <span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">PhD</span><br style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">
<span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">Life Sciences Division</span><br style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">
<span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">Lawrence Berkeley National Lab</span><br style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">
<span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">MS977, One Cyclotron Road</span><br style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">
<span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">Berkeley, CA 94720, USA</span><br style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">
<span style="color:rgb(136,136,136);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)">Email: <a href="mailto:ccbilgin@lbl.gov" target="_blank">ccbilgin@lbl.gov</a></span>
</div>