数字图像处理-图像分割

单色图像的分割算法通常基于图像亮度值的两个基本特性：不连续性和相似性。在第一种类别中，处理方法是基于亮度的突变来分割一幅图像，如图像中的边缘。在第二种类别中，主要方法是根据事先定义的准则把图像分割成相似的区域。

点、线和边缘检测

用于查找不连续的最常用方法是对整幅图像运行一人掩模。对于大小为3*3的掩模来说，该过程将计算系数和由掩模覆盖区域所包含的灰度级的乘积之和。
R = w 1 z 1 + w 2 z 2 + ⋅ ⋅ ⋅ + w 9 z 9 R = w~1~z~1~ + w~2~z~2~ + ··· + w~9~z~9~ R=w 1 z 1 +w 2 z 2 +⋅⋅⋅+w 9 z 9

点检测

嵌在常数区域（或图像中亮度基本不变的区域）中的孤立点的检测，在原理上相当直接。若
∣ R ∣ > = T |R| >= T ∣R∣>=T
则我们就可以说在掩模的中心位置已检测出了一个孤立点。其中T是一个非负阈值。
点检测在MATLAB中可以用函数imfilter来实现。最重要的需求是当掩模的中心位于一个孤立点时，掩模的相应必须最强，而在亮度不变的区域中相应为零。
若T已给出，则如下命令可以实现点检测方法

 >>g = abs(imfilter(double(f), w)) >= T

其中，f是输入图像；w是中心为8，其余为-1的3*3掩模；g是结果图像。
若T没有指定，我们可以使用g的最大值（T = max(g(: ))，然后再g中查找所有满足g>=T的点。很明显，g中不可能有比T更大的点。
点检测的另一种方法是在m *n的所有邻域中寻找一些点，这些点的最大值和最小值之差超过了T。这种方法可以使用第三章中的非线性滤波函数ordfilt2来实现

>>g = imsubtract(ordfilt2(f, m *n, ones(m, n)), ordfilt2(f, 1, ones(m, n)));
>>g = g >= T;

线检测

使用edge函数的边缘检测

到目前为止边缘检测最通用的方法是检测亮度值的不连续性。这样的不连续性是用一阶和二阶导数检测的。
二维函数f(x,y)的梯度定义为向量
梯度向量的基本性质是它指向f在(x,y)处的最大变化率方向。
在图像处理中，二阶导数通常用拉普拉斯算子来计算，但是拉普拉斯算子自身很少被直接用于做边缘检测，因为二阶导数对噪声具有无法接受的敏感性，它的幅值会产生双边缘，而且它不能检测边缘的方向。然而当与其他边缘检测技术结合使用，拉普拉斯算子是一种有效的补充方法。例如它的双边缘使得它不适合直接用于边缘检测，但是该性质可以用于边缘定位。
边缘检测的基本意图是使用如下两个准则之一在图像中找到亮度快速变化的地方：
1.找到亮度的一阶导数在幅值上比指定的阈值大的地方。
2.找到亮度的二阶导数有零交叉的地方。
IPT函数edge基于上面的讨论提供了几个导数估计器。对一些估计器来说，是有可能指定边缘检测器是否对水平和/或垂直边缘敏感的。

[g, t] = edge(f, 'method', parameters)

其中f是输入图像，method是下表的一种方法，parameters是参数。在输出中g是一个逻辑数组，其值如下决定：在f中检测到边缘的位置为1，在其他位置为0.参数t是可选的，它给出egde使用的阈值，以确定哪个梯度值足够大到可以称为边缘点。
几种检测器
Sobel检测器
其调用语法为

[g, t] = edge(f, 'sobel', T, dir)

其中f是输入图像，T是一个指定的阈值，dir指定检测边缘的首选方向：‘horizontal’、‘vertical’或’both’（默认）。输出参数是可选的，它是函数edge所用的阈值。若指定了T的值，则t=T。若T未被赋值，则函数edge会令t等于它自动确定的一个阈值，然后用于边缘检测。

