ITK系列30_ 水平集(快速步进)算法对脑部PNG图像进行二维分割

水平集分割

水平集是跟踪轮廓和表面运动的一种数字化方法。不直接对轮廓进行操作，而是将轮廓设置成一个高维函数的零水平集，这个高维函数叫做水平集函数： Ψ(X,t) 。然后水平集函数运动成为一个微分方程。在任何时候，通过从输出中提取零水平集来得到运动的轮廓。使用水平集的主要优点是可以对任何复杂的结构进行模式化和拓扑变换，比如暗中操作融合和分离。通过使用基于图像的诸如亮度均值、梯度和边缘之类的特征的微分方程的解答，水平集就可以用来对图像进行分割。在一个典型的方法中，用户对一个轮廓进行初始化并运动这个轮廓直到它符合图像中的一个解剖结构。在这个文献中，发表了许多这些基本概念的执行和变量。赛思 Sethian 已经对这一领域做了一个综述。

快速步进分割

当微分方程控制水平集运动时有一个非常简单的形式，可以使用一个快速运动算法叫做快速步进。接下来的例子阐述了 itk::FastMarchingImageFilter 的用法。这个滤波器对一个简单水平集运动问题执行一个快速行进方法。在这个例子中，微分方程中使用的速率系数是通过用户以一种图像的形式来提供的。通常作为一个梯度值的函数来计算这个图像。在文献中有几个映射是很常见的，比如负指数 exp( − x) 和倒数 1/(1+x) 。在当前的例子中我们决定使用一种 S 函数Sigmoid 函数，因为它可以提供一个控制参数的好方法，这个参数可以定制成一个很好的速率图像形状。这个映射按以下方式进行：当图像梯度很高时前面的传播速率很慢而在梯度很低的区域传播速率很快。这种方式将导致轮廓不断传播直到它到达图像中的解剖结构边缘并在那些边
缘前将速率降低下来。 FastMarchingImageFilter 的输出是一个时间交叉图，对每一个像素，它还表达了到达这个像素位置之前所花费的时间。然后在输出图像中一个门限的应用等同于在轮廓运动过程中的一个特殊时间上对轮廓取一个快照。轮廓横渡到一个特殊解剖结构边缘被认为将花费更长的时间。这将导致在接近结构边缘的时间交叉像素值发生很大的变化。通过调用一个时间范围用这个滤波器来执行分割，在这个时间范围内图像空间中一个区域将在一段长时间内包含轮廓。如图 9-15 所示显示了包含在一个分割任务中的 FastMarchingImageFilter 的应用的主要成员。它包括一个使用 itk::CurvatureAnisotropicDiffusionImageFilter 进行平滑的最初阶段。平滑后的图像作为输入传递给 itk::GradientMagnitudeRecursiveGaussianImageFilter ，然后再给
itk::SigmoidImageFilter 。 最 后 ， 将 FastMarchingImageFilter 的 输 出 传 递 给 itk::BinaryThresholdImageFilter 以便产生一个 binary mask 来表达分割对象。

接下来例子中的代码阐述了一个用来执行快速步进分割的管道的典型设置。首先，使用一个边缘保护滤波器对输入图像进行平滑；然后计算它的梯度值并传递给一个 sigmoid 滤波器。 sigmoid 滤波器的结果是图像的潜能，将用来影响微分方程的速率系数。

实例30 快速步进算法对脑部PNG图像进行二维分割

//包含用来从输入图像中去除噪声头文件
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
//这两个滤波器连在一起将产生调节描述水平集运动的微分方程中的速率系数的图像潜能。
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"

#include "itkFastMarchingImageFilter.h"
//从 FastMarchingImageFilter 得到的时间交叉图将使用 BinaryThresholdImageFilter 来进行门限处理
#include "itkBinaryThresholdImageFilter.h"
//使用 itk::ImageFileReader 和 itk::ImageFileWriter 来对图像进行读写
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"

static void PrintCommandLineUsage( const int argc, const char * const argv[] )
{
  std::cerr << "Missing Parameters " << std::endl;
  std::cerr << "Usage: " << argv[0];
  std::cerr << " inputImage  outputImage seedX seedY";
  std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue";
  std::cerr << " smoothingOutputImage gradientMagnitudeOutputImage sigmoidOutputImage" << std::endl;

  for (int qq=0; qq< argc; ++qq)
    {
    std::cout << "argv[" << qq << "] = " << argv[qq] << std::endl;
    }
}

