ITK系列29_ ITK分水岭算法对PNG图像进行二维分割

基于分水岭算法的分割

分水岭分割对图像特征使用梯度下降法和沿区域边界分析弱点 (weak points) 来将像素分类为区域。想像在一个有水流动的拓扑地形结构中，水在重力的引导下聚集到一个地势较低的盆地。随着水量的增加，水将流满整个盆地直到水流溢出到另一个盆地，这样就会将一些小盆地吞没形成大的盆地。使用局部的几何结构来形成区域 ( 集水的盆地 ) ，在图像领域中正如使用一些诸如曲率或梯度强度等特征中的局部极值来将像素连接成区域。这种技术不像其他区域分割，它几乎不需要用户定义门限，尤其适合对以不同的特征类型从不同的数据集融合而来的图像进行分割。分水岭技术在那些不仅仅产生单一图像分割的应用中也更加灵活，但是使用一个门限，或者交互式地在一个用户绘图界面的帮助下，从一个可以提取一个单一区域或区域集的 a-priori 中的一个分割层级更合适。

用来实现分水岭方法的两个常见的不同算法：顶部下降和底部上升。顶部下降，梯度下降策略是 ITK 选择的算法，因为我们想考虑这个多元微分操作的输出和问题中的 f 将会有浮动点值。底部上升策略以图像中的局部最小值处的种子开始并在离散的亮度水平上向外向上生长（等同于一个形态学操作上的次序，有时也称为形态学分水岭）。这将通过加强图像上的一系列离散灰度值限制精确度。

接下来的例子阐述了如何使用 itk::WatershedImageFilter 对图像进行预处理和分割。注意：预处理的数据将大大影响到结果的质量。典型地，通过使用边缘保护扩散滤波器对原始图像进行预处理可以得到最好的结果，比如各向异性扩散滤波器或双边图像滤波器。如 9.2.1 小节中介绍的那样，用来作为输入的高度函数，必须通过同对象边缘相关的正的更高的值来创建。对许多应用都合适的一个高度函数是使用被分割图像的梯度大小。使用 itk::VectorGradientMagnitudeAnisotropicDiffusionImageFilter 类来平滑图像并使用itk::VectorGradientMagnitudeImageFilter 类来创建高度函数。我们以包含所有预处理滤波器的头文件和分水岭滤波器的头文件开始。我们使用这些滤波器的向量形式，因为输入数据是一个彩色图像。

把这些滤波器参数转换用来做任何特殊的应用是一个实验和误差的过程。门限参数可以用来控制图像的过分割，提高门限通常将减少计算时间和产生更少更大的区域。转换参数的诀窍是考虑你将要分割的图像中的对象的范围级别。

实例29 ITK分水岭算法对PNG图像进行二维分割

