ITK系列22_ 区域增长(置信连接)算法对脑部PNG图像进行二维分割

置信连接

接下来的例子阐述了 itk::ConfidenceConnectedImageFilter 的用法。 ConfidenceConnected Image Filter 使用的标准是基于当前区域的简单统计上的。首先，算法计算包含在区域中的所有像素亮度的平均值和标准差。用户提供一个因子用来乘以标准差并定义一个平均值的范围。相邻像素中亮度值在这个范围内的将包含到这个区域中。当没有更多的像素符合这个标准时，算法将结束它的第一次迭代。使用包含在区域内的所有像素再次计算亮度值的平均值和标准差。这个平均值和标准差定义一个新的亮度范围，用来查看当前区域的邻域并评价它们的亮度是否落在这个范围内。重复这个迭代过程直到没有新的像素再加进来或者已经达到了迭代器的最大数目。下面的式子阐述了这个滤波器的使用范围：

                                                                             I(X) ∈ [m− f  σ , m+ f  σ]

其中 m 和 σ 分别是区域亮度的平均值和标准差， f 是用户提供的一个因子， I( ) 是图像，X 是区域中特殊相邻像素的位置。

                  

实例22 置信连接算法对脑部PNG图像进行二维分割

#include "itkConfidenceConnectedImageFilter.h"//包含置信连接类
//图像中存在的噪声会降低这个滤波器生长大面积区域的能力。当面对噪声图像时，通常
//是使用一个边缘保留平滑滤波器
#include "itkCastImageFilter.h"//滤波器
#include "itkCurvatureFlowImageFilter.h"//边缘保留平滑滤波器
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

//ConfidenceConnected Image Filter 使用的标准是基于当前区域的简单统计上的。首先，算法计算包含在区域中的所
//有像素亮度的平均值和标准差。用户提供一个因子用来乘以标准差并定义一个平均值的范
//围。相邻像素中亮度值在这个范围内的将包含到这个区域中。当没有更多的像素符合这个标
//准时，算法将结束它的第一次迭代。使用包含在区域内的所有像素再次计算亮度值的平均值
//和标准差。这个平均值和标准差定义一个新的亮度范围，用来查看当前区域的邻域并评价它
//们的亮度是否落在这个范围内。重复这个迭代过程直到没有新的像素再加进来或者已经达到
//了迭代器的最大数目。下面的式子阐述了这个滤波器的使用范围：
//I(X) ∈[m− fσ, m + fσ]
//其中 m 和 σ 分别是区域亮度的平均值和标准差， f 是用户提供的一个因子， I() 是图像，
//X 是区域中特殊相邻像素的位置

