<div dir="ltr">Dear ITK Community, <div><br></div><div>I am trying to familiarize myself with the new level set classes. My goal is to simulate motion under mean curvature using the new design. </div><div><br></div><div>I have put together the following, looking at the examples and tests. <a href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3872740/" style="white-space:pre-wrap">http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3872740/</a> is a great source to understand the design and tests in source directory are very helpful, thank you for the resources! </div>
<div><br></div><div>The test scenario I have is a L shaped object. I would expect the object to become somewhat ellipse and then disappear at the end of the evolution. However I get the following results attached to the email. I don't quite follow these results. Am I missing something obvious? </div>
<div><br></div><div>Here is the code and CMakeLists.txt. I ran the program with ./Motion 100 0 0 0.05</div><div><br></div><div><div>#include "itkBinaryImageToLevelSetImageAdaptor.h"</div><div>#include "itkImageFileReader.h"</div>
<div>#include "itkLevelSetContainer.h"</div><div>#include "itkLevelSetEquationPropagationTerm.h"</div><div>#include "itkLevelSetEquationAdvectionTerm2.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 "itkCastImageFilter.h"</div><div>#include "itkWhitakerSparseLevelSetImage.h"</div><div>#include "itkSpatialObjectToImageFilter.h"</div><div>#include "itkEllipseSpatialObject.h"</div>
<div><br></div><div>typedef itk::Image< float, 2 > FIT;</div><div><br></div><div>/*</div><div> * L Shape </div><div> */</div><div>void CreateSquareImage(FIT::Pointer image)</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] = 200;</div><div> size[1] = 300;</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> </div><div> // Set pixels in a square to one value</div><div> for(unsigned int r = 20; r < 160; r++)</div>
<div> {</div><div> for(unsigned int c = 30; c < 100; c++)</div><div> {</div><div> FIT::IndexType pixelIndex;</div><div> pixelIndex[0] = r;</div><div> pixelIndex[1] = c;</div>
<div> </div><div> image->SetPixel(pixelIndex, 255);</div><div> }</div><div> }</div><div> </div><div> for(unsigned int r = 20; r < 80; r++)</div><div> {</div><div> for(unsigned int c = 100; c < 200; c++)</div>
<div> {</div><div> FIT::IndexType pixelIndex;</div><div> pixelIndex[0] = r;</div><div> pixelIndex[1] = c;</div><div> </div><div> image->SetPixel(pixelIndex, 255);</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> 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> // Image Dimension</div><div> const unsigned int Dimension = 2;</div><div><br></div><div> typedef unsigned char InputPixelType;</div><div> typedef itk::Image< InputPixelType, Dimension > IIT;</div>
<div> typedef itk::Image< float, 2 > FIT;</div><div><br></div><div> FIT::Pointer input = FIT::New();</div><div> CreateSquareImage(input);</div><div> </div><div> int numberOfIterations = atoi( argv[1]);</div>
<div><br></div><div> typedef float LevelSetPixelType;</div><div> typedef itk::Image< LevelSetPixelType, Dimension > LSIT;</div><div> typedef itk::LevelSetDenseImage< LSIT > LST;</div>
<div> //typedef itk::WhitakerSparseLevelSetImage< LevelSetPixelType, 2 > LST;</div><div><br></div><div> typedef LST::OutputType LevelSetOutputType;</div><div> typedef LST::OutputRealType LevelSetRealType;</div>
<div><br></div><div> // convert a binary mask to a level-set function</div><div> typedef itk::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 itk::SinRegularizedHeavisideStepFunction< LevelSetRealType, LevelSetRealType > 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 itk::LevelSetContainer< itk::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 itk::LevelSetEquationPropagationTerm<FIT, LSContainerT> PropagationTermType;</div><div> PropagationTermType::Pointer propagationTerm = PropagationTermType::New();</div>
<div> propagationTerm->SetInput(input);</div><div> propagationTerm->SetCoefficient(atof(argv[2]));</div><div> </div><div> typedef itk::LevelSetEquationAdvectionTerm2<FIT, LSContainerT> AdvectionTermType;</div>
<div> AdvectionTermType::Pointer advectionTerm = AdvectionTermType::New();</div><div> advectionTerm->SetInput(input);</div><div> advectionTerm->SetCoefficient(atof(argv[3]));</div><div> </div><div> typedef itk::LevelSetEquationCurvatureTerm<FIT, LSContainerT> CurvatureTermType;</div>
<div> CurvatureTermType::Pointer curvatureTerm = CurvatureTermType::New();</div><div> //curvatureTerm->SetInput(input);</div><div> curvatureTerm->SetCoefficient(atof(argv[4]));</div><div> </div><div><br></div>
<div> // Create term container (equation rhs)</div><div> typedef itk::LevelSetEquationTermContainer< FIT, LSContainerT > TermContainerType;</div><div> TermContainerType::Pointer termContainer = TermContainerType::New();</div>
<div> termContainer->SetLevelSetContainer( levelSetContainer );</div><div> termContainer->SetInput( input );</div><div> //termContainer->AddTerm( 0, propagationTerm );</div><div> //termContainer->AddTerm( 1, advectionTerm );</div>
<div> termContainer->AddTerm( 0, curvatureTerm );</div><div> </div><div> // Create equation container</div><div> typedef itk::LevelSetEquationContainer< TermContainerType > 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 itk::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 itk::VTKVisualizeImageLevelSetIsoValues< FIT, LST > VisualizationType;</div>
<div> VisualizationType::Pointer visualizer = VisualizationType::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 itk::LevelSetEvolution< EquationContainerType, LST > LevelSetEvolutionType;</div><div> LevelSetEvolutionType::Pointer evolution = LevelSetEvolutionType::New();</div>
<div> evolution->SetEquationContainer( equationContainer );</div><div> evolution->SetStoppingCriterion( criterion );</div><div> evolution->SetLevelSetContainer( levelSetContainer );</div><div><br></div><div>
typedef itk::LevelSetIterationUpdateCommand< LevelSetEvolutionType, VisualizationType > 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( itk::IterationEvent(), iterationUpdateCommand );</div>
<div><br></div><div> evolution->Update();</div><div><br></div><div> return EXIT_SUCCESS;</div><div>}</div><div><br></div><div><br></div><div><br></div><div><div>cmake_minimum_required(VERSION 2.8)</div><div><br></div>
<div>project(CurvatureMotion)</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><img src="cid:ii_1475d0f98ca3f888" alt="Inline image 4" width="300" height="300"><img src="cid:ii_1475d0f97a98cb76" alt="Inline image 3" width="300" height="300"><br>
</div>
-- <br><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></div>