[Insight-users] I am learning active contour algorithm , there is problem about itkGeodesicActiveContourLevelSetImageFilter

[zhq]

Hello , Dear:
     I am learning itkGeodesicActiveContourLevelSetImageFilter , and I can segment a 2 dimension picture . And , now , I am trying to segment a 3 dimension data . But , when I segment a 2 dimension picture , it costs me about one minute . But , it costs me more than 20 minutes to segment a 3 dimension data , the data consists in 4 two-dimensional pictures  . Is my code wrong ? 
#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"
#include "itkGDCMImageIO.h"
#include "itkJPEGImageIO.h"
#include "itkChangeInformationImageFilter.h"
#include "itkImage.h"
#include "itkImageSeriesReader.h"
#include "itkNumericSeriesFileNames.h"
#include "itkImageFileWriter.h"
int main()
{
  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( -80.0 );
  thresholder->SetUpperThreshold(     2 );


  thresholder->SetOutsideValue(  0  );
  thresholder->SetInsideValue(  255 );
//  *************  reader   ******************************  //


//  typedef  itk::ImageFileReader< InternalImageType > ReaderType;
//  ReaderType::Pointer reader = ReaderType::New();
//  reader->SetFileName("C:\\Users\\zhq\\Desktop\\VolumeData.vtk");
//  reader->SetFileName("C:\\Users\\zhq\\Desktop\\BrainProtonDensitySlice.png");
//  reader->SetImageIO(itk::GDCMImageIO::New());
//  reader->SetFileName("C:\\Users\\zhq\\Desktop\\data\\SNAP_CR\\E403434298\\E403434298S1901I301.dcm");
//  reader->SetImageIO(itk::JPEGImageIO::New());
//  reader->SetFileName("C:\\Users\\zhq\\Desktop\\data\\picture8\\1.jpeg");
//  reader->SetFileName("C:\\Users\\zhq\\Desktop\\\\1.jpg");


  typedef itk::ImageSeriesReader<InternalImageType> ReaderType ;
  ReaderType::Pointer reader = ReaderType::New();
  typedef itk::NumericSeriesFileNames NameGeneratorType ; 
  NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();


  nameGenerator->SetSeriesFormat("C:\\Users\\zhq\\Desktop\\data\\SNAP_CR\\E403434298\\E403434298S1901I%3d.dcm");
  nameGenerator->SetStartIndex(301);
  nameGenerator->SetEndIndex(304);
  nameGenerator->SetIncrementIndex(1);


  reader->SetImageIO(itk::GDCMImageIO::New());
  reader->SetFileNames(nameGenerator->GetFileNames());
  reader->Update();
  std::cout<<"reader update"<<std::endl;
// *************  finish reader  **********************    // 




// *************  writer  *******************************  // 
  typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName("C:\\Users\\zhq\\Desktop\\ActiveContour.vtk");




// *************  finish  writer   ********************** //


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


  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();


  typedef itk::ChangeInformationImageFilter<InternalImageType> ChangeInformationFilterType ; 
  ChangeInformationFilterType::Pointer changeInformationFilter = ChangeInformationFilterType::New();
  changeInformationFilter->ChangeSpacingOn();
  changeInformationFilter->ChangeOriginOn();
  changeInformationFilter->ChangeDirectionOn();


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


  const double propagationScaling = 0.25;


  geodesicActiveContour->SetPropagationScaling( propagationScaling );
  geodesicActiveContour->SetCurvatureScaling( 1 );
  geodesicActiveContour->SetAdvectionScaling( 0.5 );


  geodesicActiveContour->SetMaximumRMSError( 0.02 );
  geodesicActiveContour->SetNumberOfIterations( 800 );




  changeInformationFilter->SetInput(reader->GetOutput());
  gradientMagnitude->SetInput( changeInformationFilter->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );


  geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
  geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );


  thresholder->SetInput( geodesicActiveContour->GetOutput() );
  writer->SetInput( thresholder->GetOutput() );


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


  const double alpha =  -0.1;
  const double beta  =  100;


  sigmoid->SetAlpha( alpha );
  sigmoid->SetBeta(  beta  );


  typedef FastMarchingFilterType::NodeContainer  NodeContainer;
  typedef FastMarchingFilterType::NodeType       NodeType;


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


  InternalImageType::IndexType  seedPosition;


  seedPosition[0] = 370;//661 ; //370 ; // 81;
  seedPosition[1] = 443;//559 ; //443 ; // 114;
  seedPosition[2] = 1;


  const double initialDistance = 5;


  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 );


  fastMarching->SetOutputSize(
           changeInformationFilter->GetOutput()->GetBufferedRegion().GetSize() );
  try
  {
 gradientMagnitude->Update();
 std::cout<<"gradientMagnitude update"<<std::endl;
  }
  catch( itk::ExceptionObject & excep )
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
  }


  sigmoid->Update();
  std::cout<<"sigmoid update"<<std::endl;


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


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


  writer->Update();
  std::cout<<"writer update"<<std::endl;






  std::cout << std::endl;
  std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
  std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
  std::cout << std::endl;
  std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
  std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;


  return 0;
}


    
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130811/dc3b3e95/attachment.htm>