何颖+张星阳+赵金龙
摘要：将目标图像利用具体的算法操作将待处理图像分割成几个特定的、分别具有独特性质的区域，并把目标区域和背景分离出来，分割结果的好坏直接关系着图像的进一步分析。通过区域生长法对右心室MRI图像进行分割算法研究，区域生长法要想区域分割算法能够成功的三个重要因素是，选择合理的种子点像素，设计针对特定图像合理的生长准则，制定合理的生长停止的条件。通过采用置信连接与孤立连接分割法的对比，实验验证了分割的效果比对。
关键词：右心室分割；区域生长法；医学图像；置信连接法
中图分类号：TP391.41 文献标识码：A 文章编号：1007-9416(2017)10-0115-03
1 引言
图像分割的概念主要是指：将目标图像利用具体的算法操作将待处理图像分割成几个特定的、分别具有独特性质的区域，并把目标区域和背景分离出来的分割过程，分割结果的好坏直接关系着图像的进一步分析[1]，也是在整个图像处理过程中比较困难和难以实现的一个步骤。
对于医学图像来说，其图像自身包含很多的信息，可是往往只有其中的一部分是临床诊断所关心的，往往称这一部分为感兴趣区域(Region Of Interest， ROI)，而相比较这一部分其余的部分则是对临床诊断来说并不感兴趣的区域。怎样在这些特定区域得到对临床研究有用的图像成分和解剖结构，是图像分割技术的核心问题。
近几年来，有越来越多的学者和医学研究者在研究和临床案例中发现了右心室对心血管疾病的提前预测和后续诊断治疗方面的重要作用，右心室的医学意义也越来越多地受到认可，对于右心室的分割也愈来愈受到大家的重视，然而另一方面，右心室和左心室相比有更复杂的特性，比如说高变异性、壁薄，边界区域不够明显以及和周围组织之间难以进行区分等等，这一系列特征使得右心室分割在实践操作中又增加了一定的难度。近几十年来医学图像处理中提出了一些有价值的方法，比方说，基于统计形状模型、三维模型的方法[2]和基于形变模型的方法[3]取得了不错的效果。
2 区域生长算法
区域生长算法[4]的概念是指将具有满足某种相似特性的各个像素点被划分到同一个区域来实现分割。首先针对整幅图像的多个待分割区域，在每一个待分割的区域内选择一个种子点作为区域生长的起点，然后需要按照使自己目标最优化的生长准则或者把它周围与该个像素点特征相近或者相似的像素合并到预先设置的种子像素所在的区域中，之后把这些合并完成的新的像素作为种子区域按照上述方法继续进行生长，直到遍历完整幅图像，以至于整幅图像中没有满足预先设定的条件或者准则的像素可被合并到种子区域内时，结束整个区域生长分割过程。
要想区域分割算法能够成功实施三个重要因素必不可少：(1)选择合理的种子点像素；(2)设计针对特定图像合理的生长准则；(3)制定合理的生长停止的条件。关于种子像素的选取，采取的手段既可以是通过人工交互的方法也可以通过采用全自动的方法。前者的方法不仅简单并且实用，可用于绝大多数的物体和图像，使用性更加广泛；后者方法则针对性比较单一，往往面对的是具体的问题，比如说，运用红外线来检测目标的时候，往往选择亮度最高的像素作为种子像素来使用。区域生长法[5]的算法有置信连接分割方法和孤立点连接分割方法等。
右心室的磁共振图像如图1所示，第一排显示的为右心室在心底时期获取到的磁共振图像，第二排显示的为右心室在心尖时期获取到的磁共振图像，第一列的显示状态则是处于舒张末期，第二列显示的恰恰是处于收缩末期。由图1中可以总结出：不同病人甚至是同一个病人的右心室在不同的时期其获取的右心室的磁共振图像形状上也有很大的差别，而这一系列因素增加了对右心室得到有效分割的难度。
2.1 置信连接分割法
置信连接分割法的概念是：置信连接分割法实际上本身也是区域生长分割方法的一种，根据当前区域简单统计后的结果，并把种子点的某个邻域范围内的的所有像素点的灰度的平均值和标准差计算出来，之后按照提前设置好的因子和标准差的乘积当作范围，则邻域内像素点的灰度值在这个计算出的范围之内的情况下，则把这个像素点将被划分到这个区域之内。接下来，按照上次划分好的区域，再次计算上次划分好的区域内的所有像素点灰度的平均值和标准差，按照上面方法重新计算出灰度范围，再次检验当前区域的邻域的像素点的灰度值是否在新的灰度范围之内。直至达到迭代停止的条件：某一次迭代时并没有新元素被加入到这个区域或者是达到了迭代器设置的迭代次数的上限。表1显示了用置信连接分割方法分割右心室磁共振图像需要设置的参数、得到结果的Dice重合率和分割结果，图像的标号如图2(a)、图2(b)和图2(c)。
分割结果如图2所示。
2.2 孤立连接分割法
孤立连接法也属于区域生长算法的一种，在分割过程中一般需要提前设置好三个参数，分别是：两个种子点和一个下门限值，生长区域被设置在与第一个种子点相连却孤立于第二个种子点的区域。為了检验这个算法和置信连接分割法的效果，本算法，仍然选择相同三幅图像进行操作，首先进行种子点的选取，然后设定合适的下门限灰度值分割三幅图像。表2显示了使用孤立连接法对上述三幅图像进行操作时需要用到的实验参数以及用孤立连接算法在参数下得到的分割图像的Dice重合率和分割结果，图像的标号如图3(a)、图3(b)和图3(c)。
分割结果如图3所示。
3 结语
论文对右心室磁共振图像涉及的三幅图像进行简单分割。通过置信连接法和孤立连接法的分割结果和Dice重合率发现，以上两种方法的分割结果在个别磁共振图像上的分割达到了比较好的效果，并且要求操作者对于区域生长分割算法需要手动选取种子点等操作，另外，区域生长法对于右心室磁共振图像的分割结果中不足之处在于：使用区域生长法进行分割时，分割效果容易受到图像自身噪声的影响，容易造成区域空洞。
参考文献
[1]Caudron J， Fares J， Vivier P H， et al. Diagnostic accuracy and variability of three semi-quantitative methods for assessing right ventricular systolic function from cardiac MRI in patients with acquired heart disease[J]， European radiology， 2011，21(10)：2111-2120.
[2]Petitjean C， Dacher J N. A review of segmentation methods in short axis cardiac MR images[J]， Medical Image Analysis， 2011，15(2)：169-184.
[3]Lapp R M， Lorenzo-Valdes M， Rueckert D. 3D/4D cardiac segmentation using active appearance models， non-rigid registration， and the insight toolkit[J]， MICCAI，2004： 419-426.
[4]朱峰.医学图像自动分割若干关键技术研究[D].苏州：江苏大学，2006.
[5]周振环.医学图像分割与配准[M].成都：电子科技大学出版社，2011，136-183.endprint

