[Insight-users] problems using geodesic active contours segmentation

[Paul]

Hello mailinglist,

I've got some problems using the Geodesic Active Contours Segmentation 
that is explained in the ITK Software Guide. I hope you can help me. In 
Figure 9.21 you can see the pipeline. Creating the Edge Image works 
fine. The problems start with creating the Input Levelset. First I'm not 
sure, were I have to set the seedpoints, in the object, on the contour, 
outside? I'm working with ct images of the spine (DICOM, 3D). Here is a 
part of the code. And a Link to the distance map und the segmentation 
result. I'm using MITK, a toolkit with interaction tools, like setting 
seed points. Region Growing und Fast Marching Segmentation works fine 
with it.

IMG1: http://250kb.de/f8Wg0ug  , shows edge image (black and white) and 
the distancemap ( red points)

IMG2: http://250kb.de/NZ8lIUz  , shows edge image (black and white) and 
the segmentation result ( red points)

(here I set the seed points along the horizontal and the verticle axis 
of the shown image)


CODE:

typedef float        InternalPixelType;
const   unsigned int Dimension = 3;

typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

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

typedef unsigned char CharPixelType;
typedef short ShortPixelType;

typedef itk::Image<CharPixelType, Dimension> CharImageType;
typedef itk::Image<ShortPixelType, Dimension> ShortImageType;

   typedef itk::RescaleIntensityImageFilter<itk::Image<TPixel, 
VImageDimension>, InternalImageType > FilterType;

   typename FilterType::Pointer filter = FilterType::New();
   filter->SetOutputMinimum(   -961.00 ); //  MINIMUM PIXEL VALUE OF MY 
DATA SET
   filter->SetOutputMaximum( 3106.00 ); //    MAXIMUM PIXEL VALUE OF MY 
DATA SET

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

   thresholder-> SetLowerThreshold(getLower()); // I'M USING THE VALUES 
FROM THE GUIDE, DEFAULT: -1000   (FROM SHAPE DETECTION SEGMENTATION)
   thresholder-> SetUpperThreshold(getUpper()); // I'M USING THE VALUES 
FROM THE GUIDE, DEFAULT: 0       (FROM SHAPE DETECTION SEGMENTATION)

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

   typedef 
itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType, 
InternalImageType > SmoothingFilterType;
   typename SmoothingFilterType::Pointer smoothing = 
SmoothingFilterType::New();

   typedef 
itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType, 
InternalImageType > GradientFilterType;
   typedef itk::SigmoidImageFilter<InternalImageType, InternalImageType 
 >  SigmoidFilterType;

   typename GradientFilterType::Pointer gradientMagnitude = 
GradientFilterType::New();
   typename SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

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

   typedef itk::FastMarchingImageFilter< InternalImageType, 
InternalImageType > FastMarchingFilterType;
   typename FastMarchingFilterType::Pointer fastMarching = 
FastMarchingFilterType::New();

   typedef  itk::ShapeDetectionLevelSetImageFilter< InternalImageType, 
InternalImageType > ShapeDetectionFilterType;
   typename ShapeDetectionFilterType::Pointer shapeDetection = 
ShapeDetectionFilterType::New();

   typedef  itk::GeodesicActiveContourLevelSetImageFilter< 
InternalImageType, InternalImageType >    GeodesicActiveContourFilterType;
   typename GeodesicActiveContourFilterType::Pointer 
geodesicActiveContour = GeodesicActiveContourFilterType::New();

   filter->SetInput( itkImage ); // RESCALING TO FLOAT AS INPUT FOR THE 
CurvatureAnisotropicDiffusionImageFilter
   smoothing->SetInput( filter->GetOutput() );
   gradientMagnitude->SetInput( smoothing->GetOutput() );
   sigmoid->SetInput( gradientMagnitude->GetOutput() );
   fastMarching->SetInput( filter->GetOutput() ); // MITK DOES NOT 
ACCEPT IMAGES WITH A DIMENSION OF ZERO, FOR THIS REASON I HAVE TO SET AN 
INPUT FOR THE FAST MARCHING
   geodesicActiveContour->SetInput( fastMarching->GetOutput() );
   geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
   thresholder->SetInput( geodesicActiveContour->GetOutput() );


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

   gradientMagnitude->SetSigma( getSigma() ); // I'M USING THE VALUES 
FROM THE GUIDE, DEFAULT: 1.0

   sigmoid->SetAlpha( getAlpha() ); // I'M USING THE VALUES FROM THE 
GUIDE, DEFAULT: -0.5
   sigmoid->SetBeta( getBeta() );   // I'M USING THE VALUES FROM THE 
GUIDE, DEFAULT: 2.0



   typedef typename FastMarchingFilterType::NodeContainer NodeContainer;
   typedef typename FastMarchingFilterType::NodeType  NodeType;
   typename NodeContainer::Pointer seeds = NodeContainer::New();

   itk::Index<3> seedIndex;

   NodeType node;
   double seedValue= -5.0;  // - initialDistance

   seeds->Initialize();
   int i=0;

   mitk::PointSet::PointsContainer* points = 
m_PointSet->GetPointSet()->GetPoints();
   for ( mitk::PointSet::PointsConstIterator pointsIterator = 
points->Begin();
         pointsIterator != points->End();
         ++pointsIterator )
   {
   if ( !imageGeometry->IsInside( pointsIterator.Value()) )
     {
       continue;
     }

   // convert world coordinates to image indices
     imageGeometry->WorldToIndex( pointsIterator.Value(), seedIndex);

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

   seeds->InsertElement( i, node );
   i+=1;
   }


   fastMarching->SetTrialPoints( seeds );

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

   fastMarching->SetSpeedConstant( 1.0 );

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


try
   {

   thresholder->Update();

// PUTTING THE RESULT IMAGE AS NODE IN THE DATATREE OF MITK
   mitk::Image::Pointer resultImage= mitk::ImportItkImage (thresholder-> 
GetOutput());
   mitk::DataTreeNode::Pointer newNode= mitk::DataTreeNode::New();
   newNode->SetData( resultImage );

   newNode->SetProperty("binary", mitk::BoolProperty::New(true));
   newNode->SetProperty("name", mitk::StringProperty::New("dumb 
segmentation"));
   newNode->SetProperty("color", mitk::ColorProperty::New(1.0,0.0,0.0));
   newNode->SetProperty("volumerendering", mitk::BoolProperty::New(true));
   newNode->SetProperty("layer", mitk::IntProperty::New(1));
   newNode->SetProperty("opacity", mitk::FloatProperty::New(0.5));

   this->GetDefaultDataStorage()->Add(newNode);

   mitk::RenderingManager::GetInstance()->RequestUpdateAll();
   }
catch( itk::ExceptionObject & excep )
   {
   std::cerr << "Exception caught !" << std::endl;
   std::cerr << excep << std::endl;
   }


}


best regards