[Insight-users] segmentation problems using Fast marching image filter

vijay kumar desivij at yahoo.com
Mon Apr 30 18:13:07 EDT 2007


Hi,
  I modified FastMarchingImageFilter.cxx that is found in Examples\Segmentation to work on a 3D analyze format image. After compilation, I executed it from the command line:
   
  FM_Filter3D.exe mri4itk.hdr mri4itkResult.hdr 
  1.25 -8.0 75.0 10 15 
  108 125 1 110 125 5 112 124 10 113 124 15 113 123 20 107 130 25 94 136 32
   
  The second line has values for sigma, alpha, beta, threshold and stopping time
  The third line contains indices (x,y,z) for a set of 7 points. All three lines were concatenated into one when running the program. The segmentation results don't look good. The binary image contains a large white spot in each slice. I am trying to segment the right artery which appears as a white spot on the right side in the input image. The text of my program is given below. I tried to attach the image (mri4itk.hdr, mri4itk.img) to this email, but  the email bounced. The size is approximately 3 MB. Can I send the images somehow? 
   
  Thanks,
  Vijay
   
  #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 < 11 ) {
  std::cerr << "Missing Parameters " << std::endl;
  std::cerr << "Usage: " << argv[0];
  std::cerr << " inputImage outputImage";
  std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue seedX seedY seedZ" << 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::BinaryThresholdImageFilter< InternalImageType, 
  OutputImageType > ThresholdingFilterType;
  ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
  const InternalPixelType timeThreshold = atof( argv[6] );
  
  thresholder->SetLowerThreshold( 0.0 );
  thresholder->SetUpperThreshold( timeThreshold );
  thresholder->SetOutsideValue( 0 );
  thresholder->SetInsideValue( 255 );
  
  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 > CastFilterType;
  typedef itk::CurvatureAnisotropicDiffusionImageFilter< 
  InternalImageType, 
  InternalImageType > SmoothingFilterType;
  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
  
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< 
  InternalImageType, 
  InternalImageType > GradientFilterType;
  typedef itk::SigmoidImageFilter< 
  InternalImageType, 
  InternalImageType > SigmoidFilterType;
  
  GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
  sigmoid->SetOutputMinimum( 0.0 );
  sigmoid->SetOutputMaximum( 1.0 );
  
  typedef itk::FastMarchingImageFilter< InternalImageType, 
  InternalImageType > FastMarchingFilterType;
  FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
  
  smoothing->SetInput( reader->GetOutput() );
  gradientMagnitude->SetInput( smoothing->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );
  fastMarching->SetInput( sigmoid->GetOutput() );
  thresholder->SetInput( fastMarching->GetOutput() );
  writer->SetInput( thresholder->GetOutput() );
  smoothing->SetTimeStep( 0.0625 );
  smoothing->SetNumberOfIterations( 5 );
  smoothing->SetConductanceParameter( 9.0 );
  const double sigma = atof( argv[3] );
  gradientMagnitude->SetSigma( sigma );
  const double alpha = atof( argv[4] );
  const double beta = atof( argv[5] );
  sigmoid->SetAlpha( alpha );
  sigmoid->SetBeta( beta );
  
  typedef FastMarchingFilterType::NodeContainer NodeContainer;
  typedef FastMarchingFilterType::NodeType NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();
  
  InternalImageType::IndexType seedPosition1;
  InternalImageType::IndexType seedPosition2;
  InternalImageType::IndexType seedPosition3;
  InternalImageType::IndexType seedPosition4;
  InternalImageType::IndexType seedPosition5;
  InternalImageType::IndexType seedPosition6;
  InternalImageType::IndexType seedPosition7;
   
  seedPosition1[0] = atoi( argv[8] );
  seedPosition1[1] = atoi( argv[9] );
  seedPosition1[2] = atoi( argv[10] );
  seedPosition2[0] = atoi( argv[11] );
  seedPosition2[1] = atoi( argv[12] );
  seedPosition2[2] = atoi( argv[13] );
  seedPosition3[0] = atoi( argv[14] );
  seedPosition3[1] = atoi( argv[15] );
  seedPosition3[2] = atoi( argv[16] );
  seedPosition4[0] = atoi( argv[17] );
  seedPosition4[1] = atoi( argv[18] );
  seedPosition4[2] = atoi( argv[19] );
  seedPosition5[0] = atoi( argv[20] );
  seedPosition5[1] = atoi( argv[21] );
  seedPosition5[2] = atoi( argv[22] );
  seedPosition6[0] = atoi( argv[23] );
  seedPosition6[1] = atoi( argv[24] );
  seedPosition6[2] = atoi( argv[25] );
  seedPosition7[0] = atoi( argv[26] );
  seedPosition7[1] = atoi( argv[27] );
  seedPosition7[2] = atoi( argv[28] );
  NodeType seed1;
  NodeType seed2;
  NodeType seed3;
  NodeType seed4;
  NodeType seed5;
  NodeType seed6;
  NodeType seed7;
  const double seedValue = 0.0;
  seed1.SetValue( seedValue );
  seed1.SetIndex( seedPosition1 );
  seed2.SetValue( seedValue );
  seed2.SetIndex( seedPosition2 );
  seed3.SetValue( seedValue );
  seed3.SetIndex( seedPosition3 );
  seed3.SetValue( seedValue );
  seed3.SetIndex( seedPosition4 );
  seed3.SetValue( seedValue );
  seed3.SetIndex( seedPosition5 );
  seed3.SetValue( seedValue );
  seed3.SetIndex( seedPosition6 );
  seed3.SetValue( seedValue );
  seed3.SetIndex( seedPosition7 );
  
  seeds->Initialize();
  seeds->InsertElement( 0, seed1 );
  seeds->InsertElement( 1, seed2 );
  seeds->InsertElement( 2, seed3 );
  seeds->InsertElement( 3, seed4 );
  seeds->InsertElement( 4, seed5 );
  seeds->InsertElement( 5, seed4 );
  seeds->InsertElement( 6, seed5 );
  fastMarching->SetTrialPoints( seeds );
  
  CastFilterType::Pointer caster1 = CastFilterType::New();
  CastFilterType::Pointer caster2 = CastFilterType::New();
  CastFilterType::Pointer caster3 = CastFilterType::New();
  CastFilterType::Pointer caster4 = CastFilterType::New();
  WriterType::Pointer writer1 = WriterType::New();
  WriterType::Pointer writer2 = WriterType::New();
  WriterType::Pointer writer3 = WriterType::New();
  WriterType::Pointer writer4 = WriterType::New();
  caster1->SetInput( smoothing->GetOutput() );
  writer1->SetInput( caster1->GetOutput() );
  writer1->SetFileName("FastMarchingFilterOutput1.hdr");
  caster1->SetOutputMinimum( 0 );
  caster1->SetOutputMaximum( 255 );
  writer1->Update();
  caster2->SetInput( gradientMagnitude->GetOutput() );
  writer2->SetInput( caster2->GetOutput() );
  writer2->SetFileName("FastMarchingFilterOutput2.hdr");
  caster2->SetOutputMinimum( 0 );
  caster2->SetOutputMaximum( 255 );
  writer2->Update();
  caster3->SetInput( sigmoid->GetOutput() );
  writer3->SetInput( caster3->GetOutput() );
  writer3->SetFileName("FastMarchingFilterOutput3.hdr");
  caster3->SetOutputMinimum( 0 );
  caster3->SetOutputMaximum( 255 );
  writer3->Update();
  caster4->SetInput( fastMarching->GetOutput() );
  writer4->SetInput( caster4->GetOutput() );
  writer4->SetFileName("FastMarchingFilterOutput4.hdr");
  caster4->SetOutputMinimum( 0 );
  caster4->SetOutputMaximum( 255 );
  fastMarching->SetOutputSize( 
  reader->GetOutput()->GetBufferedRegion().GetSize() );
  const double stoppingTime = atof( argv[7] );
  fastMarching->SetStoppingValue( stoppingTime );
  
  try
  {
  writer->Update();
  }
  catch( itk::ExceptionObject & excep )
  {
  std::cerr << "Exception caught !" << std::endl;
  std::cerr << excep << std::endl;
  }
  writer4->Update();
  typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
  InternalWriterType::Pointer mapWriter = InternalWriterType::New();
  mapWriter->SetInput( fastMarching->GetOutput() );
  mapWriter->SetFileName("FastMarchingFilterOutput4.mha");
  mapWriter->Update();
  InternalWriterType::Pointer speedWriter = InternalWriterType::New();
  speedWriter->SetInput( sigmoid->GetOutput() );
  speedWriter->SetFileName("FastMarchingFilterOutput3.mha");
  speedWriter->Update();
  InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
  gradientWriter->SetInput( gradientMagnitude->GetOutput() );
  gradientWriter->SetFileName("FastMarchingFilterOutput2.mha");
  gradientWriter->Update();
  return 0;
  }

       
---------------------------------
Ahhh...imagining that irresistible "new car" smell?
 Check outnew cars at Yahoo! Autos.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070430/51e14404/attachment-0001.htm


More information about the Insight-users mailing list