[Insight-users] fast marching level set image filter

[ilker hacıhaliloğlu]

hi all i am doing image segmentation( bone segmentaiton) with fast marching
level sets i have changed a little bit the code which is given in the
examples. the data that i am working is mri data in analyze format. i have
changed it from dicom to analyze.

sigma=0.98 alfa=-0.64 beta=6 stoptime=32 threshold=40
but am geting black images no segmented regions.
i have used the insightapplication demos it worked well in that examples so
what am i doing wrong i am sending the code attached


#include "itkCurvatureAnisotropicDiffusionImageFilter.h"

#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"




#include "itkImage.h"
#include "itkFastMarchingImageFilter.h"

#include "itkBinaryThresholdImageFilter.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkRescaleIntensityImageFilter.h"


int main( int argc, char *argv[] )
{
if( argc < 21 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputImage outputImage seedX seedY seedZ seedX2 seedY2
seedZ2";
std::cerr << " seedX3 seedY3 seedZ3 seedX4 seedY4 seedZ4 seedX5 seedY5
seedZ5 ";
std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue"
<< std::endl;
return 1;
}



typedef float InternalPixelType;
const unsigned int Dimension = 3;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;



typedef itk::ImageFileReader< InternalImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;



ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();

reader->SetFileName( argv[1] );
writer->SetFileName( argv[2] );



typedef itk::RescaleIntensityImageFilter<
InternalImageType,
OutputImageType > RescaleFilterType;



// Discrete CurvatureAnistropicDiffusionImageFilter

typedef itk::CurvatureAnisotropicDiffusionImageFilter<
InternalImageType,
InternalImageType > SmoothingFilterType;

SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
smoothing->SetTimeStep( 0.0625 );
smoothing->SetNumberOfIterations( 5 );
smoothing-> SetConductanceParameter( 3 );
smoothing->SetInput( reader->GetOutput() );
smoothing->Update();



RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
rescaler->SetOutputMinimum( 0 );
rescaler->SetOutputMaximum( 255 );
rescaler->SetInput( smoothing->GetOutput() );




//GradientMagnitudeRecursiveGaussianImageFilter

typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
InternalImageType,
InternalImageType > GradientFilterType;
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();

const double sigma = atof( argv[18] );

gradientMagnitude->SetSigma( sigma );
gradientMagnitude->SetInput( smoothing->GetOutput() );
gradientMagnitude->Update();

//Sigmoid filter
typedef itk::SigmoidImageFilter<
InternalImageType,
InternalImageType > SigmoidFilterType;



SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();


sigmoid->SetOutputMinimum( 0.0 );
sigmoid->SetOutputMaximum( 1.0 );

const double alpha = atof( argv[19] );
const double beta = atof( argv[20] );
sigmoid->SetAlpha( alpha );
sigmoid->SetBeta( beta );
sigmoid->SetInput( gradientMagnitude->GetOutput() );
sigmoid->Update();

//Fast marching Filter

typedef itk::FastMarchingImageFilter< InternalImageType,
InternalImageType > FastMarchingFilterType;

FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
typedef FastMarchingFilterType::NodeContainer NodeContainer;
typedef FastMarchingFilterType::NodeType NodeType;

NodeContainer::Pointer seeds = NodeContainer::New();

InternalImageType::IndexType seedPosition;
seedPosition[0] = atoi( argv[3] );
seedPosition[1] = atoi( argv[4] );
seedPosition[2] = atoi( argv[5] );
InternalImageType::IndexType seedPosition2;
seedPosition2[0] = atoi( argv[6] );
seedPosition2[1] = atoi( argv[7] );
seedPosition2[2] = atoi( argv[8] );
InternalImageType::IndexType seedPosition3;
seedPosition3[0] = atoi( argv[9] );
seedPosition3[1] = atoi( argv[10] );
seedPosition3[2] = atoi( argv[11] );
InternalImageType::IndexType seedPosition4;
seedPosition4[0] = atoi( argv[12] );
seedPosition4[1] = atoi( argv[13] );
seedPosition4[2] = atoi( argv[14] );
InternalImageType::IndexType seedPosition5;
seedPosition5[0] = atoi( argv[15] );
seedPosition5[1] = atoi( argv[16] );
seedPosition5[2] = atoi( argv[17] );
// InternalImageType::IndexType seedPosition6;
// seedPosition6[0] = atoi( argv[18] );
//seedPosition6[1] = atoi( argv[19] );
//seedPosition6[2] = atoi( argv[20] );

NodeType node;
NodeType node2;
NodeType node3;
NodeType node4;
NodeType node5;
//NodeType node6;
const double seedValue = 0.0;

node.SetValue( seedValue );
node2.SetValue( seedValue );
node3.SetValue( seedValue );
node4.SetValue( seedValue );
node5.SetValue( seedValue );
//node6.SetValue( seedValue );

node.SetIndex( seedPosition );
node2.SetIndex( seedPosition2 );
node3.SetIndex( seedPosition3 );
node4.SetIndex( seedPosition4 );
node5.SetIndex( seedPosition5 );
//node6.SetIndex( seedPosition6 );


seeds->Initialize();
seeds->InsertElement( 0, node );
seeds->InsertElement( 0, node2 );
seeds->InsertElement( 0, node3 );
seeds->InsertElement( 0, node4 );
seeds->InsertElement( 0, node5 );
//seeds->InsertElement( 0, node6 );

fastMarching->SetTrialPoints( seeds );


fastMarching->SetOutputSize(
reader->GetOutput()->GetBufferedRegion().GetSize() );


const double stoppingTime = atof( argv[22] );
fastMarching->SetStoppingValue( stoppingTime );

fastMarching->SetInput( sigmoid->GetOutput() );
fastMarching->Update();

// Binary Threshold filter


typedef itk::BinaryThresholdImageFilter< InternalImageType,
OutputImageType > ThresholdingFilterType;
ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();


const InternalPixelType timeThreshold = atof( argv[21] );

thresholder->SetLowerThreshold( 0.0 );
thresholder->SetUpperThreshold( timeThreshold );

thresholder->SetOutsideValue( 0 );
thresholder->SetInsideValue( 1 );
thresholder->SetInput( fastMarching->GetOutput() );
writer->SetInput( thresholder->GetOutput() );


try
{
writer->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception caught !" << std::endl;
std::cerr << excep << std::endl;
}



return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051118/7ca10d3b/attachment-0001.html