ITK系列33_ 阈值水平集算法对脑部MHA文件进行三维分割

桃扇骨 2023-01-22 06:55 51阅读 0赞

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

#include "itkImage.h"
    #include "itkThresholdSegmentationLevelSetImageFilter.h"
    
    #include "itkFastMarchingImageFilter.h"
    #include "itkBinaryThresholdImageFilter.h"
    #include "itkImageFileReader.h"
    #include "itkImageFileWriter.h"
    #include "itkZeroCrossingImageFilter.h"
    
    
    int main( int argc, char *argv[] )
    {
      /*if( argc < 8 )
        {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: " << argv[0];
        std::cerr << " inputImage  outputImage";
        std::cerr << " seedX seedY InitialDistance";
        std::cerr << " LowerThreshold";
        std::cerr << " UpperThreshold";
        std::cerr << " [CurvatureScaling == 1.0]";
        std::cerr << std::endl;
        return EXIT_FAILURE;
        }*/
        /*现在我们使用一个特殊的像素类型和维来定义图像类型。在这种情况下我们将使用二维
            浮点型图像*/
      typedef   float           InternalPixelType;
      const     unsigned int    Dimension = 3;
      typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
      //接下来的几行使用New()方式实例化一个ThresholdSegmentationLevelSetImageFilter：
      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.0 );
    
      thresholder->SetOutsideValue(  0  );
      thresholder->SetInsideValue(  255 );
    
      typedef  itk::ImageFileReader< InternalImageType > ReaderType;
      typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
    
      ReaderType::Pointer reader = ReaderType::New();
      WriterType::Pointer writer = WriterType::New();
      //输入图像
      reader->SetFileName( "BrainProtonDensity3Slices.mha" );
      //保存图像
      writer->SetFileName( "naoshi_Threshold_LevelSet.mha" );
    
      typedef  itk::FastMarchingImageFilter< InternalImageType, InternalImageType >
      FastMarchingFilterType;
    
      FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
    
      typedef  itk::ThresholdSegmentationLevelSetImageFilter< InternalImageType,
        InternalImageType > ThresholdSegmentationLevelSetImageFilterType;
      ThresholdSegmentationLevelSetImageFilterType::Pointer thresholdSegmentation =
        ThresholdSegmentationLevelSetImageFilterType::New();
      /*对于ThresholdSegmentationLevelSetImageFilter，缩放比例参数用来平衡从等式（9 - 3）
          中的传播(膨胀)和曲率(表面平滑)系数的影响。这个滤波器中不使用水平对流系数。使用
          SetPropagationScaling()和SetCurvatureScaling()方式来设置这些系数。在这个例子中这两个
          系数都设置为1.0*/
      thresholdSegmentation->SetPropagationScaling( 1.0 );
      /*if ( argc > 8 )
        {
        thresholdSegmentation->SetCurvatureScaling( atof(argv[8]) );
        }
      else
        {*/
        thresholdSegmentation->SetCurvatureScaling( 1.0 );
      /*  }*/
    
        thresholdSegmentation->SetMaximumRMSError( 0.02 );
        thresholdSegmentation->SetNumberOfIterations( 1200 );//设置最大迭代次数
      /*  收敛标准MaximumRMSError和MaximumIterations设置的和前面例子中的一样。现在我
            们设置上下门限U和L，以及在初始化模型时使用的等值面值*/
        //上门限值（U）
      thresholdSegmentation->SetUpperThreshold( ::atof("250") );
      //下门限值（L）
      thresholdSegmentation->SetLowerThreshold( ::atof("210") );
      thresholdSegmentation->SetIsoSurfaceValue(0.0);
    
      thresholdSegmentation->SetInput( fastMarching->GetOutput() );
      thresholdSegmentation->SetFeatureImage( reader->GetOutput() );
      thresholder->SetInput( thresholdSegmentation->GetOutput() );
      writer->SetInput( thresholder->GetOutput() );
    
      typedef FastMarchingFilterType::NodeContainer           NodeContainer;
      typedef FastMarchingFilterType::NodeType                NodeType;
    
      NodeContainer::Pointer seeds = NodeContainer::New();
    
      InternalImageType::IndexType  seedPosition;
      //设置种子点位置
      seedPosition[0] = atoi("84");
      seedPosition[1] = atoi("126");
      seedPosition[2] = atoi("2");
      //表面到种子点的最初距离
      const double initialDistance = atof("5");//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 );
      /*调用writer上的Updata()方法引发了管道的运行。通常在出现错误和抛出异议时，从一
          个try / catch模块调用updata：*/
      try
        {
        reader->Update();
        const InternalImageType * inputImage = reader->GetOutput();
        fastMarching->SetOutputRegion( inputImage->GetBufferedRegion() );
        fastMarching->SetOutputSpacing( inputImage->GetSpacing() );
        fastMarching->SetOutputOrigin( inputImage->GetOrigin() );
        fastMarching->SetOutputDirection( inputImage->GetDirection() );
        writer->Update();
        }
      catch( itk::ExceptionObject & excep )
        {
        std::cerr << "Exception caught !" << std::endl;
        std::cerr << excep << std::endl;
        return EXIT_FAILURE;
        }
    
      std::cout << std::endl;
      std::cout << "Max. no. iterations: " << thresholdSegmentation->GetNumberOfIterations() << std::endl;
      std::cout << "Max. RMS error: " << thresholdSegmentation->GetMaximumRMSError() << std::endl;
      std::cout << std::endl;
      std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
      std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;
    
      typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
      //fastMarching快速步进输出水平集
      InternalWriterType::Pointer mapWriter = InternalWriterType::New();
      mapWriter->SetInput( fastMarching->GetOutput() );
      mapWriter->SetFileName("fastMarchingImage.mha");
      mapWriter->Update();
      //阈值水平集分割的速度（传播系数P）图像
      InternalWriterType::Pointer speedWriter = InternalWriterType::New();
      speedWriter->SetInput( thresholdSegmentation->GetSpeedImage() );
      speedWriter->SetFileName("speedTermImage.mha");
      speedWriter->Update();
    
      return EXIT_SUCCESS;
    }