Prewitt检测器
这种检测器使用上图的掩模来数字化地近似一阶导数G_x和G_y。基本语法同Sobel检测器。Prewitt检测器在计算上比Sobel检测器要简单一些，但容易产生一些噪声（系数为2的Sobel检测器可以提供一种平滑作用）。

Roberts检测器
语法同Sobel检测器。最古老的边缘检测器之一。

Laplacian of Gaussian(LoG)检测器

高斯函数h ( r )是一个平滑函数，若与图像卷积，会使图像变模糊。模糊程度由σ的值决定。对r求二阶导数即是该函数的拉普拉斯算子

因为二阶导数是线性操作，所以使用LoG与一幅图像卷积，和先用平滑函数与该图像卷积，再计算卷积结果的拉普拉斯算子，结果是相同的。这是LoG检测器的关键概念。使用LoG卷积图像有两个效果：它平滑了图像（因而降低了噪声），并计算拉普拉斯算子，进而产生一幅双边缘图像。然后通过查找双边缘之间的零交叉来定位边缘。
LoG检测器的基本语法为

[g, t] = edge(f, 'log', T, sigma)

其中sigma是标准差，其他参数与前面讨论相同。sigma的默认值为2。和之前一样，函数edge会忽略不强于T的任何边缘，若T未指定，则egde会自动选择一个值。将T设为0会产生封闭的轮廓，这是LoG方法的一个常见特征。

零交叉检测器
这种检测器基于与LoG相同的概念，但是卷积是使用一个特定的滤波函数H来完成的。

[g, t] = edge(f, 'zerocross', T, H)

其他参数同LoG。

Canny检测器
该检测器是函数edge中最强大的边缘检测器。方法总结如下
Canny检测器的语法为

[g, t] = edge(f, 'canny', T, sigma)

其中T是一个向量，T = [T1,T2],它包含前述过程的步骤3中解释过的两个阈值；sigma是平滑滤波器的额标准差。如果t包含在输出参数中，那么它是一个包含该算法所用的两个阈值的二元向量。sigma的默认值为1。语法中其余参数和其他方法中的解释一样，包含未指定T值时阈值的自动计算。

使用Sobel检测器提取边缘

>> f = imread('\empire.jpg');
>> f = rgb2gray(f);
>> [gv, t] = edge(f,'sobel','vertical');这里仅取垂直方向边缘
>> imshow(gv)
>> tt =0.1132

结果

使用霍夫变换的线检测

理想情况下，之前介绍的方法应能产生位于边缘上的像素。实际上，得到的像素通常并不能完全表征边缘，因为存在噪声、不均匀照明引起的边缘断裂，以及引入杂散灰度不连续的其他效应。**因此，在执行边缘检测算法后，通常会使用连接过程来把边缘像素组装成有意义的边缘。**连接图像中线段的一种方法时霍夫变换。
这里有详解

函数hough
函数hough的默认语法为

[H, theta, rho] = hough(f)

完整语法为

[H, theta, rho] = hough(f, 'ThetaRes', val1, 'RhoRes', val2)

其中，H是霍夫变换矩阵，theta（单位为度）和rho是ρ和θ值的向量，霍夫变换是在这些值上生成的。输入f是一幅二值图像；val1是值在0到90之间的一个标量，它指定了沿θ轴的霍夫变换容器（默认为1），val2是范围0<val2<hypot(size(I,1), size(I,2))内的一个实标量，它指定了沿ρ轴的霍夫变换容器的间隔。

使用霍夫变换做峰值检测

使用霍夫变换进行线检测和连接的第一步是峰值检测，即找到具有高计数的累加器单元。因为霍夫变换参数空间中的量化和典型图像中的边缘都不是很完美的事实，霍夫变换的峰值通常会出现在多个霍夫变换单元中。克服这个问题的一个策略如下：
1.找到包含有最大的霍夫变换单元并记下它的位置。
2.把第一步中找到的最大值点的邻域中的霍夫变换单元设为零。
3.重复该步骤，直到找到需要的峰值数，或者达到一个指定的阈值时为止。

函数houghpeaks可以用于实现该策略。其默认语法为

peaks = houghpeaks(H, NumPeaks)

或

peaks = houghpeaks(H, NumPeaks, 'Threshold', val1, 'NHoodSize', val2)