int main( int argc, char *argv[] )
{
  /*if (argc != 13)
    {
    PrintCommandLineUsage(argc, argv);
    return EXIT_FAILURE;
    }*/

  /*现在我们使用一个像素类型和一个特殊维来定义图像类型。由于需要用到平滑滤波器，
   所以在这种情况下对像素使用浮点型数据类型*/
  typedef   float           InternalPixelType;
  const     unsigned int    Dimension = 2;
  typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
  //另一方面，输出图像被声明为二值的
  typedef unsigned char                            OutputPixelType;
  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
 //下面使用图像内在的类型和输出图像类型来对 BinaryThresholdImageFilter 的类型进行实
  //例化：
  typedef itk::BinaryThresholdImageFilter< InternalImageType,
                        OutputImageType    >    ThresholdingFilterType;
  ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
  //设置阈值
  const InternalPixelType  timeThreshold = atof("100");
  /*传递给 BinaryThresholdImageFilter 的上门限将定义那个从时间交叉图得来的时间快照。
  在一个理想的应用中，用户拥有使用视觉反馈交互式的选择这个门限的能力。在这里，由于
  是一个最小限度的例子，从命令行来得到这个值*/
  thresholder->SetLowerThreshold(           0.0  );
  thresholder->SetUpperThreshold( timeThreshold  );

  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("BrainProtonDensitySlice_FastMarching.png");

  typedef itk::RescaleIntensityImageFilter<
                               InternalImageType,
                               OutputImageType >   CastFilterType;
  //使用图像内在的类型对 CurvatureAnisotropicDiffusionImageFilter 类型进行实例化
    typedef   itk::CurvatureAnisotropicDiffusionImageFilter<
                               InternalImageType,
                               InternalImageType >  SmoothingFilterType;
  //然后，通过调用 New( ) 方式来创建滤波器并将结果指向一个 itk::SmartPointer
  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
  /*使用图像内在的类型来对 GradientMagnitudeRecursiveGaussianImageFilter 和
   SigmoidImageFilter 的类型进行实例化*/
  typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<
                               InternalImageType,
                               InternalImageType >  GradientFilterType;
  typedef   itk::SigmoidImageFilter<
                               InternalImageType,
                               InternalImageType >  SigmoidFilterType;
  //使用 New( ) 方式来对相应的滤波器进行实例化
  GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
  /*使用 SetOutputMinimum() 和 SetOutputMaximum() 方式来定义 SigmoidImageFilter 输出的
  最小值和最大值。在我们的例子中，我们希望这两个值分别是 0.0 和 1.0 ，以便得到一个最佳
  的速率图像来给 MarchingImageFilter 。在 6.3.2 小节详细地表达了 SigmoidImageFilter 的使用方
  法。*/
  sigmoid->SetOutputMinimum(  0.0  );
  sigmoid->SetOutputMaximum(  1.0  );
  //现在我们声明 FastMarchingImageFilter 的类型：
  typedef  itk::FastMarchingImageFilter< InternalImageType,
                              InternalImageType >    FastMarchingFilterType;
  //然后我们使用 New( ) 方式结构化这个类的一个滤波器
  FastMarchingFilterType::Pointer  fastMarching
                                              = FastMarchingFilterType::New();
  //现在我们使用下面的代码在一个管道中连接这些滤波器
  smoothing->SetInput( reader->GetOutput() );
  gradientMagnitude->SetInput( smoothing->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );
  fastMarching->SetInput( sigmoid->GetOutput() );
  thresholder->SetInput( fastMarching->GetOutput() );
  writer->SetInput( thresholder->GetOutput() );
  /*CurvatureAnisotropicDiffusionImageFilter 类需要定义一对参数。下面是二维图像中的典
  型值。然而它们可能需要根据输入图像中存在的噪声的数量进行适当的调整。在 6.7.3 小节中
  讨论了这个滤波器*/
  smoothing->SetTimeStep( 0.125 );
  smoothing->SetNumberOfIterations(  5 );
  smoothing->SetConductanceParameter( 9.0 );
  //设置Sigma(σ )
  const double sigma = atof("1.0");
  /*执行 GradientMagnitudeRecursiveGaussianImageFilter 就等同于通过一个派生的操作符紧
  跟一个高斯核的一个回旋。这个高斯的 sigma(σ) 可以用来控制图像边缘影响的范围。在 6.4.2
  小节中讨论了这个滤波器*/
  gradientMagnitude->SetSigma(  sigma  );
  //设置SigmoidAlpha(α)   SigmoidBeta(β)
  const double alpha =  atof("-0.5");
  const double beta  =  atof("3.0");
  /*SigmoidImageFilter 类需要两个参数来定义用于 sigmoid 变量的线性变换。使用 SetAlpha()
      和 SetBeta() 方式来传递这些参数。在这个例子的背景中，参数用来强化速率图像中的高低值
      区域之间的区别。在理想情况下，解剖结构中的均匀区域内的速率应该是 1.0 而在结构边缘
      附近应该迅速地衰减到 0.0 。下面是寻找这个值的启发。在梯度量值图像中，我们将沿被分
      割的解剖结构地轮廓的最小值称为 K1 ，将结构中间的梯度强度的均值称为 K2 。这两个值表
      达了我们想要映射到速率图像中间距为[0:1] 的动态范围。我们期望 sigmoid 将 K1 映射为 0.0 ，
      K2 为 1.0 。由于 K1 的值要比 K2 更大而我们期望将它们分别映射为 0.0 和 1.0 ，所以我们对 alpha(α)
      选择一个负值以便 sigmoid 函数也做一个反转亮度映射。这个映射将产生一个速率图像，其
      中水平集将在均匀区域快速行进而在轮廓上停止。 Beta(β) 的参考值是(K1 + K2) / 2 而 α 的参考值
      是(K2 - K1) / 6 必须是一个负数。在我们的简单例子中，这些值是用户从命令行提供的。用户
      可以通过留心梯度量值图像来估计这些值*/
  sigmoid->SetAlpha( alpha );
  sigmoid->SetBeta(  beta  );
  /*FastMarchingImageFilter 需要用户提供一个轮廓扩张的种子点。用户实际上不仅要提供
      一个种子点而是一个种子点集。一个好的种子点集将增大分割一个复杂对象而不丢失数据的
      可能性。使用多种子也会减少访问整个对象所需要的时间同时减少前面访问的区域边缘的漏
      洞的风险。例如，当分割一个拉长的对象时，在其中一端设置一个种子将导致传播到另一端
      会需要很长一段时间。沿着轴线设置几个种子无疑将是确定整个对象尽早被捕获的最佳策
      略。水平集的一个重要的性质就是它们不需要任何额外辅助来融合几个起点的自然能力。使
      用多种子正是利用了这个性质。
      这些种子储存在一个容器中。在 FastMarchingImageFilter 的特性中将这个容器的类型定
      义为 NodeContainer ：*/
  typedef FastMarchingFilterType::NodeContainer           NodeContainer;
  typedef FastMarchingFilterType::NodeType                NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();
  
  InternalImageType::IndexType  seedPosition;
  //设置种子点位置
  seedPosition[0] = atoi("99");
  seedPosition[1] = atoi("114");
  //节点作为堆栈变量来创建，并使用一个值和一个 itk::Index 位置来进行初始化
  NodeType node;
  const double seedValue = 0.0;

  node.SetValue( seedValue );
  node.SetIndex( seedPosition );
  //节点列表被初始化，然后使用 InsertElement( ) 来插入每个节点：
  seeds->Initialize();
  seeds->InsertElement( 0, node );
  //现在使用 SetTrialPoints( ) 方式将种子节点集传递给 FastMarchingImageFilter ：
  fastMarching->SetTrialPoints(  seeds  );
  /*FastMarchingImageFilter 需要用户指定作为输出产生的图像的大小。使用 SetOutputSize()
      来完成。注意：这里的大小是从平滑滤波器的输出图像得到的。这个图像的大小仅仅在直接
      或间接地调用了这个滤波器的 Updata() 之后才是有效的。*/
  //writer 上的 Updata( ) 方法引发了管道的运行。通常在出现错误和抛出异议时 , 从一个
  //smoothingOutputImage
  //使用边缘保护平滑滤波器平滑滤波后的图像
  try
    {
    CastFilterType::Pointer caster1 = CastFilterType::New();
    WriterType::Pointer writer1 = WriterType::New();
    caster1->SetInput( smoothing->GetOutput() );
    writer1->SetInput( caster1->GetOutput() );
    //
    writer1->SetFileName("smoothingOutputImage.png");
    caster1->SetOutputMinimum(   0 );
    caster1->SetOutputMaximum( 255 );
    writer1->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }
  // gradientMagnitudeOutputImage
  //滤波后图像的梯度强度图像
  try
    {
    CastFilterType::Pointer caster2 = CastFilterType::New();
    WriterType::Pointer writer2 = WriterType::New();
    caster2->SetInput( gradientMagnitude->GetOutput() );
    writer2->SetInput( caster2->GetOutput() );
    writer2->SetFileName("gradientMagnitudeOutputImage.png");
    caster2->SetOutputMinimum(   0 );
    caster2->SetOutputMaximum( 255 );
    writer2->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }
  // sigmoidOutputImage
  //梯度强度的 sigmoid 函数图像
  try
    {
    CastFilterType::Pointer caster3 = CastFilterType::New();
    WriterType::Pointer writer3 = WriterType::New();
    caster3->SetInput( sigmoid->GetOutput() );
    writer3->SetInput( caster3->GetOutput() );
    // FastMarchingImageFilter 分割处理产生的结果图像
    writer3->SetFileName("sigmoidOutputImage.png");
    caster3->SetOutputMinimum(   0 );
    caster3->SetOutputMaximum( 255 );
    writer3->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  
  fastMarching->SetOutputSize(
           reader->GetOutput()->GetBufferedRegion().GetSize() );
  //设置stoppingTime，比阈值（门限值高一点）
  const double stoppingTime = atof("103");
  /*由于表示轮廓的起点将从头到尾不断传播，一旦到达一个特定的时间就希望停止这一过
      程。这就允许我们在假定相关区域已经计算过的假设下节省计算时间。使用 SetStopping
      Value() 方式来定义停止过程的值。原则上，这个停止值应该比门限值高一点*/
  fastMarching->SetStoppingValue(  stoppingTime  );
  //  try / catch 模块调用 updata 。
  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    return EXIT_FAILURE;
    }

  try
    {
    CastFilterType::Pointer caster4 = CastFilterType::New();
    WriterType::Pointer writer4 = WriterType::New();
    caster4->SetInput( fastMarching->GetOutput() );
    writer4->SetInput( caster4->GetOutput() );
    writer4->SetFileName("FastMarchingFilterOutput4.png");
    caster4->SetOutputMinimum(   0 );
    caster4->SetOutputMaximum( 255 );
    writer4->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;

  InternalWriterType::Pointer mapWriter = InternalWriterType::New();
  mapWriter->SetInput( fastMarching->GetOutput() );
  mapWriter->SetFileName("FastMarchingFilterOutput4.mha");
  mapWriter->Update();
 
  InternalWriterType::Pointer speedWriter = InternalWriterType::New();
  speedWriter->SetInput( sigmoid->GetOutput() );
  speedWriter->SetFileName("FastMarchingFilterOutput3.mha");
  speedWriter->Update();

  InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
  gradientWriter->SetInput( gradientMagnitude->GetOutput() );
  gradientWriter->SetFileName("FastMarchingFilterOutput2.mha");
  gradientWriter->Update();

  return EXIT_SUCCESS;
}

现在我们使用目录 Examples/Data 中提供的 BrainProtonDensitySlice.png 作为输入图像来运行这个例子。我们可以通过选取适当的位置作为种子来很便捷地对主要的解剖学结构进行分割。表 9-6 所示表达了一些结构使用的参数。

使用的终止值都为 100 。如图 9-16 所示表达了图 9-15 所示中阐述的流程的中间输出。从左到右分别是：各向异性扩散滤波器的输出，平滑后图像的梯度强度和最后作为 FastMarchingImageFilter 的速率图像使用的梯度强度的 sigmoid 。

 

                                         

        输入图像             使用边缘保护平滑滤波器平滑滤波后的图像   滤波后图像的梯度强度图像     梯度强度的 sigmoid 函数图像

基于 FastMarchingImageFilter 分割处理产生的结果：

                                    

       左脑室分割图像                          右脑室分割图像                           白质分割图像                          灰质分割图像

注意：图像中的灰质部分并没有被完全分割，这就阐述了在对被分割的解剖结构并没有占据图像的延伸区域进行分割时水平集方法的弱点。这在当结构的宽度和梯度滤波器所产生的衰减带的大小做比较时显得尤为真实。对这个限制可行的一个工作区是沿着拉长的对象使用多种子分布，然而，白质相对于灰质分割来说是一个价值不高的工作，可能需要一个比这个简单例子中使用的更精细的方法。

 

ITK系列目录：

1 ITK图像数据表达之图像

2 ITK图像处理之图像滤波

3 ITK图像处理之图像分割

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