输入三维图像（BrainProtonDensity3Slices.mha）：

## ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 1][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 2][] ##

**切片1 切片2 切片3**

输出三维脑室图像（naoshi\_Threshold\_LevelSet.mha）：

## ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 3][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 4][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 5][] ##

**切片1 切片2 切片3**

fastMarching快速步进输出水平集（fastMarchingImage.mha）：

## ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 6][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 7][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 8][] ##

**切片1 切片2 切片3**

阈值水平集分割的速度（传播系数P）图像（speedTermImage.mha）：

## ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 9][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 10][] ![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 11][] ##

**切片1 切片2 切片3**

**ITK系列目录：**

**[1 ITK图像数据表达之图像][1 ITK]**

**[2 ITK图像处理之图像滤波][2 ITK]**

**[3 ITK图像处理之图像分割][3 ITK]**

**注：例程配套素材见系列目录**

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20201203105054261.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/2020090817513321.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 2]: /images/20221020/301e2c47c18c4386a328d1b6306703e1.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 3]: https://img-blog.csdnimg.cn/20201203110539347.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 4]: https://img-blog.csdnimg.cn/20201203110554263.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 5]: /images/20221020/4fdc605ad1f245379e0cf5cf1d59fbb0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 6]: https://img-blog.csdnimg.cn/20201203110705966.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 7]: https://img-blog.csdnimg.cn/20201203110723350.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 8]: /images/20221020/ffea3577c7eb437d9efc37ea98dbb8a6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 9]: https://img-blog.csdnimg.cn/20201203110925614.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 10]: https://img-blog.csdnimg.cn/20201203110945200.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 11]: /images/20221020/7cc62eee6c094cc99e706e10b7019b19.png
[1 ITK]: https://blog.csdn.net/qq_32809093/article/details/108307568
[2 ITK]: https://blog.csdn.net/qq_32809093/article/details/108882952
[3 ITK]: https://blog.csdn.net/qq_32809093/article/details/108327734