[Insight-users] 3D signed Distance Map using fast marching image filter

[Prathyusha Kamarajugadda]

Hello,

    I need as 3D distance map (signed) to initialize a level set for
Geodesic Active Contour Segmentation. I want to implement using fast
marching filter (as it says in the itk software tutorial). My code is
working fine for 2D distance map but when I increased the dimension to 3,
it is showing a very absurd output. It is showing intensity values of
2.59619e+033 (inside - black region) and 1.70141e+038 (outside - white
region).
I attached my code and a snapshot of my 3D output.

Please let me know if anyone has done the same in 3D.


Prathyusha Kamarajugadda,
Junior Undergraduate,
Bachelor of Technology,
Electrical Engineering Department,
National Institute of Technology, Warangal.


Prathyusha Kamarajugadda,
Junior Undergraduate,
Bachelor of Technology,
Electrical Engineering Department,
National Institute of Technology, Warangal.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130708/7787b7be/attachment-0001.htm>
-------------- next part --------------
#include "itkFastMarchingImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkImage.h"
static void PrintCommandLineUsage( const int argc, const char * const argv[] )
{
  std::cerr << "Missing Parameters " << std::endl;
  std::cerr << "Usage: " << argv[0];
  std::cerr << " outputImage seedX seedY seedZ initialDistance";
  std::cerr << "TimeThreshold StoppingValue" << std::endl;

  for(int qq=0; qq< argc; ++qq)
    {
    std::cout << "argv[" << qq << "] = " << argv[qq] << std::endl;
    }
  return;
}

int main( int argc, char *argv[] )
{
  if( argc != 8 )
    {
    PrintCommandLineUsage(argc, argv);
    return EXIT_FAILURE;
    }

  std::cout << argv[2] <<" " << argv[3] << " " << argv[4] <<std::endl;


  typedef   float           InternalPixelType; 
  const     unsigned int    Dimension = 3;
  typedef itk::Image< InternalPixelType, Dimension> ImageType;
  ImageType::Pointer image = ImageType::New();
  ImageType::SizeType  size;

  size[0]  = 64;  
  size[1]  = 64;  
  size[2]  = 32;

  ImageType::IndexType start;
  start[0] = 0; // first index on X
  start[1] = 0; // first index on Y
  start[2] = 0; // first index on Z

  ImageType::RegionType region;
  region.SetSize( size );
  region.SetIndex( start );
  image->SetRegions( region );
  image->Allocate();

  //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< ImageType,
                        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< ImageType > ReaderType;
  typedef  itk::ImageFileWriter<  ImageType  > WriterType;

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

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

  typedef itk::RescaleIntensityImageFilter<
                               ImageType,
                               OutputImageType >   CastFilterType;
  typedef  itk::FastMarchingImageFilter< ImageType,
                              ImageType >    FastMarchingFilterType;
  FastMarchingFilterType::Pointer  fastMarching
                                              = FastMarchingFilterType::New();
  typedef FastMarchingFilterType::NodeContainer NodeContainer;
  typedef FastMarchingFilterType::NodeType NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();

  ImageType::IndexType  seedPosition;

  seedPosition[0] = atoi( argv[2] );
  seedPosition[1] = atoi( argv[3] );
  seedPosition[2] = atoi( argv[4] );
  const float initialdistance = atof( argv[5] );

  NodeType node;
  const double seedValue = - initialdistance;
  node.SetValue( seedValue );
  node.SetIndex( seedPosition );

  std::cout << node.GetIndex() << std::endl;

  seeds->Initialize();
  seeds->InsertElement( 0, node );
  fastMarching->SetTrialPoints( seeds );
  fastMarching->SetSpeedConstant( 1.0 );

  fastMarching->SetOutputSize( size );
  //fastMarching->SetInput( image );
  const double stoppingTime = atof( argv[7] );
  fastMarching->SetStoppingValue(  stoppingTime  );

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


 try
    {
    CastFilterType::Pointer caster4 = CastFilterType::New();
    WriterType::Pointer writer4 = WriterType::New();
    caster4->SetInput( fastMarching->GetOutput() );
    //writer4->SetInput( caster4->GetOutput() );
	writer4->SetInput( fastMarching->GetOutput() );
    writer4->SetFileName("FastMarchingFilterOutput4.nii");
    caster4->SetOutputMinimum(   0 );
    caster4->SetOutputMaximum( 255 );
    writer4->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cout << "ExceptionObject caught !" << std::endl;
    std::cout << err << std::endl;
    return EXIT_FAILURE;
    }

  return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3D distance map output.jpg
Type: image/jpeg
Size: 189535 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130708/7787b7be/attachment-0001.jpg>