其中NumPeaks是指定数量的峰值。peaks是容纳这些峰值的行列坐标的一个Q*2矩阵；Q的范围从0到NumPeaks。H是霍夫变换矩阵。参数val1是一个非负的标量，它指定H中的哪些值被认为是峰值；val1可以从0到Inf变化，默认值是0.5 *max(H(: ))。参数val2是一个奇整数的两元素向量，它指定围绕峰值的邻域的大小。识别出峰值后，该邻域中的元素被置零。默认值是个两元素向量，该向量大于等于size(H)/50的最小奇数值组成。

使用霍夫变换做线检测和连接

一旦在霍夫变换中识别出一组候选的峰值后，剩下的工作就是确定是否存在与这些峰值相关联的有意义线段，以及这些线段的起点和终点。
函数houghlines采用如下策略
1.将像素位置旋转90-θ，以便它们大概位于一条垂直线上。
2.按旋转的x值来对这些像素位置分排序。
3.使用函数diff找到裂口。忽视小裂口；这将合并被小空白分离的相邻线段。
4.返回比最小值长的线段信息。
该函数的默认语法为

lines = houghlines(f, theta, rho, peaks)

或

lines = houghlines(f, theta, rho, peaks, 'FillGap', val1, 'MinLength', val2)

其中theta和rho是函数hough的输出，peaks是函数houghpeaks的输出。输出lines是一个结构数组，该数组的长度等于所找到的线段数。该结构的每个元素识别一条线，并有如下字段：
1.point1，一个两元素向量[r1,c1]，它指定线段终点的行、列坐标。
2.point2，一个两元素向量[r2,c2]，它指定线段其他终点的行、列坐标。
3.theata，与线相关的霍夫变换容器的角度（单位为度）。
4.rho，与线相关的霍夫变换容器的ρ轴的位置。
其他参数如下：val1是一个正的标量，它指定与相同霍夫变换容器相关联的两条线段间的距离。当两条线段间的距离小于指定值时，函数houghlines将把这两条线段聚合为一条线段（默认距离为20像素）。参数val2是一个正的标量，他指定聚合后的线是保留还是丢弃。比val2指定的值要短的线丢弃。

使用霍夫变换做检测和连接
首先使用Canny检测器提取图像边缘

>> f = imread('\empire.jpg');
>> f = rgb2gray(f);
>> [g_best,t] = edge(f,'canny',[0.05,0.15],1.5);
>> figure,imshow(g_best)

结果如下
接着使用hough函数和houghpeaks函数寻找线段和5个峰值，并标记出来

>> [H, theta,rho] = hough(g_best,'ThetaResolution',0.2);%霍夫变换
>> imshow(H,[],'XData',theta,'YData',rho,'InitialMagnification','fit')
>> axis on, axis normal
>> xlabel('\theta'),ylabel('\rho')
>> peaks = houghpeaks(H,5);%霍夫峰值检测，取5个峰值
>> hold on
>> plot(theta(peaks(:,2)),rho(peaks(:,1)),'linestyle','none','marker','s','color','w')

结果如下
最后使用houghlines函数连接线段，并在原图上标粗显示出来。

>> lines = houghlines(g_best,theta,rho,peaks);%霍夫线检测和连接
>> figure, imshow(g_best), hold on
>> for k = 1:length(lines)
xy = [lines(k).point1;lines(k).point2];
plot(xy(:,1),xy(:,2), 'LineWidth', 4, 'color', [.8 .8 .8]);
end

结果如下

阈值处理

灰度阈值处理的成功与分隔直方图模式的峰谷的宽度和深度直接相关。因此影响峰谷性质的关键因素是：1.峰间的可分性（峰分得越开，分离模式的机会越好）；2.图像中的噪声内容（模式随噪声增加而加宽）；3.物体和背景的相对大小；4光源的额均匀性；5图像反射特性的均匀性。

全局阈值处理

