[Insight-users] I am learning GeodesicActiveContourImageFilter , and I have some problem , please help me

[zhq]

Hi ,everyone 
     Using the code which included in folder Example/Segmentation , GeodesicActiveContourImageFilter , I can segment a jpg format picture . But if I want to segment a dicom file , a error happen . Please help me !  My code is shown as follow . 
    Thank you in advance !
zhq


#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"


int main()
{
  typedef   float           InternalPixelType;
  const     unsigned int    Dimension = 2;
  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 );


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


  typedef  itk::ImageFileReader< InternalImageType > ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
//  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->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.png");




// *************  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::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
                InternalImageType >    GeodesicActiveContourFilterType;
  GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
                                     GeodesicActiveContourFilterType::New();
//  geodesicActiveContour->Update();
//  std::cout<<"geodesicActiveContour new updata"<<std::endl;


  const double propagationScaling = 0.4;


  geodesicActiveContour->SetPropagationScaling( propagationScaling );
//  geodesicActiveContour->Update();
//  std::cout<<"geodesicActiveContour SetPropagationScaling updata"<<std::endl;
  geodesicActiveContour->SetCurvatureScaling( 1.0 );
  geodesicActiveContour->SetAdvectionScaling( 1.0 );
//  geodesicActiveContour->Update();
//  std::cout<<"geodesicActiveContour SetAdvectionScaling updata"<<std::endl;


  geodesicActiveContour->SetMaximumRMSError( 0.02 );
  geodesicActiveContour->SetNumberOfIterations( 800 );
//  geodesicActiveContour->Update();
//  std::cout<<"geodesicActiveContour SetNumberOfIterations updata"<<std::endl;


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


  geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
  geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
//  geodesicActiveContour->Update();
//  std::cout<<"geodesicActiveContour SetFeatureImage updata"<<std::endl;


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


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


  const double alpha =  -0.2;
  const double beta  =  10;


  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 ; // 81;
  seedPosition[1] = 443 ; // 114;


  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(
           reader->GetOutput()->GetBufferedRegion().GetSize() );




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


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


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


  geodesicActiveContour->Update();
  std::cout<<"geodesicActiveContour update"<<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/20130810/27a36cf3/attachment-0001.htm>