图像分割算法详解

图像分割

边缘检测

边缘检测的基本步骤是：

1. 平滑滤波：由于梯度计算容易受到噪声的影响，因此第一步是用滤波去除噪声。但是，降低噪声的平滑能力越强，边界强度的损失就越大。

2. 锐化滤波：为了检测边界，必须确定某点领域中灰度的变化。锐化操作加强了存在有意义的灰度局部变化位置的像素点。

3. 边缘判定：在图像中存在许多梯度不为零的点，但是对于特定应用，不是所有点都有意义。这就要求操作者根据具体情况选择和去除处理点，具体的方法包括二值化处理和过零检测等。

4. 边缘连接：将间断的边缘连接成为有意义的完整边缘，同时去除假边缘。主要方法有：Hough变换。

几个常用的梯度算子：

Robert算子利用局部差分算子寻找边缘，边缘定位精度比较高，但是容易丢失一部分边缘，同时由于图像没有经过平滑处理，因此不具备抑制噪音的能力。该算子对于具有陡峭边缘且含有少量噪音的图像效果比较好。

Sobel算子和Prewitt算子都考虑了邻域信息，相当于对图像先做加权平滑处理，然后再做微分运算，所不同的是平滑部分的权值有些不同，因此对噪声具有一定的抑制能力，但不能完全排除检测结果中出现的虚假边缘。虽然两个算子边缘定位效果不错，但检测出来的边缘容易出现多像素宽度。

高斯–拉普拉斯算子是一个二阶导数，对噪声具有无法接受的敏感度，而且其幅值很产生双边缘，而且边缘方向的不可检测性也是拉普拉斯算子的缺点之一，因此一般不以原始形式用于边缘检测，为了弥补其缺陷，再运用拉普拉斯算子一般先进行高斯低通滤波。

Canny边缘检测算子

上面几种算法都是基于微分方法的边缘检测算法，它们都只有再图像不含噪音或者首先通过平滑去噪的前提下才可以正常使用。在图像检测中，抑制噪声和边缘精确定位是无法同时满足的，一些边缘检测算法通过平滑滤波去噪的同时，也增加了边缘定位的不确定性，而提高边缘检测算子对边缘的敏感性的同时，也提高了对噪声的灵敏度。Canny算子力图在抗噪声干扰和精确定位之间寻求最佳的折中方案。

Canny对边缘检测质量进行分析，提出了3个准则：

1. 信噪比
2. 定位精度
3. 单边缘响应准则

Canny边缘检测的基本思想就是首先对图像选择一定的Gauss滤波器进行平滑滤波，然后采用非极值抑制技术进行处理得到最后的边缘图像。其步骤如下：

1. 用高斯滤波器平滑图像
2. 用一阶偏导的有限差分来计算梯度的幅值和方向
3. 对梯度幅值进行非极大抑制。仅仅得到全局的梯度并不足以确定边缘，为了确定边缘，必须保留局部梯度最大的点，而抑制非极大值，即将非局部极大值点置为零以得到细化的边缘。
4. 用双阈值算法检测和连接边缘。

霍夫变换

霍夫（hough）变换是一个非常重要的检测间断点边界形状的方法。它通过将图像坐标空间变换到参数空间，来实现直线和曲线的拟合。

Hough变换直线检测的MATLAB实现

检测的时候需要三个步骤：

1. 利用hough()函数执行霍夫变换，得到霍夫矩阵
2. 利用houghpeaks()函数在霍夫矩阵中寻找峰值点
3. 利用houghlines()函数在之前两步结果的基础上得到原二值图像的直线信息。

[H,theta,rho] = hough(BW)
[H,theta,rho] = hough(BW,Name,Value,...)

1. BW是边缘检测后的二值图像
2. H是变换得到的Hough矩阵
3. theta与rho为分别对应于Hough矩阵每一列和每一行的theta和rho值组成的向量

peaks = houghpeaks(H,numpeaks)
peaks = houghpeaks(___,Name,Value,...)

1. H是由hough()函数得到的Hough矩阵
2. numpeaks是要寻找的峰值数目，默认值为1
3. peaks是一个Q*2的矩阵，每行的两个元素分别为某一峰值点在Hough矩阵的行，列索引，Q为找到的峰值点的数目

lines = houghlines(BW,theta,rho,peaks)
lines = houghlines(___,Name,Value,...)