选取阈值的一种方法是图像直方图的视觉检测。另一种方法是反复实验来选区不同的阈值，直到观察者认为产生了较好的结果为止。
通常首选的方法是使用一种能基于图像自动地选择阈值的算法。下面的迭代过程就是这样一种方法：
1.为全局阈值选择一个初始估计T。
2.使用T分割图像。这会产生两组像素：所有灰度值大于T的像素组成的G₁，由所有灰度值小于等于T的像素组成的G₂。
3.分别计算区域G₁和G₂中像素的平均灰度值m₁和m₂。
4.计算一个新的阈值：
T = 1 / 2 ( m 1 + m 2 ) T = 1/2(m~1~ + m~2~) T=1/2(m 1 +m 2 )
5.重复步骤2到步骤4，直到后续迭代中T的差小于一个预定的ΔT为止。
6.使用函数im2bw分割图像：

g = im2bw(f, T/den)

其中den是一个整数（例如，对一幅8比特图像，den是255），它把比值T/den的最大值标定为1。
若初始阈值在图像中最小和最大灰度之间选择（平均图像灰度就是T的一个较好初始选择），那么算法会在有限步数内收敛。

使用Otsu方法进行全局阈值处理

类间方差最大化的思想是方差越大，越接近正确分割图像的阈值。注意，这种最佳测度完全基于直接由图像直方图得到的参数。此外，因为k是区间[0,L-1]内的一个整数，所以找到找到最大的方差非常简单：只需要选L个可能k值中的一个，并在每步中计算方差。然乎选择给出最大方差值的k。这个k就是最佳阈值。如果最大值不唯一，那么所用的阈值就是所找到的所有最佳k值的平均值。
工具箱函数graythrash可以计算Otsu阈值。语法为

[T, SM] = graythresh(f)

其中f是输入图像，T是产生的阈值，它被归一化到区间[0,1]，因此在使用前应缩放到合适的范围，SM是可分性测度，其值越高灰度分成两类的可分性就越高。

使用图像平滑改进全局阈值处理
噪声会把简单的阈值处理问题变为不能解决的问题。当不能再源头降低噪声且阈值处理时所选择的分割方法时，增强性能的一种常用技术时在阈值处理前先对图像进行平滑。
在分割且平滑后的图像中，物体和背景间的边界稍微有点失真，这是由平滑过图像的边界模糊导致的。事实上，对图像平滑得越强烈，分割结果中预期的边界误差就越大。

使用边缘改进全局阈值处理
1.使用前述的任意一种检测器由图像f计算一幅边缘图像。边缘图像可以是梯度或拉普拉斯算子的绝对值。
2.指定一个阈值T。
3.使用T对边缘图像进行阈值处理，产生一幅二值图像g_T。这幅图像在步骤4中用做一幅标记图像，以便从f中选取对应于“强”边缘像素的像素。
4.仅使用f中对应于g_T中1值像素位置的像素来计算直方图。
5.使用步骤4的直方图，例如采用Otsu方法来全局地分割f。
通常指定对应于一个百分位的T值，一般为一个高值（如第90个百分位），以便边缘图像中只有少量像素用于阈值的计算。自定义函数percentile2i可用于这一目的。该函数计算一个灰度值I，它对应于一个指定的百分位P。

I = percentile2i(h, P)

其中h是图像的直方图，P是区间[0,1]内的百分位值。输出I是对应于第P个百分位的灰度级。

某些图像中，背景占据了绝大部分位置，即0像素数量非常大，不利于观察直方图或做后续处理，此时可以将直方图中第一个元素排除在外，即有必要去除0对直方图的贡献。

h = imhist(f);
h(1) = 0;

基于局部统计的可变阈值处理
当照明背景非常不均匀时，全局阈值处理通常会失败。该问题的一种解决方法是试图估计明暗函数，并用它来补偿不均匀的灰度模式，然后使用上述方法之一对图像做全局阈值处理。在不规则光照情形下，或者子啊由多个主要物体灰度的情况下（此时全局阈值处理也有困难），进行补偿的另一种方法是采用可变阈值处理。这种方法基于(x,y)邻域中像素的一种或多种指定特性，在每个点(x,y)处计算一个阈值。
局部标准差计算函stdfilt

g = stdfilt(f,nhood)