1.介绍

本项目使用ubuntu16.04，ubuntu的其他版本应该也是可以的；

点云库PCL某某这本书其实是以官方的英文手册翻译成的，英文链接见以下地址：

https://pcl.readthedocs.io/projects/tutorials/en/latest/region_growing_segmentation.html#region-growing-segmentation

但是想做相关的算法开发，只看书中的内容还是不够明白，只能找来官方文档读了，以下基本是翻译再加一些我自己的理解，如果有不清楚的还是建议自己去啃英文说明。

使用pcl:: regiong类实现区域增长算法。该算法的目的是在光滑性约束条件下合并相近的点，实现区域的分割，适用于不同平/曲面的分割。该算法的输出是不同曲面单元的集合，其中每个单元可以认为是一个光滑表面的点的总成。该算法基于点之间的法线角度比较完成划分。

2.算法实现

因为生长点是根据曲率进行寻找的，选择曲率最小的点，所以首先根据曲率将点进行分类。由于选择的曲率是最小的点，所以生长点附近一般均为平面或者相对较平坦的面。

将生长点添加进种子的集合。然后对种子里的每个生长点进行生长运算。

种子点的运算过程如下：

测试每个临近点与种子点法向量的夹角，如果角度小于临界值，则添加进当前的种子区域；
	
测试每个临近点的曲率，如果曲率小于临界值，则将该点添加进当前区域的种子列表；
	