1. BW是边缘检测后的二值图像
2. theta与rho是由hough()函数返回的
3. peaks是由houghpeaks()函数返回的
4. lines是一个结构体数组，数组长度是找到的直线条数，而每一个数组元素（直线段结构体）的内部结构

clc
clear all
I  = imread('circuit.tif');

%旋转图像并寻找边缘
rotI = imrotate(I,33,'crop');
BW = edge(rotI,'canny');

%执行hough变换并显示Hough矩阵
[H,T,R] = hough(BW);
imshow(H,[],'XData',T,'YData',R,...
            'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;

%在Hough矩阵中寻找前5个大于Hough矩阵中最大值0.3倍的峰值
P  = houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));
x = T(P(:,2)); y = R(P(:,1));    %由行列索引转换成实际坐标
plot(x,y,'s','color','white');   %在Hough矩阵图像中标出峰值位置

%找到并绘制直线
lines = houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);   %合并距离小于5的线段，丢弃所有长度小于7的直线段
figure, imshow(rotI), hold on
max_len = 0;
for k = 1:length(lines)     %依次标出各条直线
   xy = [lines(k).point1; lines(k).point2];
   plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

   % 绘制线段端点
   plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
   plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

   % 确定最长线段
   len = norm(lines(k).point1 - lines(k).point2);
   if ( len > max_len)
      max_len = len;
      xy_long = xy;
   end
end
%高亮显示最长线段
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');

阈值分割

迭代选择阈值法

迭代式阈值选择方法的基本思想是：开始选择一个阈值作为初始估计值，然后按照某种规则不断更新这一估计值，直到满足给定的条件为止。这个过程关键是选择怎么样的迭代规则。一个好的迭代规则必须能够快速收敛，又能够在每一个迭代过程中产生的由于上一次迭代的结果。

最小均方误差法

这种方法通常以图像中的灰度为模式特征，假设各模式的灰度是独立分布的随机变量，并假设图像中待分割的模式服从一定的概率分布。一般来说，采用的是正态分布。

首先假设一幅图像仅包含两个主要灰度区域，前景和背景。令z表示灰度值，p(z)表示灰度值概率密度函数的估计值。假设概率密度函数中的参数一个对应于背景的灰度值，一个对应于图像的前景即对象的灰度值。

即图像中的像素只能是属于前景或者是背景，没有第三种情况，现在选定一个阈值T，将图像上的像素进行分类。采用最小均方误差发的目的是选择T时，使对一个给定像素进行使前景还是背景的分类时出错的概率最小，

最大类间方差法

在对图像进行阈值分割时，选定的分割阈值应该使前景区域的平均灰度，背景区域的平均灰度与整幅图的平均灰度之间差别最大，这种差异用区域的方差来表示。由此otsu在1978年提出的最大方差法。该算法在判决分析最小二乘法原理的基础上推到出来的。

区域分割

区域生长是根据事先定义的准则将像素或者子区域聚合成更大区域的过程。其基本思想是从一组生长点开始（生长点可以是单个像素，也可以为某个区域），将与该点该生长点性质相似的相邻像素或者区域与生长点合并，形成新的生长点，重复此过程直到不再生长为止。生长点和相邻区域的相似性判据可以是灰度值，纹理，颜色等多种图像信息。

MATLAB的实现方式：

function J=regionGrow(I)

if isinteger(I)
    I=im2double(I);
end

figure,imshow(I),title('原始图像')
[M,N]=size(I);
[y,x]=getpts;      %获得区域生长起始点
x1=round(x);       %横坐标取整
y1=round(y);       %纵坐标取整
seed=I(x1,y1);     %将生长起始点灰度值存入seed中
J=zeros(M,N);      %做一个全零与原图像等大的图像矩阵J，作为输出图像矩阵
J(x1,y1)=1;        %将J中与所取点相对应为止的点设置为白
sum=seed;          %存储符合区域生长条件的点的灰度值的和
suit=1;            %存储符合区域生长条件的点的个数
count=1;           %记录每次判断一点周围八点符合条件的新点的数目
threshold=0.15;    %阈值，注意需要和double类型存储的图像相符合