其中f是输入图像，nhood是由0和1组成的数组，其中非零元素指定用于计算局部标准差所用的邻域。nhood的尺寸在每一个维度上必须是奇数，默认值是ones(3)。
局部均值计算函数localmean

function mean = localmean(f,nhood)
%计算局部均值
%Output MEAN is an array the same size as F containing
%the local mean ateach point
if nargin == 1nhood = ones(3) / 9;
elsenhood = nhood / sum(nhood(:));
end
mean = imfilter(tofloat(f), nhood, 'replicate');

基于局部均值m和局部标准差σ的可变局部阈值的通用公式是
T x y = a σ x y + b m x y T~xy~ = aσ~xy~ + bm~xy~ T xy =aσ xy +bm xy
式中，a和b都是非负常数，xy是下标。另一种有用的形式是
T x y = a σ x y + b m G T~xy~ = aσ~xy~ + bm~G~ T xy =aσ xy +bm G
m_G是全局图像均值，G是下标。
下面的函数使用以上公式执行局部阈值处理

function g = localthresh(f, nhood, a, b, meantype)
%根据局部均值和局部或全局标准差计算阈值
%The segmented image is given by
%       1 if (F > A*SIG) AND (F > B*MEAN)
% G =
%       0 otherwise
% If MEANTYPE = 'local'(the default),then MEAN is an array of
%local means. if MEANTYPE = 'global', then MEAN is the global
%(image) mean, a scalar.
f = tofloat(f);SIG = stdfilt(f, nhood);if nargin == 5 && strcmp(meantype, 'global')MEAN = mean2(f);
elseMEAN = localmean(f, nhood);
endg = (f > a*SIG) & (f > b*MEAN);

当背景接近于常数且所有物体的灰度高于或低于背景灰度时，选择全局均值一般会给出较好结果。

基于区域的分割

分割的约束条件

区域生长

区域生长时根据预先定义的生长准则将像素或子区域组合为更大区域的过程。基本方法是从一组“种子”点开始，将与种子性质相似的那些邻域像素附加到每个种子上来形成这些生长区域（如特定范围的灰度或颜色）。
当先验知识不可用时，这一过程是在每个像素处计算一组相同的特性，然后再生长处理过程期间，利用这组特性把像素分配到各个区域。如果计算的结果显示了一簇值，则那些性质靠近这些簇的中心的像素，可以作为种子使用。
函数regiongrow可以实现基本的区域生长

[g, NR, SI, TI] = regoingrow(f, S, T)

function [g,NR,SI,TI] = regiongrow(f,S,T)
%通过区域生长的方式分割图像
f = tofloat(f);
%If S ia a scalar, obtain the same image
if numel(S) == 1SI = f == S;S1 = S;
else%S is an array. Eliminate duplicate, connected seed locations%to reduce the number of loop executions in the following%section of code.SI = bwmorph(S, 'shrink', Inf);S1 = f(SI); % Array of seed values.
endTI = false(size(f));
for K = 1:length(S1)seedvalue = S1(K);S = abs(f - seedvalue) <= T; % Re-use variable S.TI = TI | S;
end
%Use function imreconstruct with SI as the maker image to
%obtain the regions corresponding to each seed in S. Function
%bwlabel assigns a different integer to each connected region.
[g,NR] = bwlabel(imreconstruct(SI, TI));

区域分离与聚合

刚刚讨论的过程是由一组种子点来生长区域。另一种方法是首先将图像细分为一组任意的不相交区域，然后聚合/或分离这些区域，以便满足约束。
以上内容可以总结为如下过程：
1.将满足条件P(R_i)=FALSE的任何区域R_i分离为4个不相交的象限区域。
2.当不可能进一步分离时，聚合满足条件P(R_i∪R_j)=TRUE的任意两个相邻区域R_i和R_j 。
3.当无法进一步聚合时，停止操作。

工具箱中提供的四叉树分解函数是qtdecomp

Z = qtedcomp(f, @split_test, parameters)