第三步的意思其实是在找到曲率小于临界值的点后添加进种子列表，原来的种子点就要删除了，但是这一个地方有一个问题，是否需要原有的种子周围临近点全部搜索完毕，这个应该是要搜索完毕的；但是如果发现多个曲率小于临界值的点，是否需要都添加？我的理解是都需要添加；
注：曲面的一阶微分的分量包括法向量，曲面的二阶微分的分量包括曲率；从某种意义上说法向量的夹角是类似与曲率的。但是为什么区域增长算法通过夹角生成区域，而通过曲率生成种子列表？这个后面会有描述。

待种子列表中为空时，说明该区域生长完毕；

这里贴出算法伪代码：





 

 

关于伪代码这里就不解释了，只要花时间认真去看，相信是可以理解的。其中的\是去除的意思。

代码如下：

#include <iostream>
#include <vector>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/search/search.h>
#include <pcl/search/kdtree.h>
#include <pcl/features/normal_3d.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/filters/passthrough.h>
#include <pcl/segmentation/region_growing.h>

int
main (int argc, char** argv)
{
//加载pcd文件
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("region_growing_tutorial.pcd", *cloud) == -1)
  {
    std::cout << "Cloud reading failed." << std::endl;
    return (-1);
  }
//通过NormalEstimation类计算法线
  pcl::search::Search<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
  pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
  normal_estimator.setSearchMethod (tree);
  normal_estimator.setInputCloud (cloud);
  normal_estimator.setKSearch (50);
  normal_estimator.compute (*normals);
//在z轴方向滤波，剔除0~1.0范围之外的点
// PassThrough是滤波类，得到的indices是各点的索引值
  pcl::IndicesPtr indices (new std::vector <int>);
  pcl::PassThrough<pcl::PointXYZ> pass;
  pass.setInputCloud (cloud);
  pass.setFilterFieldName ("z");
  pass.setFilterLimits (0.0, 1.0);
  pass.filter (*indices);
//设置一个簇中包含点的最大值门限与最小值门限，通过区域生长算法生成的所有的簇，小于50个点或者大//于1000000个点都被丢弃，不包含在簇内
  pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
  reg.setMinClusterSize (50);
  reg.setMaxClusterSize (1000000);
//算法需要k近邻搜索，搜索近邻点的数量限制为30
  reg.setSearchMethod (tree);
  reg.setNumberOfNeighbours (30);
  reg.setInputCloud (cloud);
  //reg.setIndices (indices);
  reg.setInputNormals (normals);
//此处两个初始化的值是最关键的，直接影响分割结果的好坏；
// setSmoothnessThreshold是法线的夹角门限
// setCurvatureThreshold是曲率门限，对于满足发现夹角的点，如果也满足曲率门限
//则将会被添加进种子列表内；
  reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
  reg.setCurvatureThreshold (1.0);
//下面执行区域生长算法，返回一个簇的数组
  std::vector <pcl::PointIndices> clusters;
  reg.extract (clusters);
//显示一些信息
  std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;
  std::cout << "First cluster has " << clusters[0].indices.size () << " points." << std::endl;
  std::cout << "These are the indices of the points of the initial" <<
    std::endl << "cloud that belong to the first cluster:" << std::endl;
  int counter = 0;
  while (counter < clusters[0].indices.size ())
  {
    std::cout << clusters[0].indices[counter] << ", ";
    counter++;
    if (counter % 10 == 0)
      std::cout << std::endl;
  }
  std::cout << std::endl;
//显示分割后的图像
  pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
  pcl::visualization::CloudViewer viewer ("Cluster viewer");
  viewer.showCloud(colored_cloud);
  while (!viewer.wasStopped ())
  {
  }
//
  return (0);
}


需要将region_growing_tutorial.pcd替换成你想分析的pcd文件，或者自己去网站下载相应的pcd文件；

3.编译与执行

CMakeList.txt内容如下：

cmake_minimum_required(VERSION 2.8 FATAL_ERROR)

project(region_growing_segmentation)

find_package(PCL 1.5 REQUIRED)

include_directories(${PCL_INCLUDE_DIRS})
link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})

add_executable (region_growing_segmentation region_growing_segmentation.cpp)
target_link_libraries (region_growing_segmentation ${PCL_LIBRARIES})


将CMakeList.txt与cpp文件放在同一个文件夹，然后依次执行

cmake CMakeList.txt

make