while count>0
    s=0;           %记录判断一点周围八点时，符合条件的新点的灰度值之和
    count=0;
    for i=1:M
        for j=1:N
            if J(i,j)==1
                if (i-1)>0 & (i+1)<(M+1) & (j-1)>0 & (j+1)<(N+1)   %判断此点是否为图像边界上的点
                    for u=-1:1
                        for v=-1:1
                            if J(i+u,j+v)==0 & abs(I(i+u,j+v)-seed)<=threshold & 1/(1+1/15*abs(I(i+u,j+v)-seed))>0.8
                                J(i+u,j+v)=1;
                                %判断是否尚未标记，并且为符合阈值条件的点
                                %符合以上两个条件即将其在J中与之位置对应的点设置为白
                                count=count+1;
                                s=s+I(i+u,j+v);      %此点的灰度加入s中
                            end
                        end
                    end
                end
            end
        end
    end
    suit=suit+count;       %将n加入符合点数计数器中
    sum=sum+s;             %将s加入符合点的灰度值总合中
    seed=sum/suit;         %计算新的灰度平均值
end

           

调用的方法是：

I=imread('coins.png');
if ndims(I)>=3
    I=rgb2gray(I);
end
J=regionGrow(I);
figure;
imshow(J);

上面程序运行之后，会弹出一个包含原图像的窗口，用户可以用鼠标在其中选取一个种子并按下“Enter”键之后，

区域分裂与合并及其MATLAB实现

区域生长是从一组生长点开始的，另以中国方法是在开始的时候将图像分割为一系列任意不相交的区域，然后将它们合并或者拆分以满足限制条件，这就是区域分裂与合并。通过分裂，可以将不同特征区域分离开，而通过合并，将相同特征的区域合并起来。

区域分裂的MATLAB实现

在MATLAB中，和区域分裂相关的3个主要函数是：qtdecomp(),qtgetblk(),qtsetblk()。

S = qtdecomp(I)
S = qtdecomp(I,threshold)
S = qtdecomp(I,threshold,mindim)
S = qtdecomp(I,threshold,[mindim maxdim])
S = qtdecomp(I,fun)

1. I是输入的灰度图像
2. threshold是分割成的子块中允许的阈值，默认值为0，如果子块中最大的元素和最小元素的差值小于该阈值就认为满足一致性条件；对于double型矩阵，threshold将直接作为阈值；而对于uinit8和uinit16类型矩阵，threshold将被乘以255和65535以作为实际阈值；而对图像而言，threshold的取值范围是0到1；
3. [mindim maxdim]是尺度阈值；mindim参数可以屏蔽函数对尺度上小于mindim的子块的处理，而不论这个子块是否满足一致性条件，如果参数形式是[mindim maxdim]，则表示不产生小于mindim尺度的子块，也不保留大于maxdim尺度的子块，此时maxdim/mindim必须是2的整数次幂。
4. S是一个稀疏矩阵，在每个子块的左上角给出子块的大小

qtgetblk()函数。在得到稀疏矩阵S之后，利用IPT函数qtgetblk()可进一步获得四叉树分解后所有指定大小的子块像素及位置信息。

[vals,r,c] = qtgetblk(I,S,dim)
[vals,idx] = qtgetblk(I,S,dim)

1. I是输入的灰度图像
2. 稀疏矩阵S是I经过qtdecomp()函数处理的输出结果
3. dim是指定的子块大小
4. vals是dim*dim k的三位矩阵，包含I 中所有符合条件的子块数据，其中k为符合条件的dimdim大小的子块的个数，vals(:,:,i)表示符合条件的第i个子块的内容；
5. r和c均为列向量，分别表示图像I中符合条件子块左上角的纵坐标和横坐标。

qtsetblk()函数

在将图像划分为子块后，还需要使用函数qtsetblk()将四叉树分解所得到的子块中符合条件的部分全部替换为指定的子块。函数的语法如下：

J = qtsetblk(I,S,dim,vals)

matlab实现

I1=imread('rice.png');
imshow(I1);
%选择阈值为0.2，对原始图像进行四叉树分解
S=qtdecomp(I1,0.2);
%原始的稀疏矩阵转换为普通矩阵，使用full函数
S2=full(S);
figure
imshow(S2);
%记录子块数目的列向量
ct=zeros(6,1);
%分别获取不同大小块的信息，子块内容保存在三维数组vals1-vals16中，子块数目保存在ct向量中
for i=1:6
    [vals{i},r,c]=qtgetblk(I1,S2,2^(i-1));
    ct(i)=size(vals{i},3);
end