int main( int argc, char *argv[] )
{
  //if( argc < 5 )
  //  {
  //  std::cerr << "Missing Parameters " << std::endl;
  //  std::cerr << "Usage: " << argv[0];
  //  std::cerr << " inputImage  outputImage seedX seedY " << std::endl;
  //  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;

  typedef itk::CastImageFilter< InternalImageType, OutputImageType >
    CastingFilterType;
  CastingFilterType::Pointer caster = CastingFilterType::New();

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

  //使用图像类型作为模板参数来对平滑滤波器进行实例化
  typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >
    CurvatureFlowImageFilterType;
  //调用 New( ) 方式来创建滤波器并将接指向 itk::SmartPointer
  CurvatureFlowImageFilterType::Pointer smoothing =
                         CurvatureFlowImageFilterType::New();

  /*现在我们声明区域生长滤波器的类型，本例中使用 ConnectedThresholdImageFilter*/
  typedef itk::ConfidenceConnectedImageFilter<
            InternalImageType, InternalImageType> ConnectedFilterType;
  //然后我们使用 New( ) 方式构造这种类的一个滤波器
  ConnectedFilterType::Pointer confidenceConnected = ConnectedFilterType::New();
  /*现在到了创建一个简单的线性管道的时候。在管道的开头添加一个 reader 件头并在末尾
  添加一个 cast filter 和 writer 。由于只有仅仅一小部分图像文件格式支持浮点型数据类型，
  所以使用 cast filter 将浮点型数据类型转换成整型*/
  smoothing->SetInput( reader->GetOutput() );//滤波
  confidenceConnected->SetInput( smoothing->GetOutput() );//区域增长置信连接
  caster->SetInput( confidenceConnected->GetOutput() );//类型转换，浮点型数据类型转换成整型
  writer->SetInput( caster->GetOutput() );
  
  //滤波
  /*CurvatureFlowImageFilter 需要定义两个参数。下面是一个二维图像的常见值。当然它们
  也需要根据输入图像存在的噪声的数量进行适当的调整*/
  smoothing->SetNumberOfIterations( 5 );
  smoothing->SetTimeStep( 0.125 );
  

  //区域增长-置信连接
  /*ConfidenceConnectedImageFilter 需要定义两个参数。第一个，定义亮度范围大小的因子
  f 。乘法因子太小将限制当前区域中很多有相似亮度的像素；乘法因子太大将放松接受条件
  并导致区域过度生长。值太大就会使得这些区域生长到邻近区域，而这些邻近区域实际上可
  能属于独立的解剖结构*/
  confidenceConnected->SetMultiplier( 2.5 );//定义亮度范围大小的因子f
  /*迭代器的数目是基于被分割的解剖学结构亮度的均匀性之上指定的。均匀性高的区域仅
  需要一对迭代器。带有 ramp 效果的区域，如带有不同一领域的 MRI 图像，就需要更多的迭代
  器。实际上，精心选择乘法因子似乎比迭代器的数目更加重要。然而，务必紧记这个算法毫
  无疑问需要集中于一个稳定的区域。在这个算法上使用过多的迭代器将很有可能使这个区域
  吞没整个图像*/
  confidenceConnected->SetNumberOfIterations( 5 );//迭代器数目（迭代次数）
  /*这个滤波器的输出是一个二值图像，这个二值图像除了分割出的区域外到处都是零值像
  素。区域中的亮度值是由 SetReplaceValue() 方式来选择的*/
  confidenceConnected->SetReplaceValue( 255 );

  /*这个算法的实例化需要用户提供一个种子点。将这个点选在被分割的解剖学结构的典型
  区域是很便捷的。种子区域周围的一个小的邻域将用来计算包含在标准里的最初的平均值和
  标准差。种子是以一种 itk::Index 的形式传递给 SetSeed() 方式的*/
  InternalImageType::IndexType  index;
  //设置种子点坐标位置
  index[0] = atoi( "107" );
  index[1] = atoi( "69" );
  confidenceConnected->SetSeed( index );
  /*种子区域周围邻域最初的大小是使用 SetInitialNeighborhoodRadius() 方式来定义的。这
  个邻域将定义成一个拥有 2r + 1 个像素的 N（二维图像为2，三维为3） 维矩形区域，其中 r 是作为最初邻域范围的传递值*/
  confidenceConnected->SetInitialNeighborhoodRadius( 2 );//设置最初邻域范围的传递值r  即5像素矩阵
  /*writer 上的 Updata() 方法引发了管道的运行。通常在出现错误和抛出异议时, 从一个
  try / catch 模块调用 updata*/
  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }
  
  return EXIT_SUCCESS;
}

                            

                        输入图像                         分割出的白质组织                  分割出的脑室组织                  分割出的灰质组织

注意：灰质部分并未被完全分割，这就是区域生长方法在对被分割的解剖学结构在图像空间上没有相似的统计分布时进行分割的弱点。你可以通过使用不同的迭代次数进行分割，以达到扩展可接受区域的目的。

 

ITK系列目录：

1 ITK图像数据表达之图像

2 ITK图像处理之图像滤波

3 ITK图像处理之图像分割

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