Hi all,<br><br>I would like to use the very usefull ITK classes to segment mammographic images with level set methods.<br>In order to learn these classes, I run properly the examples provided with the library and explained in the software guide.<br>
Therefore, I've tried to apply the filter ShapeDetectionLevelSetImageFilter to a very simple circular model,<br>but I cannot obtain any reasonable segmentation. <br><br>Hereafter, I report the code used (sorry for the long list of code lines!! :) ).<br>
The application accepts a list of arguments, precisely :<br>1_size_model 2_amplitude_model 3_sigma_smoothing 4_propagationscale 5 curvaturescale 6_levelset_initialvalue<br>I tried a lot of different argument lists, among them you can try the following one<br>
101 10 1 1 1 10<br><br>I would very appreciate to receive any comment, suggestion to solve the problem.<br><br>Thank u in advance,<br><br>Silvano<br><br><br>#include <itkIndex.h><br>#include <itkImage.h><br>#include <itkImageFileWriter.h><br>
#include <itkCastImageFilter.h><br>#include <itkTimeProbesCollectorBase.h><br>#include <itkGradientAnisotropicDiffusionImageFilter.h><br>#include <itkFastMarchingImageFilter.h><br>#include <itkShapeDetectionLevelSetImageFilter.h><br>
#include <itkGradientMagnitudeRecursiveGaussianImageFilter.h><br>#include <itkBoundedReciprocalImageFilter.h><br><br>typedef itk::Image<float,2> TImage;<br><br>TImage::Pointer model2D_flat( float sizeM, float amp )<br>
{<br> TImage::Pointer model = TImage::New();<br> <br> int hsize = (int)((sizeM-1.0)/2.0);<br> <br> std::cout << "size model image = " << sizeM;<br> std::cout << ", half size model image = " << hsize << std::endl;<br>
<br> TImage::IndexType start;<br> start[0] = 0; <br> start[1] = 0; <br><br> TImage::SizeType size;<br> size[0] = sizeM;<br> size[1] = sizeM;<br> <br> TImage::RegionType region;<br> region.SetSize( size );<br>
region.SetIndex( start );<br> <br> model->SetRegions( region );<br> model->Allocate();<br> <br> float radius2 = static_cast<float>(sizeM)/4.0;<br> radius2 = radius2 * radius2;<br> for ( int y=-hsize; y <= hsize; ++y )<br>
{<br> for ( int x=-hsize; x <= hsize; ++x )<br> {<br> TImage::IndexType index;<br> index[0] = (x+hsize);<br> index[1] = (y+hsize);<br> if ( (x*x+y*y) < radius2 )<br>
model->SetPixel( index, amp );<br> else<br> model->SetPixel( index, amp/2.0 );<br> }<br> }<br> <br> return model;<br>}<br><br><br>int main(int argc, char **argv)<br>
{<br> double sizeM = atof( argv[1] );<br> double amplitude = atof( argv[2] );<br> double sigma = atof( argv[3] );<br> double propscale = atof( argv[4] );<br> double curvscale = atof( argv[5] );<br> double isovalue = atof( argv[6] );<br>
<br> std::cout << "size= " << sizeM << std::endl;<br> std::cout << "amplitude= " << amplitude << std::endl;<br> std::cout << "sigma= " << sigma << std::endl;<br>
std::cout << "propscale= " << propscale << std::endl;<br> std::cout << "curvscale= " << curvscale << std::endl;<br> std::cout << "isovalue= " << isovalue << std::endl;<br>
<br> TImage::Pointer model = model2D_flat( sizeM, amplitude );<br> <br> itk::TimeProbesCollectorBase timeCollector;<br> timeCollector.Start("Total");<br> <br> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< TImage, TImage > GradientFilterType;<br>
GradientFilterType::Pointer gradient = GradientFilterType::New();<br> gradient->SetSigma( sigma );<br> gradient->SetNormalizeAcrossScale(false);<br> <br> typedef itk::BoundedReciprocalImageFilter< TImage, TImage > BoundFilterType;<br>
BoundFilterType::Pointer bounder = BoundFilterType::New();<br> <br> typedef itk::FastMarchingImageFilter< TImage, TImage > FastMarchingFilterType;<br> FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();<br>
<br> typedef itk::ShapeDetectionLevelSetImageFilter< TImage, TImage > ShapeDetectionFilterType;<br> ShapeDetectionFilterType::Pointer shapeDetection = ShapeDetectionFilterType::New();<br> <br> gradient->SetInput( model );<br>
bounder->SetInput( gradient->GetOutput() );<br> shapeDetection->SetInput( fastMarching->GetOutput() );<br> shapeDetection->SetFeatureImage( bounder->GetOutput() );<br> <br> typedef FastMarchingFilterType::NodeContainer NodeContainer;<br>
typedef FastMarchingFilterType::NodeType NodeType;<br> NodeContainer::Pointer seeds = NodeContainer::New();<br><br> NodeType node;<br> const double seedValue = -isovalue;<br> <br><br> itk::Index<2> pos; <br>
int hsize = (int)((sizeM-1.0)/2.0);<br> pos[0] = hsize; pos[1] = hsize;<br> node.SetValue( seedValue );<br> node.SetIndex( pos );<br><br> seeds->Initialize();<br> seeds->InsertElement( 0, node );<br>
<br> fastMarching->SetTrialPoints( seeds );<br> fastMarching->SetOutputSize( model->GetLargestPossibleRegion().GetSize() );<br> fastMarching->SetOutputSpacing( model->GetSpacing() );<br> fastMarching->SetOutputOrigin ( model->GetOrigin() );<br>
fastMarching->SetSpeedConstant( 1.0 );<br> <br> shapeDetection->SetPropagationScaling( propscale );<br> shapeDetection->SetCurvatureScaling( curvscale );<br> shapeDetection->SetMaximumRMSError( 0.02 );<br>
shapeDetection->SetNumberOfIterations( 1000 );<br><br> typedef itk::ImageFileWriter<TImage> WriterType;<br> WriterType::Pointer iwriter = WriterType::New();<br><br> timeCollector.Start("process");<br>
iwriter->SetInput( model );<br> iwriter->SetFileName( "model.nii" );<br> iwriter->Update();<br> <br> iwriter->SetInput( fastMarching->GetOutput() );<br> iwriter->SetFileName( "fastMarchOutput.nii" );<br>
iwriter->Update();<br> <br> iwriter->SetInput( shapeDetection->GetOutput() );<br> iwriter->SetFileName( "shapeOutput.nii" );<br> iwriter->Update();<br> <br> iwriter->SetInput( bounder->GetOutput() );<br>
iwriter->SetFileName( "edge.nii" );<br> iwriter->Update();<br> <br> timeCollector.Stop("process");<br> <br> std::cout << "No. elpased iterations: " << shapeDetection->GetElapsedIterations() << std::endl;<br>
std::cout << "RMS change: " << shapeDetection->GetRMSChange() << std::endl;<br><br> timeCollector.Stop("Total");<br> std::cout << std::endl;<br> std::cout << "Computing Time (sec): " << std::endl;<br>
timeCollector.Report();<br>}<br><br><br>