<div dir="ltr">Hi Cagatay,<div><br></div><div>I just came back from vacation... Have you solved this problem?</div><div>It would be great if we could make the equivalent code in the old framework to compare and fix what could be wrong. Have you done something like this?<br>
</div><div><br></div><div>I'll start looking at it...</div><div><br></div><div>Best,</div><div>Arnaud</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Jul 22, 2014 at 9:58 AM, Cagatay Bilgin <span dir="ltr"><<a href="mailto:bilgincc@gmail.com" target="_blank">bilgincc@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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" target="_blank">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"><span class="HOEnZb"><font color="#888888"><br>
</font></span></div><span class="HOEnZb"><font color="#888888">
-- <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>
</font></span></div></div>
<br>_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" target="_blank">http://public.kitware.com/mailman/listinfo/insight-users</a><br>
<br></blockquote></div><br></div>