ITK 实例17 阈值水平集算法对脑部MHA文件进行三维分割

  1 #include "itkImage.h"
  2 #include "itkThresholdSegmentationLevelSetImageFilter.h"
  3  
  4 #include "itkFastMarchingImageFilter.h"
  5 #include "itkBinaryThresholdImageFilter.h"
  6 #include "itkImageFileReader.h"
  7 #include "itkImageFileWriter.h"
  8 #include "itkZeroCrossingImageFilter.h"
  9  
 10  
 11 int main( int argc, char *argv[] )
 12 {
 13   /*if( argc < 8 )
 14     {
 15     std::cerr << "Missing Parameters " << std::endl;
 16     std::cerr << "Usage: " << argv[0];
 17     std::cerr << " inputImage  outputImage";
 18     std::cerr << " seedX seedY InitialDistance";
 19     std::cerr << " LowerThreshold";
 20     std::cerr << " UpperThreshold";
 21     std::cerr << " [CurvatureScaling == 1.0]";
 22     std::cerr << std::endl;
 23     return EXIT_FAILURE;
 24     }*/
 25     /*现在我们使用一个特殊的像素类型和维来定义图像类型。在这种情况下我们将使用二维
 26         浮点型图像*/
 27   typedef   float           InternalPixelType;
 28   const     unsigned int    Dimension = 3;
 29   typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
 30   //接下来的几行使用New()方式实例化一个ThresholdSegmentationLevelSetImageFilter：
 31   typedef unsigned char                            OutputPixelType;
 32   typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
 33   typedef itk::BinaryThresholdImageFilter<InternalImageType, OutputImageType>
 34                                                    ThresholdingFilterType;
 35  
 36   ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
 37  
 38   thresholder->SetLowerThreshold( -1000.0 );
 39   thresholder->SetUpperThreshold(     0.0 );
 40  
 41   thresholder->SetOutsideValue(  0  );
 42   thresholder->SetInsideValue(  255 );
 43  
 44   typedef  itk::ImageFileReader< InternalImageType > ReaderType;
 45   typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
 46  
 47   ReaderType::Pointer reader = ReaderType::New();
 48   WriterType::Pointer writer = WriterType::New();
 49   //输入图像
 50   reader->SetFileName( "BrainProtonDensity3Slices.mha" );
 51   //保存图像
 52   writer->SetFileName( "naoshi_Threshold_LevelSet.mha" );
 53  
 54   typedef  itk::FastMarchingImageFilter< InternalImageType, InternalImageType >
 55   FastMarchingFilterType;
 56  
 57   FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
 58  
 59   typedef  itk::ThresholdSegmentationLevelSetImageFilter< InternalImageType,
 60     InternalImageType > ThresholdSegmentationLevelSetImageFilterType;
 61   ThresholdSegmentationLevelSetImageFilterType::Pointer thresholdSegmentation =
 62     ThresholdSegmentationLevelSetImageFilterType::New();
 63   /*对于ThresholdSegmentationLevelSetImageFilter，缩放比例参数用来平衡从等式（9 - 3）
 64       中的传播(膨胀)和曲率(表面平滑)系数的影响。这个滤波器中不使用水平对流系数。使用
 65       SetPropagationScaling()和SetCurvatureScaling()方式来设置这些系数。在这个例子中这两个
 66       系数都设置为1.0*/
 67   thresholdSegmentation->SetPropagationScaling( 1.0 );
 68   /*if ( argc > 8 )
 69     {
 70     thresholdSegmentation->SetCurvatureScaling( atof(argv[8]) );
 71     }
 72   else
 73     {*/
 74     thresholdSegmentation->SetCurvatureScaling( 1.0 );
 75   /*  }*/
 76  
 77     thresholdSegmentation->SetMaximumRMSError( 0.02 );
 78     thresholdSegmentation->SetNumberOfIterations( 1200 );//设置最大迭代次数
 79   /*  收敛标准MaximumRMSError和MaximumIterations设置的和前面例子中的一样。现在我
 80         们设置上下门限U和L，以及在初始化模型时使用的等值面值*/
 81     //上门限值（U）
 82   thresholdSegmentation->SetUpperThreshold( ::atof("250") );
 83   //下门限值（L）
 84   thresholdSegmentation->SetLowerThreshold( ::atof("210") );
 85   thresholdSegmentation->SetIsoSurfaceValue(0.0);
 86  
 87   thresholdSegmentation->SetInput( fastMarching->GetOutput() );
 88   thresholdSegmentation->SetFeatureImage( reader->GetOutput() );
 89   thresholder->SetInput( thresholdSegmentation->GetOutput() );
 90   writer->SetInput( thresholder->GetOutput() );
 91  
 92   typedef FastMarchingFilterType::NodeContainer           NodeContainer;
 93   typedef FastMarchingFilterType::NodeType                NodeType;
 94  
 95   NodeContainer::Pointer seeds = NodeContainer::New();
 96  
 97   InternalImageType::IndexType  seedPosition;
 98   //设置种子点位置
 99   seedPosition[0] = atoi("84");
100   seedPosition[1] = atoi("126");
101   seedPosition[2] = atoi("2");
102   //表面到种子点的最初距离
103   const double initialDistance = atof("5");//5
104  
105   NodeType node;
106  
107   const double seedValue = - initialDistance;
108  
109   node.SetValue( seedValue );
110   node.SetIndex( seedPosition );
111  
112   seeds->Initialize();
113   seeds->InsertElement( 0, node );
114  
115   fastMarching->SetTrialPoints(  seeds  );
116  
117   fastMarching->SetSpeedConstant( 1.0 );
118   /*调用writer上的Updata()方法引发了管道的运行。通常在出现错误和抛出异议时，从一
119       个try / catch模块调用updata：*/
120   try
121     {
122     reader->Update();
123     const InternalImageType * inputImage = reader->GetOutput();
124     fastMarching->SetOutputRegion( inputImage->GetBufferedRegion() );
125     fastMarching->SetOutputSpacing( inputImage->GetSpacing() );
126     fastMarching->SetOutputOrigin( inputImage->GetOrigin() );
127     fastMarching->SetOutputDirection( inputImage->GetDirection() );
128     writer->Update();
129     }
130   catch( itk::ExceptionObject & excep )
131     {
132     std::cerr << "Exception caught !" << std::endl;
133     std::cerr << excep << std::endl;
134     return EXIT_FAILURE;
135     }
136  
137   std::cout << std::endl;
138   std::cout << "Max. no. iterations: " << thresholdSegmentation->GetNumberOfIterations() << std::endl;
139   std::cout << "Max. RMS error: " << thresholdSegmentation->GetMaximumRMSError() << std::endl;
140   std::cout << std::endl;
141   std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
142   std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;
143  
144   typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
145   //fastMarching快速步进输出水平集
146   InternalWriterType::Pointer mapWriter = InternalWriterType::New();
147   mapWriter->SetInput( fastMarching->GetOutput() );
148   mapWriter->SetFileName("fastMarchingImage.mha");
149   mapWriter->Update();
150   //阈值水平集分割的速度（传播系数P）图像
151   InternalWriterType::Pointer speedWriter = InternalWriterType::New();
152   speedWriter->SetInput( thresholdSegmentation->GetSpeedImage() );
153   speedWriter->SetFileName("speedTermImage.mha");
154   speedWriter->Update();
155  
156   return EXIT_SUCCESS;
157 }