一、点、线、边缘检测
1、背景知识
数字函数的导数可用差分来定义。
我们按如下方式得到一维函数 f(x) 在点 x 处的导数的近似：
对于用于 一阶 导数的任何近似，要求：
（1）在恒定灰度区域必须为零；
（2）在灰度台阶或斜坡开始处必须不为零；
（3）在沿灰度斜坡点处也必须不为零。对于用于 二阶 导数的任何近似，要求：
（1）在恒定灰度区域必须为零；
（2）在灰度台阶或斜坡开始处和结束处必须不为零；
（3）在沿灰度斜坡点处也必须为零。
相关结论
- 一阶导数通常在图像中产生较粗的边缘；
- 二阶导数对精细细节，如细线、孤立点和噪声有较强的相应；
- 二阶导数在灰度斜坡和灰度台阶过度处产生双边缘相应；
- 二阶导数的符号可用于确定边缘的过渡是从亮到暗还是从暗到亮。
2、点检测
拉普拉斯算子：
二值图像：
方法1：拉普拉斯模板过滤，即二阶导数
w = [-1 -1 -1; -1 8 -1; -1 -1 -1]
g = abs(imfilter(tofloat(f), w)) >= T;
示例
w = [-1 -1 -1; -1 8 -1; -1 -1 -1];
g = abs(imfilter(tofloat(f), w));
T = max(g(:));
g = g >= T;
imshow(g)
方法2：在所有大小为 m x n 的领域中，找到其最大像素和最小像素值的差大于指定 T 值的点。
g = ordfilt2(f, m*n, ones(m,n)) - ordfit2(f, 1, ones(m,n));
g = g >= T;
2、线检测
类似于，点检测过程，w 模板为
% 水平
w = [-1 -1 -1, 2 2 2, -1 -1 -1];
% +45
w = [2 -1 -1, -1 2 -1, -1 -1 2];
% 垂直
w = [-1 2 -1, -1 2 -1, -1 2 -1];
% -45
w = [-1 -1 2, -1 2 -1, 2 -1 -1];
3、使用函数 edge 检测边缘
<基本思想>
1、寻找灰度的一阶导数的幅度大于某个指定阈值的位置。
2、寻找灰度的二阶导数有零交叉的位置。
<edge函数>
提供了几个基于上述规则的边缘估计器，通用语法为
[g, t] = edge(f, 'method', parameters)
% method: 'sobel', 'prewitt', 'roberts', 'log', 'zerocross', 'canny'
<Sobel边缘检测器>
使用函数 imfilter 和下面左侧模板对图像 f 滤波，然后再使用另一个模板对 f 滤波，对每幅滤波后的图像的像素值取平方，将两个结果相加，并计算它们的平方根。
-1 -2 -1 -1 0 1
0 0 0 -2 0 2
1 2 1 -1 0 1
sobel 检测器的通用语法如下
[g, t] = edge(f, 'sobel', T, dir)
% dir: 'horizontal'、‘vertical’、'both'
% T: 一个指定的阈值
<Prewitt边缘检测器>
Prewitt 检测器与 Sobel 检测器相比，计算上要简单一些，但产生的结果中噪声可能会稍微大一些。
-1 -1 -1 -1 0 1
0 0 0 -1 0 1
1 1 1 -1 0 1
[g, t] = edge(f, 'prewitt', T, dir)
<Roberts边缘检测器>
与其他边缘检测器相比，Roberts 较为古老，功能有限（它是非对称的，不能检测多种 45度倍数的边缘），但在简单和速度上较为优势
-1 0 0 -1
0 1 -1 0
[g, t] = edge(f, 'roberts', T, dir)
<LoG边缘检测器>
使用高斯拉普拉斯算子与一幅图像卷积，与先用平滑函数与该图像卷积，再计算卷积结果的拉普拉斯算子，结果是相同的。
使用高斯拉普拉斯算子有两个效果：他平滑了图像（降低了噪声），并计算拉普拉斯算子，进而产生一幅双边缘图像。然后通过查找双边缘之间的零交叉来定位边缘。
[g, t] = edge(f, 'log', T, sigma)
<零交叉检测器>
这种检测器基于与 LoG 方法相同的概念，但卷积是使用一个特定的滤波器函数 H 来完成的。
[g, t] = edge(f, 'zerocross', T, H)
<Canny边缘检测器>
Canny 边缘检测器是函数 edge 中最强大的边缘检测器。方法总结如下：
1、使用具有指定标准差的一个高斯滤波器来平滑图像，目的是去除噪声。
2、在每个点处计算局部梯度和 边缘方向
，前三个检测器都可以用来计算这些导数。
3、应用非最大抑制（non-maximum suppression）技术来消除边误检（本来不是但检测出来是）
4、应用 滞后阈值处理 方法，划分出强边缘像素和弱边缘像素
5、通过将 8 连接的弱像素执行到强像素来执行边缘连接
Canny 边缘检测器的语法为
[g, t] = edge(f, 'canny', T, sigma)
T 是一个向量， T = [T1, T2]，sigma 是平滑滤波器的标准差，默认为1。
二、使用霍夫变换进行线检测
霍夫变换概念
霍夫变换最简单的是检测直线。我们知道，直线的方程表示可以由斜率和截距表示（这种表示方法，称为斜截式），如下所示：
如果用参数空间表示,即用斜率和截距就能表示一条直线。但是这样会参数问题，垂直线的斜率不存在（或无限大），这使得斜率参数mm的值接近于无限。
解决方法之一是，使用法线来表示直线：
其中 p 是原点到直线上最近点的距离，θ 是 x轴与连接原点和最近点直线之间的夹角。如图所示。
霍夫变换坐标系中，轴分别是 p 轴 和 θ 轴。
pθ 坐标系中每一条正弦曲线代表 xy 坐标系中经过某一点的所有直线。
pθ 坐标系中两条正弦曲线的交点代表 xy 坐标系中对应两点的一条直线。
霍夫变换矩阵
对于 xy 平面中的每个非背景点，求出每个 θ 轴上允许的细分值对应的 p 值。之后，将 p 值四舍五入到 p 轴上最接近的允许单元值。然后，相应的累加器单元加1。最终，累加器数组则成为 霍夫变换矩阵。
图像处理工具箱提供了 3 个与霍夫变换有关的函数。
函数 hough 实现霍夫变换的概念。
函数 houghpeaks 寻找霍夫变换中的峰值（高计数累加器单元）。
函数 houghlines 基于前两个函数的结果，提取原始图像中的线段。
hough
[H, theta, rho] = hough(f, 'ThetaRes', val1, 'RhoRes', val2)
H 是霍夫变换矩阵
theta 和 rho 是 p 和值的向量。
val1 是值在 0 到 90 之间的一个标量，它指定了沿
val2 是一个实表量，它指定了沿 p 轴的霍夫变换容器的间隔，默认为1
houghpeaks
peaks = houghpeaks(H, NumPeaks)
peaks 是容纳这些峰值的行列坐标矩阵。
H 是霍夫变换矩阵。
NumPeaks 是指定数量的峰值。
houghlines
lines = houghlines(f, theta, rho, peaks)
lines 是一个结构数组，该数组的长度等于所找到的线段数。该结构的每个元素识别一条线，并具有以下字段：
point1，point2， 线段两端终点的行列坐标
theta，与线相关的霍夫变换容器的角度
rho，与线相关的霍夫变换容器的 p 轴的位置
输入中 theta 和 rho 是函数 hough 的输出。
peaks 是函数 houghpeaks 的输出。
示例
[H, theta, rho] = hough(f);
peaks = houghpeaks(H, 5);
lines = houghlines(f, theta, rho, peaks);
三、阈值处理
图像阈值处理在图像分割的许多应用中处于核心地位。
1、基本的全局阈值处理
适用于观察的直方图主要分为两种灰度级。
首选方法是使用一种能基于图像数据 自动选择 的阈值算法。过程如下：
1、为全局阈值选择一个初始估计值 T。
2、使用 T 分割图像。这会产生两组像素：由所有灰度值大于 T 的像素组成的 G1，由所有灰度值小于等于 T 的像素组成的 G2。
3、分别计算区域 G1 和 G2 中像素的平均灰度值 m1 和 m2。
4、计算一个新的阈值：
5、重复步骤 2 到步骤 4，直到后续迭代中 T 的差小于一个预定义的值为止。
6、使用函数 im2bw 分割图像
示例
count = 0;
T = mean2(f);
done = false;
while ~done
count = count + 1;
g = f > T;
Tnext = 0.5*(mean(f(g)) + mean(f(~g)));
done = abs(T - Tnext) < 0.5;
T = Tnext;
end
g = im2bw(f, T/255);
2、使用 Otsu 方法进行最佳全局阈值处理
类间方差最大化的思想是方差越大，越接近正确分割图像的阈值。Otsu 可以找到这个最佳阈值，如果最大值不唯一，则找到所有最佳阈值的平均值。
工具箱中 graythresh 可计算 Otsu 阈值，语法为
[T, SM] = graythresh(f)
T 是产生的阈值，它被归一化到区间 [0,1]
SM 是可分性测度
示例
[T, SM] = graythresh(f);
g = im2bw(f, T);
下面自定义函数给定图像 直方图 的 T 和 SM
function [T, SM] = otsuthresh(h)
h = h/sum(h);
h = h(:);
i = (1:numel(h))';
P1 = cumsum(h);
m = cumsum(i.*h);
mG = m(end);
sigSquared = ((mG*P1 - m).^2)./(P1.*(1 - P1) + eps);
maxSigsq = max(sigSquared);
T = mean(find(sigSquared == maxSigsq));
T = (T - 1)/(numel(h) - 1);
SM = maxSigsq / (sum(((i - mG).^2).*h) + eps);
3、使用图像平滑改进全局阈值处理
噪声会把简单的阈值处理问题变为不能解决的问题。当不能再源头降低噪声且阈值处理是所选择的分割方法时，增强性能的一种常用技术是在阈值处理前先对图像进行平滑。
示例
w = fspecial('average', 5);
fa = imfilter(fn, w, 'replicate');
Ta = graythresh(fa);
ga = im2bw(fa, Ta);
4、使用边缘改进全局阈值处理
当物体比背景小得多时，他们对直方图的贡献可以或略不计，使用边缘信息可以改进这种情况。算法如下：
1、使用一种方法计算一幅边缘图像。边缘图像可以是梯度或拉普拉斯算子绝对值。
2、指定一个阈值 T，对边缘图像进行阈值处理，产生一幅二值图像 g(x, y)，以便从 f(x, y) 中选取对应于 强 边缘像素的像素。
3、仅使用 f(x,y) 中对应于 g(x,y) 中 1 值像素的位置的像素来计算直方图。
4、对直方图采用 Otsu 方法来全局地分割 f(x,y)。
通常指定对应于一个百分位的 T 值，它通常设置一个高值（如第90个百分位），以便边缘图像中只有很少的像素用于阈值的计算。
自定义函数 percentile2i 可实现这一目的。该函数计算一个灰度值 I，它对应于一个指定的百分位 P。语法为
I = percentile2i(h, P)
h 是图像的直方图
P 是区间 [0, 1] 内的百分位值
I 是对应于第 P 个百分位的灰度级
function I = percentile2i(h, P)
if P < 0 || P > 1
error('The percentile must be in the range [0, 1].')
end
h = h/sum(h);
C = cumsum(h);
idx = find(C >= P, 1, 'first');
I = (idx - 1)/(numel(h) - 1);
示例
- 使用基于梯度的边缘信息改进全局阈值处理
f = tofloat(imread('a.tif'));
sx = fspecial('sobel');
sy = sx'
gx = imfilter(f, sx, 'replicate');
gy = imfilter(f, sy, 'replicate');
grad = sqrt(gx .* gx + gy .* gy);
grad = grad/max(grad(:));
h = imhist(grad);
Q = percentile2i(h, 0.999);
markerImage = grad > Q;
fp = f .* markerImage;
hp = imshist(fp);
hp(1) = 0;
T = otsuthresh(hp);
g = im2bw(f, T);
- 使用拉普拉斯算子边缘信息来改进全局阈值处理
f = tofloat(imread('a.tif'));
w = [-1 -1 -1; -1 8 -1; -1 -1 -1];
lap = abs(imfilter(f, w, 'replicate'));
lap = lap/max(lap(:))
h = imhist(lap);
Q = percentile2i(h, 0.995);
markerImage = lap > Q;
fp = f .* markerImage;
hp = imshist(fp);
hp(1) = 0;
T = otsuthresh(hp);
g = im2bw(f, T);
5、基于局部统计的可变阈值处理
在不规则光照情形下，或者在有多个主要物体灰度的情况下（此时全局阈值处理有困难），进行补偿的一种方法是采用可变阈值处理。
令 和
分别表示包含在一幅图像中以坐标 （x,y）为中心的领域的一组像素的标准差和均值。
- 计算局部标准差，使用函数 stdfilt：
g = stdfilt(f, nhood)
nhood 是由 0 和 1 组成的数组，其中非零元素指定用于计算局部标准差所用的领域。nhood 的尺寸在每个维度上必须是奇数，默认值是 ones(3)。
- 计算局部均值，使用如下自定义函数：
function mean = localmean(f, nhood)
if nargin == 1
nhood = ones(3) / 9;
else
nhood = nhood / sum(nhood(:));
end
mean = imfilter(tofloat(f), nhood, 'replicate');
令 ，或者将
替换为全局图像均值
.
则分割后的图像计算如下：
下面自定义函数可执行上式进行局部阈值处理
function g = localthresh(f, nhood, a, b, meantype)
f = tofloat(f);
SIG = stdfilt(f, nhood);
if nargin == 5 && strcmp(meantype, 'global')
MEAN = mean2(f);
else
MEAN = localmean(f, nhood);
end
g = (f > a*SIG) & (f > b*MEAN);
示例
g = localthresh(f, ones(3), 30, 1.5, 'global');
6、使用移动平均的图像阈值处理
局部阈值处理方法的一个特殊情况是沿一幅图像的扫描线来计算移动平均。这一实现在文本处理中相当有用。
四、基于区域的分割
1、区域生长
区域生长是根据预先定义的生长准则将像素或子区域组合为更大的区域的过程。
基本方法是从一组 “种子” 点开始，将与种子性质相似的那些邻域像素附加到每个种子上来形成这些生长区域。
下面这个自定义函数 regiongrow 用来实现基本的区域生长，语法为
[g, NR, SI, TI] = regiongrow(f, S, T)
S： 可以是一个数组（与 f 大小相同）或一个标量。若 S 是一个数组，则它在种子点的所有坐标处必须包含 1，而在其他地方包含 0。如果 S 是一个标量，则它定义一个灰度值，那么 f 中具有该值的所有点都会成为种子点。
T： 可以是一个数组或一个标量。若 T 是一个数组，则它对 f 中的每个位置包含一个阈值。若 T 是一个标量，则它定义一个全局阈值。S 和 T 所有值必须标定到区间 [0, 1]，且与输入图像的类无关。
g： 是分割后的图像。
NR： 是所找到的区域的数量。
SI： 是包含种子点的一幅图像。
TI： 是包含处理连接性前就已通过阈值测试的像素的一幅图像。
function [g, NR, SI, TI] = regiongrow(f, S, T)
f = tofloat(f);
if numel(S) == 1
SI = f == S;
S1 = S;
else
SI = bwmorph(S, 'shrink', Inf);
S1 = f(SI);
end
TI = false(size(f));
for K = 1:length(S1)
seedvalue = S1(K);
S = abs(f - seedvalue) <= T;
TI = TI | S;
end
[g, NR] = bwlabel(imreconstruct(SI, TI));
2、区域分离与聚合
另一种方法是首先将一幅图像细分为一组任意的不相交区域，然后聚合或分离这些区域。过程如下：
1、令 R 表示整幅图像区域，选择一个属性 P。将满足
的任何区域
分离为 4 个不相交的象限区域。
2、当不可能进一步分离时，聚合满足条件的任意两个相邻区域
。
3、当无法进一步聚合时，停止操作。
工具箱中实现四叉树分解的函数是 qtdecomp。
Z = qtdecomp(f, @split_test, parameters)
Z 是包括四叉树结构的一个稀疏矩阵。若 Z(k, m) 是非零的，则 Z(k, m) 是分解的一个块的左上角，且该块的大小是 Z(k, m) 。
split_test 用于确定一个区域是否进行分离。
函数 qtgetblk 可得到四叉树分解中的实际四分区域像素值。
[vals, r, c] = qtgetblk(f, Z, m)
vals 是一个数组，包含 f 的四叉树分解中的尺寸为 m*m 的块的值。
Z 是由 qtdecomp 返回的稀疏矩阵。
r, c 是包含块的左上角的行坐标和列坐标的向量。
下面是一个基本的分离、聚合函数，说明 qtdecomp 的用法。
g = splitmerge(f, mindim, @predicate)
g 是输出图像，输出图像每个连接区域都使用不同的整数标注。
mindim 定义分解中允许的最小块尺寸，必须是 2 的正整数次幂。
flag = predicate(region)
如果在 region 区域中的像素满足函数中由代码定义的属性，那么函数就必须写成返回 true 的形式，否则，flag 的值就必须是false。
function g = splitmerge(f, mindim, fun)
Q = 2^nextpow2(max(size(f)));
[M, N] = size(f);
f = padarray(f, [Q-M, Q-N], 'post');
Z = qtdecomp(f, @split_test, mindim, fun);
Lmax = full(max(Z(:)));
g = zeros(size(f));
MARKER = zeros(size(f));
for K = 1:Lmax
[vals, r, c] = qtgetblk(f, Z, K);
if ~isempty(vals)
for I = 1:length(r)
xlow = r(I); ylow = c(I);
xhigh = xlow + K - 1; yhigh = ylow + K - 1;
region = f(xlow:xhigh, ylow:yhigh);
flag = fun(region);
if flag
g(xlow:xhigh, ylow:yhigh) = 1;
MARKER(xlow, ylow) = 1;
end
end
end
end
g = bwlabel(imreconstruct(MARKER, g));
g = g(1:M, 1:N);
%---------------------------------------------%
function v = split_test(B, mindim, fun)
k = size(B, 3);
v(1:k) = false;
for I = 1:k
quadregion = B(:, :, I);
if size(quadregion, 1) <= mindim
v(I) = false;
continue
end
flag = fun(quadregion);
if flag
v(I) = true;
end
end
%---------------------------------------------%
function flag = predicate(region)
sd = std2(region);
m = mean2(region);
flag = (sd > 10) & (m > 0) & (m < 125);
五、使用分水岭变换的分割
分水岭变换将找到灰度图像中的集水盆地和脊线。集水盆地就是我们要识别的物体或区域。
使用距离变换的分水岭分割
二值图像的距离变换是每个像素到最接近非零像素的距离。
首先，把图像转换为二值图像。
g = im2bw(f, graythresh(f));
对图像求补，计算其距离变换。
gc = ~g;
D = bwdist(gc);
接着，计算负距离变换的分水岭变换。
L = watershed(-D);
L 的零值即分水岭的脊线像素。
w = L == 0;
最后，原始二值图像和图像 w 的 “补” 的逻辑 “与” 操作可完成分割。
g2 = g & ~w;
使用梯度的分水岭分割
计算图像的幅度梯度
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));
计算梯度的分水岭变换，并找到分水岭脊线。
L = watershed(g2);
wr = L == 0;
将分水岭脊线叠加到原始图像中。
f2 = f;
f2(wr2) = 255;
标记符控制的分水岭分割
由于噪声和梯度的其他局部不规则性，常会导致过分分割。一种方法是基于标记符修改梯度图像。
前面步骤与上面一样
h = fspecial('sobel');
fd = tofloat(f);
g = sqrt(imfilter(fd, h, 'replicate') .^ 2 + imfilter(fd, h, 'replicate') .^ 2);
计算图像中所有局部最小区域的位置
rm = imregionalmin(g);
获取图像内部标记符
im = imextendedmin(f, 2)
fim = f;
fim(im) = 175;
2 是高度阈值。
获取图像的外部标记符
Lim = watershed(bwdist(im));
em = Lim == 0;
通过将局部最小区域叠加在内部和外部标记的位置来修改梯度图像
g2 = imimposemin(g, im | em);
计算修改后的梯度图像的分水岭变换
L2 = watershed(g2);
f2 = f;
f2(L2 == 0) = 255;
