九、图像分割

一、点，线和边缘检测

1，点检测

嵌在图像的恒定或近似恒定区域中的孤立点的检测，原理上非常简单。使用图中的模板时，若|R|>=T（T是一个非负阈值），则此时在模板中心检测处了一个孤立点。

这种点检测的方法可以使用工具像函数imfilter以及图中的模板来实现。重要的要求是当孤立点位于模板中心时，模板的响应最强，而在恒定灰度区域中响应为0。
若T已给i出，则是实现的命令为

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

其中，f是输入图像，w是一个合适的点检测模板，g是包含点检测点的图像。输出图像g是logical类，它的值是0和1。

2，线检测

更复杂一些的是线检测。

上图中第一个模板在一幅图像上移动，那么它对一个像素宽的水平线的响应最为强烈。对于恒定背景，当线通过模板的中间一行时，会产生最大的响应。第二，三，四个模板依次是对45度线响应，垂直响应，负45度线响应最好。

令R1,R2,R3,R4分别表示模板的响应，如果在图像的某个点处|RI|>|Rj|,则该点与模板i所支持的方向中的一条线可能更相关。

检测指定方向的线

f = imread('Fig1004(a).tif');     
w = [2 -1 -1;-1 2 -1;-1 -1 2];          % +45°方向检测线
g = imfilter(double(f),w);
gtop = g(1:120,1:120);                  % 左上角区域
gtop = pixeldup(gtop,4);                % 通过复制像素将图像扩大gtop*4倍
gbot = g(end-119:end,end-119:end);      % 右下角区域
gbot = pixeldup(gbot,4);
g1 = abs(g);                             % 检测图的绝对值
T = max(g1(:));
g2 = g1>=T;

subplot(3,2,1);imshow(f);title('(a)连线模板图像');
subplot(3,2,2);imshow(g,[]);title('(b)+45°线处理后的结果');
subplot(3,2,3);imshow(gtop,[]);title('(c)(b)中左上角的放大效果');
subplot(3,2,4);imshow(gbot,[]);title('(d)(b)中右下角的放大效果');
subplot(3,2,5);imshow(g1,[]);title('(e)(b)的绝对值');
subplot(3,2,6);imshow(g2);title('(f)满足g>=T的所有点');

实验分析：图b中比灰度背景更暗的阴影对应于负值。在45°方向上有两条主要线段，一条在左上方，另一条在右下方。图像中线段亮度高是因为只有少数像素宽，对于一个像素宽的元件，模板响应最强。图中孤立的点也是对模板具有强烈响应的那些点。

3，使用函数edge检测边缘

虽然点检测和线检测在任何图像分割的讨论中都很重要，但是目前为止，边缘检测最通用的方法是检测灰度值的不连续性。这种不连续性使用一阶和二阶导数来检测。
边缘检测的基本思想是使用如下两个准则之一来找到图像中灰度快速变化的位置：
1，寻找灰度的一阶导数的幅度大于某个指定阈值的位置。
2，寻找灰度的二阶导数由零交叉的位置。

图像处理工具函数edge提供了几个基于上述规则的边缘估计器。对于一些估计器，能否作为边缘检测器取决于对水平，垂直是否敏感，还是对两者都敏感。该函数的通用语法:

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

函数edge可用边缘检测器

Sobel边缘检测器
一阶导数数字上可以近似为差分。Sobel算子模板如下图。


在（x,y）处响应大于T（T是一个指定的阈值），那么在该位置的像素是边缘像素。
Sobel检测器的通用语法：

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

其中，f是输入图像，T是一个指定的阈值，dir指定所检测边缘的首选方向：‘horizontal’,'vertical’或‘both’（默认值）。g是一幅逻辑图像，检测到边缘的位置其值为1，其他位置为0。输出中的参数t是可以选择的，它是被edge使用的阈值。

prewitt边缘检测
prewitt边缘检测器的模板如图


与sobel算子相比，计算过程相对简单，但产生的结果中噪声可能会稍微大一些。

Roberts边缘检测器


Roberts检测器是最古老且最简单的边缘检测器之一。与其他检测器相比，功能有限，很少使用。