其中f是输入图像，Z是包含四叉树结构的一个稀疏矩阵。若Z(k,m)是非零的，贼(k,m)是分解的一个块的左上角，且该块的大小是Z(k,m)。函数split_test用于确定一个区域是否进行分离，parameters是split_test所要求的其他参数。
要得到四叉树分解中的实际四分区域像素值，我们使用函数qtgetblk

[vals, r, c] = qtgetblk(f, z, m)

其中vals是一个数组，它包含f的四叉树分解中的尺寸为m*m的块的值，Z是由qtdecomp返回的稀疏矩阵。参数r和c是包含块的左上角的行坐标和列坐标的向量。

使用分水岭变换的分割

在地理学中，分水岭是一个山脊，该山脊通过不同的水系来区分排水区域。集水盆地是把水排入河流或水库的地理区域。分水岭变换把这些概念应用到灰度图像处理中，从而解决许多图像分割问题。
理解分水岭变换要求我们把灰度图像视为一个拓扑表面，表面中f(x,y)的值被解释为高度。如果雨水降落到该表面上，那么雨水明显会流入集水盆地中。正好落在分水岭脊线上的雨水会等概率地流到集水盆地中。分水岭变换将找到灰度图像中的集水盆地和脊线。

分水岭分割方法的主要应用是从背景中提取几乎一致（类似水滴）的对象。那些在灰度级上变化较小的区域的梯度值也较小。因此，实际上，我们经常可以见到分水岭分割方法于图像的梯度有更大的关系，而不是图像本身。有了这样的方法，汇水盆地的局部最小值就可以与对应与所关注的对象的小的梯度值练习起来了。

水坝构造

水坝的构造是以二值图像为基础的，这种图像属于二维整数空间Z² 。构造水坝分离二元点集的最简单方法是使用形态学膨胀。
若两个汇水盆地在第n步成功聚合成一个连通区域，设此时的连通区域为q，则通过膨胀构建水坝的两个条件为：
1.膨胀受到q的约束，即在膨胀的过程中结构元素的中心只能定位在q中；
2.在引起聚合的那些点上不执行膨胀（成为单一的连通分量）。
在满足以上两个条件的情况下，执行膨胀的结果会在两个汇水盆地之间构建一个单像素宽的水坝。
通过这一方法得到的水坝是连通分量，就是我们希望得到的分割边界。就是说，这种方法消除了分割线间断的问题。

使用距离变换的分水岭分割

二值图像的距离变换的含义是：从每个像素到最接近非零值像素的距离。
，如图。
该方法的完整过程为

f = imread('/houses.png');
g = im2bw(f,graythresh(f));
g1  = imdilate(g, strel('square',10));%膨胀
gc = ~g1;%图像求补
D = bwdist(gc);%计算距离变换
L = watershed(-D);%计算负距离变换的分水岭变换
w = L == 0;%取出脊线
g2 = g & ~w;%对原始图像和w的补求逻辑与，即可完成分割

经过形态学处理的二值图像
距离变换图像
取出L中的脊线
合并原始二值图像和脊线图像

距离变换使用函数bwdist来计算

D = bwdist(f)

计算负距离变换的分水岭变换的函数是watershed

L = watershed(A, conn)

其中，L是标记矩阵，A是一个输入数组，conn指定连接性，默认为8连接。L中的正整数对应于集水盆地，零值指出分水岭脊线像素。

如上图，该方法存在过分割现象，是基于分水岭的分割方法的常见问题，以下两个方法可以克服这种现象。

使用梯度的分水岭分割

使用分水岭变换进行分割前，常用梯度幅度来对灰度图像进行预处理。梯度幅度图像沿物体的边缘有较高的像素值，而在其他位置则有较低的像素值。理想情况下，分水岭变换可得到沿物体边缘的分水岭脊线。
由开篇的知识，梯度是两个偏导数平方和的开平方，我们可以用线性滤波方法得到梯度图像

h = fspecial('sobel');
fd = tofloat(f);
g = sqrt(imfilter(fd,h,'replicate') .^ 2 + imfilter(fd,h,'replicate') .^ 2);

有时需要对梯度图像进行平滑操作，这里使用形态学的开-闭滤波方法

g2 = imclose(imopen(g, ones(3,3)), ones(3,3));

标记控制符的分水岭分割