#include <iostream>
#include "itkVectorGradientAnisotropicDiffusionImageFilter.h"
#include "itkVectorGradientMagnitudeImageFilter.h"
#include "itkWatershedImageFilter.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVectorCastImageFilter.h"
#include "itkScalarToRGBPixelFunctor.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 conductanceTerm diffusionIterations lowerThreshold outputScaleLevel gradientMode " << std::endl;
    return EXIT_FAILURE;
    }*/
  /*现在我们声明图像和像素类型用来势力化滤波器。所有的滤波器都需要实数类型的像素
  类型来正常工作。预处理阶段直接使用向量类型的数据而分割时使用浮点型标量数据。使用
  itk::VectorCastImageFilter 来将图像从 RGB 像素类型转换为数字向量类型*/
  typedef itk::RGBPixel< unsigned char >       RGBPixelType;
  typedef itk::Image< RGBPixelType, 2 >        RGBImageType;
  typedef itk::Vector< float, 3 >              VectorPixelType;
  typedef itk::Image< VectorPixelType, 2 >     VectorImageType;
  typedef itk::Image< itk::IdentifierType, 2 > LabeledImageType;
  typedef itk::Image< float, 2 >               ScalarImageType;
  //使用上面创建的类型来声明各种图像域处理滤波器，并最终将它们用在处理过程中
  typedef itk::ImageFileReader< RGBImageType >   FileReaderType;
  typedef itk::VectorCastImageFilter< RGBImageType, VectorImageType >
                                                 CastFilterType;
  typedef itk::VectorGradientAnisotropicDiffusionImageFilter<
                        VectorImageType, VectorImageType >
                                                 DiffusionFilterType;
  typedef itk::VectorGradientMagnitudeImageFilter< VectorImageType >
                                                 GradientMagnitudeFilterType;
  typedef itk::WatershedImageFilter< ScalarImageType >
                                                 WatershedFilterType;
  
  typedef itk::ImageFileWriter<RGBImageType> FileWriterType;

  FileReaderType::Pointer reader = FileReaderType::New();
  reader->SetFileName("VisibleWomanEyeSlice.png");

  CastFilterType::Pointer caster = CastFilterType::New();
  /*接下来我们实例化这些滤波器并设置它们的参数。第一步在图像预处理管道中使用一个
  各向异性扩散滤波器来扩散输入彩色图像。对于这种滤波器类， CFL 条件要求对二维图像的
  time step 不能超过 0.25 ，对三维不能超过 0.125 。迭代器的数量和 conductance term 将从命令行
  得到。参见 6.7.3 小节将得到更多关于 ITK 各向异性扩散滤波器的信息*/
  DiffusionFilterType::Pointer diffusion = DiffusionFilterType::New();
  //设置迭代次数
  diffusion->SetNumberOfIterations( atoi("10") );
  //设置参数 conductance
  diffusion->SetConductanceParameter( atof("2.0") );
  diffusion->SetTimeStep(0.125);
  //对向量类型图像 ITK 梯度大小滤波器可以随意地选择几个参数，这里我们仅允许那些主
  //要成分分析的可用和不可用
  GradientMagnitudeFilterType::Pointer
    gradient = GradientMagnitudeFilterType::New();
  //设置主要组成部分
  gradient->SetUsePrincipleComponents(atoi("on"));
  /*最后我们设置分水岭滤波器。它有两个参数。水平 Level 控制分水岭深度，而门限
  Threshold 控制输入的最低门限值。这两个参数都作为输入图像中最大深度的一个百分比
  (0.0 - 1.0) 来设置的*/
  WatershedFilterType::Pointer watershed = WatershedFilterType::New();
  //设置level
  watershed->SetLevel( atof("0.2") );
  /*分水岭滤波器的输出是一个无符号长整型 lable 图像，其中 lable 表示在一个特殊分割区域
  中一个像素的成员关系。这种形式对视觉是不现实的，所以为了这个例子的目的，我们将把
  它转换为 RGB 像素。 RGB 图像具有可以被保存为 png 文件和使用标准视图软件观看的优点。
  itk::Functor::ScalarToRGBPixelFunctor 类是可以将一个标量值转换成一个 itk::RGBPixel 的一
  个特殊功能对象。将这个算符写进 itk::UnaryFunctorImageFilter 来创建一个图像滤波器，用来
  把标量图像转换成 RGB 图像*/
  //设置阈值
  watershed->SetThreshold( atof("0.01") );
  
  typedef itk::Functor::ScalarToRGBPixelFunctor<unsigned long>
    ColorMapFunctorType;
  typedef itk::UnaryFunctorImageFilter<LabeledImageType,
    RGBImageType, ColorMapFunctorType> ColorMapFilterType;
  ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New();
  
  FileWriterType::Pointer writer = FileWriterType::New();
  //保存的文件名
  writer->SetFileName("VisibleWomanEyeSlice_FSL1.png");
  //这种滤波器连接到一个单一的管道，在每个结尾都有 readers 和 writers
  caster->SetInput(reader->GetOutput());
  diffusion->SetInput(caster->GetOutput());
  gradient->SetInput(diffusion->GetOutput());
  watershed->SetInput(gradient->GetOutput());
  colormapper->SetInput(watershed->GetOutput());
  writer->SetInput(colormapper->GetOutput());

  try
    {
    writer->Update();
    }
  catch (itk::ExceptionObject &e)
    {
    std::cerr << e << std::endl;
    return EXIT_FAILURE;
    }

  return EXIT_SUCCESS;
}

中图是以参数 conductance=2.0、迭代=10、阈值=0.001、level=0.15、主要组成部分=off 生成的。生成右图使用的参数为 conductance=2.0、迭代=10、阈值=0.01、level=0.2、主要组成部分=off生成：

                                                                        

                      初始图像                                          分割图1                                               分割图2

分水岭算法计算复杂性的一个问题是授权许可问题。 ITK 实现的大部分复杂性是在产生层级，这一阶段的处理时间和最初分割中的集水盆地的数目是非线性关系，这就意味着一个图像中包含的信息数量远比图像中的像素数目重要。对一个很大但比较均匀的输入图像进行分割可能比对一个很小但细节很突出的输入图像进行分割所需要的时间更短。

 

ITK系列目录：

1 ITK图像数据表达之图像

2 ITK图像处理之图像滤波

3 ITK图像处理之图像分割

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