LOG检测器
考虑高斯函数：

这个公式表示的是平滑函数。如果此函数和图像进行卷积，则图像会变得模糊，且模糊的程度是由σ决定的。该函数的拉氏算子是；

由于求二阶导数是线性操作，所以用这个函数卷积这幅图像与先用平滑函数对图像卷积，再对结果进行拉普拉斯计算的结果是一样的。这样会得到两种效果：平滑图像，减少了噪声，计算拉普拉斯，从而产生双边缘图像，然后在双边缘之间定位由发现的零交叉组成的边缘。
LOG检测器通用的语法是：

[g,t] = edge(f,'log',T,sigma)  %T是阈值，sigma是标准差，默认为2

零交叉检测器
这种算子与LoG方法类似，不同之处在于卷积使用特殊的滤波函数H来完成的。

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

Canny边缘检测器
这种算子是edge函数中最强的边缘检测算子：
（1）使用指定了标准差σ的高斯滤波器来平滑，用来减少噪声；
（2）局部梯度和边缘方向在每个点都进行计算。边缘点被定义为梯度方向上局部强度最大的点。
（3）在（2）中决定的边缘点在梯度幅度图像上给出脊。然后算法追踪所有脊的顶部，设置所有的不在脊的顶部的像素为零。因此在输出中给出一条细线，这是非最大值抑制处理。脊像素使用称为“滞后阈值”的技术进行阈值处理，这种技术以使用两个阈值为基础，即T1和T2，且T1<T2。值大于T2的脊像素称为强边缘像素，T1和T2之间的脊像素称为弱边缘像素。
（4）算法通过将8连接的弱像素点执行到强像素点来执行边缘连接。
调用语法;

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

使用Sobel边缘检测器

f = imread('Fig1006(a).tif');     
subplot(3,2,1),imshow(f),title('(a)原始图像');
[gv, t] = edge(f,'sobel','vertical');
subplot(3,2,2),imshow(gv),title('(b)Sobel模板处理后结果');
gv = edge(f, 'sobel', 0.15, 'vertical');
subplot(3,2,3),imshow(gv),title('(c)使用指定阈值的结果');
gboth = edge(f, 'sobel', 0.15);
subplot(3,2,4),imshow(gboth),title('(d)指点阈值确定垂直边缘和水平边缘的结果');
w45 = [-2 -1 0;-1 0 1; 0 1 2];
g45 = imfilter(double(f), w45, 'replicate');
T = 0.3 * max(abs(g45(:)));
g45 = g45 >= T;
subplot(3,2,5),imshow(g45),title('(e)-45°方向边缘');
f45= [0 1 2;-1 0 1;-2 -1 0];
h45= imfilter(double(f), f45,'replicate');
T = 0.3 * max(abs(h45(:)));
h45 = h45 >= T;
subplot(3,2,6),imshow(h45),title('(f)+45°方向边缘'); 

Sobel,LOG,Canny边缘检测器的比较

f = imread('Fig1006(a).tif');     
[g_sobel_default,ts] = edge(f,'sobel');
[g_log_default,tlog] = edge(f,'log');
[g_canny_default,tc] = edge(f,'canny');

g_sobel_best = edge(f,'sobel',0.05);
g_log_best = edge(f,'log',0.003,2.25);
g_canny_best = edge(f,'canny',[0.04 0.10],1.5);

subplot(3,2,1);imshow(g_sobel_default);title('(a)默认sobel');
subplot(3,2,2);imshow(g_sobel_best);title('(b)最佳sobel');
subplot(3,2,3);imshow(g_log_default);title('(c)默认LoG');
subplot(3,2,4);imshow(g_log_best);title('(d)最佳LoG');
subplot(3,2,5);imshow(g_canny_default);title('(e)默认canny');
subplot(3,2,6);imshow(g_canny_best);title('(f)最佳canny');

实验总结：
同样的一幅图像，使用默认的sobel边缘检测得到的图像中线条是比较其他两种是比较少的，而且获取的边缘细节不是很完整，缺失了很多重要的细节 。但是在最佳的sobel边缘检测中，所得到的图像是比较完整的，线条细节基本上都被获取并显示出来，相比较其他两种方法，所得到的边缘细节是最完整的。在canny边缘检测中，默认参数的canny检测是最好的，所获得细节是最多的，其边缘基本都可以完整的检测到。

