ITK系列32_ 阈值水平集算法对脑部PNG图像进行二维分割

我不是女神ヾ 2023-01-22 06:55 58阅读 0赞

## **阈值水平集分割** ##

**itk::ThresholdSegmentationLevelSetImageFilter** 是对**阈值连接成员分割在水平集框架上的一个拓展**。目标是**定义一个亮度值的范围**来**对相关的组织类型继续分类**，**然后求出对那个亮度范围基于水平集等式上的传播系数**。使用水平集方法，进化的表面的平滑可以被用来阻止在连接成员方案中常见的“漏泄”。依照下面的公式可以从 FeatureImage 输入 g（x） 同 UpperThreshold 值 U 和 LowerThreshold 值 L 计算出下等式的传播系数 P(x) ：

![20201018095927171.png][] （9-4）

**上式阐述了传播系数函数P(x),在g中的亮度值在L和H之间是正数值P，而在亮度范围之外是负数值P**。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70][]

**从公式（9-4）得到的基于阈值水平集分割的传播系数**

ThresholdSegmentationFilter需要两个输入，第一个是itk::Image形式的**一个最初的水平集**，第二个是**一个特征图像**。对于许多应用，这个滤波器需要很少甚至不需要处理它的输入。平滑输入图像并不一定非要产生合理的解答，在一些情况下也是许可的。如下图所示展示了图像处理途径结构。**使用快速行进滤波器来生成最初的表面**。**分割滤波器的输出传递给一个itk::BinaryThresholdImageFilter来创建一个被分割对象的二值表示**。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 1][]

现在我们使用和9.1.1小节中**itk::ConnectedThresholdImageFilter**例子中使用的相同的数据和参数来运行这个例子。我们将使用一个数值为**5**的值来作为**表面到种子点的最初距离**。这个算法对这个初始化感觉相对迟缓。把图9-26所示中的结果和图9-1所示的结果相比较。**注意：表面的平滑约束是如何阻止分割泄露进入两个脑室的，还有灰质部分的分割局部化到了一个更小的部分中。**如表9-9所示为分割使用的参数。  
![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 2][]

### **实例32 阈值水平集****算法****对脑部PNG****图像进行二维分割** ###

#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 = 2;
      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( "BrainProtonDensitySlice.png" );
      //保存图像
      writer->SetFileName( "naoshi_Threshold_LevelSet.png" );
    
      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("81");
      seedPosition[1] = atoi("112");
      //表面到种子点的最初距离
      const double initialDistance = atof("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;
    }

**左脑室分割结果（种子点（81,112））：**

![20201018103539376.png][] ![20201018103516390.png][] ![20201018103737197.png][] ![20201018104428834.png][]![20201018104800442.png][]

**输入图像 左脑室分割图像 快速步进生产输入水平集 传播系数P(x)图像 打印输出**

**白质分割结果1（取左边种子点（60,116））：**

![20201018110237873.png][] ![20201018110339369.png][] ![20201018110632299.png][] ![20201018110143649.png][]

**左白质水平集分割图像 快速步进生产输入水平集 传播系数P(x)图像 打印输出**

**白质分割结果2（取右边种子点（60,116））：**

![20201018111644405.png][] ![20201018111720214.png][] ![20201018111810350.png][]

**右白质水平集分割图像 快速步进生产输入水平集 传播系数P(x)图像**

**右脑室分割结果（种子点（129,116））：**  **灰质分割结果（种子点（107,69）、（94,190））：**

![20201018105336302.png][] ![20201018113020554.png][] ![20201018113425630.png][]

结论：此方法对脑室分割效果较好，但对于灰质和白质，效果不是很好

**ITK系列目录：**

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

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

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

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

[20201018095927171.png]: /images/20221020/72c49feebd19450b9fa075fe159131a7.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70]: /images/20221020/ae0f3f939fa1444a82e7a2e99860a201.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 1]: /images/20221020/590a6283a0914ee799052605760f3ac1.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyODA5MDkz_size_16_color_FFFFFF_t_70 2]: /images/20221020/0164189d7f794cf4a489b1b616fa9299.png
[20201018103539376.png]: https://img-blog.csdnimg.cn/20201018103539376.png
[20201018103516390.png]: https://img-blog.csdnimg.cn/20201018103516390.png
[20201018103737197.png]: https://img-blog.csdnimg.cn/20201018103737197.png
[20201018104428834.png]: https://img-blog.csdnimg.cn/20201018104428834.png
[20201018104800442.png]: /images/20221020/f531d05a46a048d081da5b29339dfc63.png
[20201018110237873.png]: https://img-blog.csdnimg.cn/20201018110237873.png
[20201018110339369.png]: https://img-blog.csdnimg.cn/20201018110339369.png
[20201018110632299.png]: https://img-blog.csdnimg.cn/20201018110632299.png
[20201018110143649.png]: /images/20221020/5d7aa6221b7d4ebd88544d62bfa9e5b0.png
[20201018111644405.png]: https://img-blog.csdnimg.cn/20201018111644405.png
[20201018111720214.png]: https://img-blog.csdnimg.cn/20201018111720214.png
[20201018111810350.png]: /images/20221020/97acc75f7db745d7b511b6e3c0f31763.png
[20201018105336302.png]: https://img-blog.csdnimg.cn/20201018105336302.png
[20201018113020554.png]: https://img-blog.csdnimg.cn/20201018113020554.png
[20201018113425630.png]: /images/20221020/6be01e873700405c9a6c45a740a23649.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