[Insight-users] exception with GeodesicActiveContours 3D segmentation

[Marco Gheza]

Hi,
I'm trying to use GeodesicActiveContoursSegmentation with a 3D volume. I
want to use a PET gradient image as the feature image and the CT image as
the image to be segmented.
This is my code:

#include "itkImage.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"

#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"


int main( int argc, char *argv[] )
{
  if( argc < 7 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImageCT  inputImagePT outputImage";
    std::cerr << " seedX seedY seedZ";
    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();

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

  thresholder->SetOutsideValue(  0  );
  thresholder->SetInsideValue(  255 );


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

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

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


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

///filter in order to obtain the feature image///
  typedef   itk::CurvatureAnisotropicDiffusionImageFilter<
                               InternalImageType,
                               InternalImageType >  SmoothingFilterType;
  typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<
                               InternalImageType,
                               InternalImageType >  GradientFilterType;
  typedef   itk::SigmoidImageFilter<
                               InternalImageType,
                               InternalImageType >  SigmoidFilterType;

  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
  GradientFilterType::Pointer  gradientMagnitude =
GradientFilterType::New();
  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

  smoothing->SetTimeStep( 0.0625 );
  smoothing->SetNumberOfIterations(  5 );
  smoothing->SetConductanceParameter( 9.0 );

  const double sigma = 1.0;
  gradientMagnitude->SetSigma(  sigma  );

  const double alpha =  -0.5;
  const double beta  =  3.0;
  sigmoid->SetAlpha( alpha );
  sigmoid->SetBeta(  beta  );
  sigmoid->SetOutputMinimum(  0.0  );
  sigmoid->SetOutputMaximum(  1.0  );
///filter in order to obtain the feature image///


///level set///
  typedef  itk::FastMarchingImageFilter<
                              InternalImageType,
                              InternalImageType >    FastMarchingFilterType;

  typedef  itk::GeodesicActiveContourLevelSetImageFilter<
InternalImageType,
                InternalImageType >    GeodesicActiveContourFilterType;

  FastMarchingFilterType::Pointer  fastMarching =
FastMarchingFilterType::New();
  GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
                                     GeodesicActiveContourFilterType::New();

  const double propagationScaling = 2.0;

  geodesicActiveContour->SetPropagationScaling( propagationScaling );
  geodesicActiveContour->SetCurvatureScaling( 1.0 );
  geodesicActiveContour->SetAdvectionScaling( 1.0 );

  geodesicActiveContour->SetMaximumRMSError( 0.02 );
  geodesicActiveContour->SetNumberOfIterations( 800 );
///level set///


///utility for the fast marching///
  typedef FastMarchingFilterType::NodeContainer  NodeContainer;
  typedef FastMarchingFilterType::NodeType       NodeType;

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

  InternalImageType::IndexType  seedPosition;

  seedPosition[0] = atoi( argv[4] );
  seedPosition[1] = atoi( argv[5] );
  seedPosition[1] = atoi( argv[6] );

  const double initialDistance = 5.0;

  NodeType node;

  const double seedValue = - initialDistance;

  node.SetValue( seedValue );
  node.SetIndex( seedPosition );

  seeds->Initialize();
  seeds->InsertElement( 0, node );

  fastMarching->SetTrialPoints(  seeds  );

  fastMarching->SetSpeedConstant( 1.0 );
///utility for the fast marching///


  // 1° step
  smoothing->SetInput( readerPT->GetOutput() );
  gradientMagnitude->SetInput( smoothing->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );

  // 2° step
  geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
  geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );

  // 3° step
  thresholder->SetInput( geodesicActiveContour->GetOutput() );
  writer->SetInput( thresholder->GetOutput() );


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


unfortunatly, when i try to run it, i obtain this output message:

Exception caught !

itk::ExceptionObject (009EEB20)
Location: "__thiscall itk::ImageConstIterator<class itk::Image<class
itk::FixedArray<float,3>,3> >:
ImageConstIterator(const class itk::Image<class itk::FixedArray<float,3>,3>
*,const class itk::Imag
Region<3> &)"
File: c:\insighttoolkit-3.16.0\code\common\itkImageConstIterator.h
Line: 182
Description: itk::ERROR: Region ImageRegion (009EE8DC)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [128, 128, 287]
 is outside of buffered region ImageRegion (024035E0)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [16, 16, 16]


i don't know why. Please help me, thanks, bye,

Marco
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100503/3b823c73/attachment.htm>