[Insight-users] CURVES Level Set Segmentation, Fast Marching Gives Empty Volume

[baris kandemir]

Hi ITK users,

I am trying to implement CURVES method(Lorigo et. al., 2001). I working
with 3D data set and trying to modify the example for 2D. I believe I got
the feature image properly. However I can't get the fast marching method
work. So i can't get an initial curve. What do you suggest? Thank you in
advance.
----------------------------------------------------------------------------
Here is my code:

#include "itkImage.h"
#include "itkCurvesLevelSetImageFilter.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"
#include "itkCurvesLevelSetFunction.h"
#include "itkCurvesLevelSetImageFilter.h"
#include <math.h>
#include "itkExpNegativeImageFilter.h"

int main (int argc, char *argv[])
{

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

  typedef itk::ImageFileReader<InternalImageType> ReaderType;
  typedef itk::ImageFileWriter<ExternalImageType> WriterType;
  typedef itk::RescaleIntensityImageFilter<InternalImageType,
ExternalImageType> CasterType;
  typedef itk::RescaleIntensityImageFilter<InternalImageType,
InternalImageType> RescaleFilterType;
  typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<
                               InternalImageType,
                               InternalImageType >  GradientFilterType;
  typedef   itk::CurvatureAnisotropicDiffusionImageFilter<
                               InternalImageType,
                               InternalImageType >  SmoothingFilterType;
  typedef   itk::SigmoidImageFilter<
                               InternalImageType,
                               InternalImageType >  SigmoidFilterType;
  typedef   itk::ExpNegativeImageFilter<InternalImageType,
InternalImageType> ExponentialType;
  typedef
itk::FastMarchingImageFilter<InternalImageType,InternalImageType>
FastMarchingType;

  RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();
  CasterType::Pointer caster = CasterType::New();
  GradientFilterType::Pointer GradientFilter = GradientFilterType::New();
  ExponentialType::Pointer ExponentialFilter = ExponentialType::New();
  FastMarchingType::Pointer fastmarchingfilter = FastMarchingType::New();

  reader->SetFileName(argv[1]);
  GradientFilter->SetSigma(1.2);
  ExponentialFilter->SetFactor(1.0);
  GradientFilter->SetInput(reader->GetOutput());
  rescaleFilter->SetInput(GradientFilter->GetOutput());
  rescaleFilter->SetOutputMinimum(0);
  rescaleFilter->SetOutputMaximum(255);
//   InternalImageType::Pointer gradientimage = InternalImageType::New();
  InternalImageType::Pointer FeatureImage = InternalImageType::New();
  caster->SetInput(rescaleFilter->GetOutput());
  writer->SetFileName("gradientimage.mhd");
  writer->SetInput(caster->GetOutput());
  writer->Update();
  ExponentialFilter->SetInput(rescaleFilter->GetOutput());
  FeatureImage = ExponentialFilter->GetOutput();
  caster->SetInput(FeatureImage);
  writer->SetFileName("featureimage.mhd");
  writer->SetInput(caster->GetOutput());
  writer->Update();

//   Obtaining the initial curve
  typedef FastMarchingType::NodeContainer  NodeContainer;
  typedef FastMarchingType::NodeType       NodeType;

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

  InternalImageType::IndexType  seedPosition;

  seedPosition[0] =  110;
  seedPosition[1] =  110;
  seedPosition[2] =  90;
//   seedPosition[3] =  80;
//   seedPosition[4] =  250;
//   seedPosition[5] =  240;
  const double initialDistance = 10;
  NodeType node;
  const double seedValue =  initialDistance;
  node.SetValue( seedValue );
  node.SetIndex( seedPosition );
  seeds->Initialize();
  seeds->InsertElement( 0, node );
  seeds->InsertElement( 1, node );
  seeds->InsertElement(2,node);
  fastmarching->SetStoppingValue(  150 );
  fastmarchingfilter->SetTrialPoints(  seeds  );
  fastmarchingfilter->SetSpeedConstant( 1.0 );
  caster->SetInput(fastmarchingfilter->GetOutput());
  writer->SetInput( caster->GetOutput() );
  writer->SetFileName("fastmarchingoutput.mhd");
  writer->Update();



  return 0;
}

-- 
Baris KANDEMIR

Bogazici University
Student of  Electrical and Electronics Engineering
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120216/1e402e1f/attachment.htm>