使用ITK的itkGeodesicActiveContourLevelSetImageFilter时，报错，求解释，求改正

[yingzijuntuan]

是对dicom图像处理时才出现这种问题，对jpg或者png能够分割，代码如下：

#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->SetFileName("C:\\Users\\zhq\\Desktop\\\\1.jpg");
  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();

  const double propagationScaling = 0.4;

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

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

  gradientMagnitude->SetInput( reader->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.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;//661 ; //370 ; // 81;
  seedPosition[1] = 443;//559 ; //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->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;
}

[aimi19930326]

使用ScaleTransform改变输入图像各方向的spacing值，设置比例使其与轮廓模型默认的一致，轮廓模型默认的是1.000000.

[aimi19930326]

这是我的处理过程，你看看有借鉴价值没，我的是三维的。 InternalImageType::SizeType size; size[0] = 512; size[1] = 512; size[2] = 50; typedef itk::ScaleTransform<double, 3> TransformType; TransformType::Pointer scaleTransform = TransformType::New(); itk::FixedArray<float, 3> scale; scale[0] = 1.505883; // newWidth/oldWidth scale[1] = 1.505883; scale[2] = 1.000000; scaleTransform->SetScale(scale); itk::Point<float, 3> center; center[0] = size[0]/ 2; center[1] = size[1] / 2; center[2] = size[2] / 2; scaleTransform->SetCenter(center); typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType; ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New(); resample->SetInput(sToc->GetOutput()); resample->SetSize(size); resample->SetTransform(scaleTransform); resample->Update(); InternalImageType::Pointer inputimage = resample->GetOutput();