二、使用霍夫变换进行线检测

理想情况下，通过以上方法能产生位于边缘上的像素。实际上，得到的像素通常不能完全地表征边缘，因为存在噪声，不均匀照明引起的边缘断裂。因此，在执行边缘检测算法后，通常会使用连接过程来把边缘像素组装成有意义的边缘，其中之一是霍夫变换。
图像处理工具箱提供了三个与霍夫变换有关的函数。
函数hough实现霍夫变换的概念；
函数houghpeaks寻找霍夫变换中的峰值；
函数houghlines基于前两个函数的结果，提取原始图像中的线段。
函数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,2),size(I,2))的实际量，指定沿ρ轴的霍夫变换的间距（默认为1）。
霍夫变换的说明

f=zeros(101,101);
>> f(1,1)=1;f(101,1)=1;f(1,101)=1;f(101,101)=1;f(51,51)=1;
>> H=hough(f);
>> imshow(H,[])

f 是一个简单的合成图像如下

测试结果：

[H, theta, rho] = hough(f);  
imshow(H, [], 'XData', theta, 'YData', rho ,'InitialMagnification', 'fit')  
axis on, axis normal  
xlabel('\theta'), ylabel('\rho') 

三条曲线（直线也可以考虑为曲线）在±45°处的交点指出：f中有两组三个共线的点。两条曲线在(ρ,θ)=(0,-90)、(-100,-90)、(0,0)和(100,0)处的交点指出：有4组位于垂直线和水平线上的共线点。

函数houghpeaks
语法：

peaks = houghpeaks(H,NumPeaks)

完整语法

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

“…”指出来自默认语法和 peaks 的输入是持有峰值行和列坐标的 Q×2 大小的矩 阵。Q的范围是 0 到NumPeaks，H是霍夫变换矩阵。参数val1是非负的标量，指定了 H中 的什么值被考虑为峰值；val1可以从 0 到Inf变化，默认值是0.5*max(H( : ))。参数val2 是奇整数的两元素矢量，指定量围绕峰值的邻域大小。

函数houghlines
当峰值在霍夫变换中被识别出来，则可用此函数来决定线的起始点和终点。
语法

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

完整语法：

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

theta和rho是来自函数hough的输出，peaks是函数houghpeaks的输出。输出lines是结构数组，长度等于找到的线段数。结构中的每个元素可以看成一条线，并含有下列字段：

point1：两元素向量[r1, c1]，指定了线段终点的行列坐标。
point2：两元素向量[r2, c2]，指定了线段其他终点的行列坐标。
theta：与线相关的霍夫变换的以度计量的角度。
rho：与线相关的霍夫变换的ρ轴位置。
val1：是正的标量，指定了与相同的霍夫变换相关的两条线段的距离。当两条线段之间的距离小于指定的值时，函数houghlines把线段合并为一条线段(默认的距离是20个像素)。
val2：是正的标量，指定合并的线是保留还是丢弃。如果合并的线比val2指定的值短，就丢弃(默认值是40)。

使用霍夫变换做检测和连接
在这个例子中，我们用函数hough、houghpeaks和houghlines寻找二值图像f的一组线段。首先，我们用比默认值更好的角间距(用0.2代替1.0)计算和显示霍夫变换：

[H, theta, rho] = hough(f, 'ThetaResolution', 0.2);  
imshow(H, [], 'XData', theta, 'YData', rho, 'InitialMagnification', 'fit') 
axis on, axis normal  
xlabel('\theta'), ylabel('\rho') 
peaks = houghpeaks(H, 5);  
hold on  
plot(theta(peaks(:, 2)), rho(peaks(:, 1)), ...  
         'linestyle', 'none', 'marker', 's', 'color', 'w');title('(a)带有已选择的5个峰值位置的霍夫变换');  
lines = houghlines(f, theta, rho, peaks);  
figure, imshow(f), 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]); title('(b)与霍夫变换峰值对应的线段 ');   
end