./region_growing_segmentation

执行之后就显示出了分割后的图，不同区域通过不同颜色分割，其中的红色点是为分区的点，

  reg.setMinClusterSize (50);

  reg.setMaxClusterSize (1000000);

红色点区域上述规划的小于50或者大于1000000的情况

4.问题总结

关于依据曲率选择种子的问题

根据官方的文档可以看出，是在通过法线夹角筛选后的点中选择满足曲率要求的种子点的。同时，种子点也是加进种子列表的，其原来的点并没有删除，只有当原来的点的规定半径区域的点运算完毕后才会将其删除，这一点官方文档描述的不是很清楚，但是伪代码写的是对的。

关于初始种子的选取问题

另一篇文章有描述，是将所有点按照曲率值进行排序，找到最小曲率点，将其添加到种子点集，链接如下：

https://cloud.tencent.com/developer/article/1475904

 

这一篇是补充前一篇的。相比前一篇的第二种实现方法，这里使用了均值、标准差随迭代更新的策略，且根据区域的性质，对三通道取不同的权值。修正后的程序为：

clear all;close all;clc;
img=imread('parrot.png');
subplot(3,1,1),
imshow(img);
[m,n,~]=size(img);
seeds=[400,200,1;230,600,2;440,670,3];  %第3列为label
labelNum=3;
newSeeds=seeds;
hold on;
plot(seeds(:,2)',seeds(:,1)','gs','linewidth',1);
title('原始图像及种子位置');
flag=1;
ycc=double(rgb2ycbcr(img));
Y=ycc(:,:,1);
Cb=ycc(:,:,2);
Cr=ycc(:,:,3);
img=double(img);
R=img(:,:,1);
G=img(:,:,2);
B=img(:,:,3);
mask=zeros(m,n);
for i=1:size(seeds,1)
    mask(seeds(i,1),seeds(i,2))=1;
end

avgs=zeros(labelNum,3);%初始均值
stds=repmat([30,10,10],[size(seeds,1),1]); %初始标准差
w=[5,5,5;3,2.2,2.8;3,3,3.5];  %各区域三通道权重

while(flag)
    temp=[];
    
    for i=1:labelNum
        index=find(seeds(:,3)==i);
        selected=seeds(index,1:2);
        indexS=m*(selected(:,2)-1)+selected(:,1);
        yS=Y(indexS);
        CbS=Cb(indexS);
        CrS=Cr(indexS);
        rgbSamp=[yS,CbS,CrS];
        avgs(i,:)=sum(rgbSamp,1)/length(index);
        if size(rgbSamp,1)>1
            stds(i,:)=std(rgbSamp,1);
        end
    end
    stds=max(5,stds);
    for i=1:size(newSeeds,1)
        rowS=newSeeds(i,1);
        colS=newSeeds(i,2);
        label=newSeeds(i,3);
        rlow=max(1,rowS-1);
        rhigh=min(m,rowS+1);
        clow=max(1,colS-1);
        chigh=min(n,colS+1);
        for row=rlow:rhigh
            for col=clow:chigh
                if mask(row,col)==1  %已经添加过的排除掉
                    continue;
                end
                diffycc=reshape(ycc(row,col,:),1,size(ycc,3))-avgs(label,:);
                judge=abs(diffycc)<w(label,:).*stds(label,:);
                if sum(judge)==3 % R、G、B、Y均在3倍均方根范围内
                    temp=[temp;row,col,label];
                    mask(row,col)=1;
                end
            end
        end
    end
    newSeeds=temp;
    seeds=[seeds;newSeeds];
    if isempty(temp)
        flag=0;
    end
end
imlabel=zeros(m,n);
for i=1:size(seeds,1)
    imlabel(seeds(i,1),seeds(i,2))=seeds(i,3);
end
mask=cat(3,mask,mask,mask);
sgImg=mask.*img;
imlabel=0.5*mat2gray(imlabel)+0.5;
subplot(3,1,2);
imshow(uint8(sgImg));
title('分割结果');
subplot(3,1,3);
imshow(imlabel);
title('种子点生长获得的掩模');


分割结果：




可以看到，这种方法比第一种方法获得了更好的分割效果。 

 

转载：      有梦想的田园犬 https://blog.csdn.net/AmbitiousRuralDog/article/details/80267519