实验三 图像分割与边缘检测

1、实验三 图像分割与边缘检测上课老师:赵欢喜 实验指导:吴磊 一实验目的1. 理解图像分割的基本概念；2. 理解图像边缘提取的基本概念；3. 掌握进行边缘提取的基本方法；4. 掌握用阈值法进行图像分割的基本方法。二实验基本原理图象边缘检测图像理解是图像处理的一个重要分支，研究为完成某一任务需要从图像中提取哪些有用的信息，以及如何利用这些信息解释图像。边缘检测技术对于处理数字图像非常重要，因为边缘是所要提取目标和背景的分界线，提取出边缘才能将目标和背景区分开来。在图像中，边界表明一个特征区域的终结和另一个特征区域的开始，边界所分开区域的内部特征或属性是一致的，而不同的区域内部的特征或属性是不同的，

2、边缘检测正是利用物体和背景在某种图像特性上的差异来实现的，这些差异包括灰度，颜色或者纹理特征。边缘检测实际上就是检测图像特征发生变化的位置。图象边缘检测必须满足两个条件：一能有效地抑制噪声；二必须尽量精确确定边缘的位置由于噪声和模糊的存在，检测到的边界可能会变宽或在某些点处发生间断，因此，边界检测包括两个基本内容：首先抽取出反映灰度变化的边缘点，然后剔除某些边界点或填补边界间断点，并将这些边缘连接成完整的线。边缘检测的方法大多数是基于方向导数掩模求卷积的方法。导数算子具有突出灰度变化的作用，对图像运用导数算子，灰度变化较大的点处算得的值比较高，因此可将这些导数值作为相应点的边界强度，通过设置门

3、限的方法，提取边界点集。一阶导数与是最简单的导数算子，它们分别求出了灰度在x和y方向上的变化率，而方向上的灰度变化率可以用相应公式进行计算;对于数字图像，应该采用差分运算代替求导。一幅数字图像的一阶导数是基于各种二维梯度的近似值。图像f(x,y)在位置(x,y)的梯度定义为下列向量： (3-4)在边缘检测中，一般用这个向量的大小，用表示 (3-5)函数f在某点的方向导数取得最大值的方向是，方向导数的最大值是称为梯度模。利用梯度模算子来检测边缘是一种很好的方法，它不仅具有位移不变性，还具有各向同性。为了运算简便，实际中采用梯度模的近似形式。 或者 传统的边缘检测算法通过梯度算子来实现的，在求边缘

4、的梯度时，需要对每个象素位置计算。在实际中常用小区域模板卷积来近似快速计算，简单有效，即梯度算子一般采用滤波算子的形式来完成，因此应用很广泛。模板是N*N的权值方阵，经典的梯度算子模板有：Sobel模板、Prewitt模板、Roberts模板、Laplacian模板等。具体模板请见书。拉普拉斯高斯（LoG）算法是一种二阶边缘检测方法。它通过寻找图像灰度值中二阶微分中的过零点（Zero Crossing）来检测边缘点。其原理为，灰度级变形成的边缘经过微风算子形成一个单峰函数，峰值位置对应边缘点；对单峰函数进行微分，则峰值处的微分值为0，峰值两侧符号相反，而原先的极值点对应二阶微分中的过零点，通过

5、检测过零点即可将图像的边缘提取出来。(a)原图 (b)边缘检测后的图 (c) 阈值处理后的图图3-1 检测具有-45度边缘的图例1点检测原理：常数灰阶区域中的某孤立点对某种模板的响应绝对值肯定是最强烈的。最常用的模板有：Matlab 实现方法：g = abs(imfilter(double(f),w) >= T; where w is a appropriate point detection mask which satisfies the above condition.实例：the detection of isolated bright point in the dark gra

6、y area of the northeast quadrant.(image size: 675*675)实例代码：f=imread('Fig1002(a)(test_pattern_with_single_pixel).tif');w = -1 -1 -1; -1 8 -1; -1 -1 -1;g = abs(imfilter(double(f), w);T = max(g(:);g = g>= T;subplot(121);imshow(f);title('the original image');subplot(122);imshow(g(1:en

7、d-400, 400:end);title('the isolated point detected (only a part displayed)');=2线 (通常假定一个象素厚度) 检测 原理与上同，典型模板有（主要方向性）：实例：-450 方向线的检测： 3边沿检测方法：使用一阶或者二阶导数。对一节导数，关键问题是怎样估计水平和垂直方向的梯度Gx 和Gy，二阶导数通常使用Laplacian算子计算，但是Laplacian算子很少单独用来检测边缘，因为其对噪声非常敏感，而且其结果会产生双边沿，加大了边缘检测的困难。然而，如果Laplacian算子能与其他边缘检测算法相结合

8、，如边缘定位算法，则其是一个强有力的补充。通常两个标准用来测度图像强度的迅速变化：(1) 找出强度的一阶导数值大于某个事先阈值标准的位置；(2) 找出图像二阶导数的跨零点。IPT 工具箱函数edge 提供了几种基于上面两种标准的估计器：其语法为：g, t = edge(f, method, parameters);这里 method 参数包括这几种类型的边缘检测子：Sobel, Prewitt, Roberts, Laplacian of a Gaussian (LoG), Zero crossings and Canny，前三种的模板见下图： 另一个强有力的边缘检测器：Canny Edge

9、Detector (Canny 1986)，其算法的基本步骤如下：(1) First, the image is smoothed using a Gaussian filter with a specified standard deviation s(2) The local gradient, g(x, y) = Gx2+Gy21/2, and edge direction, q(x, y) = tan-1(Gy /Gx), are computed at each point. Any of the first three techniques can be used to compu

10、ter the Gx and Gy. An edge point is defined to be a point whose strength is locally maximum in the direction of the gradient.(3) The edge points give rise to ridges in the gradient magnitude image. The algorithm then tracks along the top of these ridges and sets to zero all pixels that are not actua

11、lly on the ridge top so as to give a thin line, a process known as nonmaximal suppression. The ridge pixels are the thresholded using thresholds, T1 and T2, with T1 < T2. Ridge pixels with values greater than T2 are said to be “strong” edge pixels and pixels between T1 and T2 “weak” edge pixels.(

12、4) Finally, the algorithm performs edge linking by incorporating the weak pixels that are 8-connected to strong pixels.注意：Edge function does not compute edges at ±450. To compute edges we need to specify the mask and use imfilter. 4Hough变换In practice, the resulting pixels produced by the method

13、s discussed in the previous sections seldom characterize an edge completely because of noise, breaks from nonuniform illumination, and other effects that introduce spurious discontinuities. Hough Transform is one type of linking procedure to find and link line segments for assembling edge pixels int

14、o meaningful edges.About the principle of Hough transform, please refer to page 586 in textbook.Instance of Hough transform:% constructing an image containing 5 isolated foreground pixels in several locaitons:f = zeros(101, 101);f(1, 1) = 1, f(101, 1) = 1, f(1, 101) = 1, f(101, 101) = 1, f(51, 51) =

15、 1;H, theta, rho = hough(f); % hough transformimshow(theta, rho, H, , 'notruesize');axis on, axis normal;xlabel('theta'), ylabel('rho');图象分割图像分割是将图像划分成若干个互不相交的小区域的过程， 小区域是某种意义下具有共同属性的像素的连通集合。如不同目标物体所占的图像区域、 前景所占的图像区域等。连通是指集合中任意两个点之间都存在着完全属于该集合的连通路径。单色（灰度）图像的分割通常是基于图像强度的两个基本特征

16、：灰阶值的不连续性和灰度区域的相似性。第一类方法主要是基于图像灰阶值的突然变换（如边缘）来分割图像，而第二类方法主要是把图像的某个子区域与某预定义的标准进行比较，以二者之间的相似性指标为指导来划分图像区域：如阈值化技术、面向区域的方法、形态学分水岭分割算法等。1.双峰法先给出原图的直方图，再定出阈值(门限)T，一般取两个峰值间的谷值。2.P参数法这种方法用于目标所占图像面积已知的情况。设目标在最简单图像f(i , j) 中所占的面积s0与图像面积s之比为P = s0/ s，则背景所占面积比为1-P = (s - s0) / s。一般来说，低灰度值为背景，高灰度值为目标。如果统计图像f(i ,

17、j)灰度值不大于某一灰度t的像元数和图像总像元数之比为1-p时，则以t为阈值。3自适应全局阈值（单阈值）算法步骤如下：1、 初始化阈值T (一般为原图像所有像素平均值)。2、 用T分割图像成两个集合：G1 和G2，其中G1包含所有灰度值小于T的像素，G2包含所有灰度值大于T的像素。3、 计算G1中像素的平均值m1及G2中像素的平均值m2。4、 计算新的阈值：T （m1m2）/2 。5、 如果新阈值跟原阈值之间的差值小于一个预先设定的范围，停止循环，否则继续24步。全局单阈值分割只适用于很少的图像。对一般图像采用局部阈值法或多阈值法会得到更好的效果4. 最大类间方差法 (OTSU)设有M-1个阈

18、值：0k1k2KM-1L-1。将图像分割成M个灰度值的类Cj， （Cjkj-1+1, , kj；j=1, 2, , M ; k0=0, kM=L），则各类Cj的发生概率j和平均值j为 (3-1) (3-2)式中, (0)=0，(0)=0。 由此可得各类的类间方差为 (3-3)将使上式的2值为最大的阈值组(k1, k2， , kM1)， 作为M值化的最佳阈值组。若取M为2，即分割成2类，则可用上述方法求出二值化的阈值。三实验提示1.MATLAB的图像处理工具箱中提供的edge函数可以实现检测边缘的功能，其语法格式如下：BW = edge (I,sobel)BW = edge (I,sobel,

19、thresh,direction)BW = edge (I, Roberts)BW = edge (I, log)BW = edge (I, log, thresh, sigma)BW = edge (I, canny)BW = edge (I,canny, thresh, sigma)这里BW = edge(I,'sobel')采用Sobel算子进行边缘检测。BW = edge(I,'sobel',direction)可以指定算子方向，即：direction= horizontal ，为水平方向；direction= vertical ， 为垂直方向；dire

20、ction= both ， 为水平和垂直两个方向。BW = edge(I,'canny',thresh,sigma)和BW = edge(I,'log',thresh,sigma)分别为用canny算子和拉普拉斯高斯算子进行边缘检测。选项log和canny例：用三种算子进行边缘检测。I=imread('eight.tif');imshow(I)BW1=edge(I,'roberts');figure ,imshow(BW1),title('用Roberts算子')BW2=edge(I,'sobel'

21、);figure,imshow(BW2),title('用Sobel算子 ')BW3=edge(I,'log');figure,imshow(BW3),title('用拉普拉斯高斯算子') 2也可以通过滤波方式实现图像的图象边缘，如拉普拉斯模板（有好几个）：i=imread('w01.tif')I=doubole(i);h=0,1,0;1,-4,0;0,1,0j= conv2(I,h,same);k=I-j;imshow(k, );3I>T 结果为一个逻辑值矩阵，I中大于T的值对应的位置为1（真），其余位置为0（假）。I(I

22、> T) 表示I 中所有大于T的值组成的向量。255* (I>T) + 0*(I<=T) 可将I中大于T的像素值设为255，小于等于T的像素值设为0。四、练习1 点、线和边缘检测111 点检测 点检测模板w：-1-1-1-18-1-1-1-1检测方法： g=abs(imfilter(double(f), w)>=T练习1f=imread(moon.tif);w=-1 -1 -1; -1 8 -1; -1 -1 -1;g=abs(imfilter(double(f), w);T=max(g(:);T=T*0.5;g=g>=T;imshow(f); figure, i

23、mshow(g);线检测水平模板、+45度模板、垂直模板、-45度模板。练习2f=imread(circbw.tif);imshow(f);w=2 -1 -1; -1 2 -1; -1 -1 2;g=abs(imfilter(double(f), w);figure,imshow(g);使用edge函数的边缘检测语法：g,t=edge(f, method, parameter)说明：g是一个逻辑数组，其值为：在f中检测到边缘的位置为1，其他位置为零；t是edge是用的阈值；method为边缘监测器方法，可选为： sobel, prewit, roberts, log(LoG), zerocos

24、s, canny等；parameter包含两部分：T为指定的阈值，第二部分为dir（检测边缘的首选方向： horizontal, vertical, both），或sigma（标准方差），或H（指定的滤波函数）。练习3f=imread(rice.tif);imshow(f);gsobel,t=edge(f, sobel);figure, imshow(gsobel);glog,t=edge(f, log);figure, imshow(glog);gcanny,t=edge(f, canny);figure, imshow(gcanny);22 使用Hough变换的线检测练习4设计与实现一个基

25、于Hough变换的直线检测器。23 阈值处理231 全局阈值处理语法：T=graythresh(f)说明：T是阈值，归一化为0至1之间的值。232 局部阈值处理通过一个形态学顶帽算子并对得到的结果使用graythresh来计算。练习5f=imread(moon.tif);imshow(f);T=graythresh(f);g=f>=T;figure, imshow(g);24 基于区域的分割241 区域生长 242 区域分裂和合并练习6设计与实现一个基于区域生长的分割程序。25 使用分水岭变换的分割练习7f=imread(cell.tif);imshow(f);g=im2bw(f, gr

26、aythresh(f);figure,imshow(g);gc=g;D=bwdist(gc);L=watershed(-D);w=L= =0;g2=g&w;figure,imshow(g2);26 分割后处理语法：BW2 = bwfill(BW1,c,r,n)说明：填充二进制图像的背景色。（形态学处理）练习8BW1 =1 0 0 0 0 0 0 01 1 1 1 1 0 0 01 0 0 0 1 0 1 01 0 0 0 1 1 1 01 1 1 1 0 1 1 11 0 0 1 1 0 1 01 0 0 0 1 0 1 01 0 0 0 1 1 1 0BW2 = bwfill(BW1

27、,3,3,8)I = imread('blood1.tif');BW3 = im2bw(I);BW4 = bwfill(BW3,'holes');imshow(BW3)figure, imshow(BW4)语法：bwareaopen说明：二进制图像区域打开，清除小物体。五. 部分参考程序和参考结果1房屋轮廓描绘代码：f = imread('Fig1006(a)(building).tif');gv, t = edge(f, 'sobel', 'vertical'); % using threshold comput

28、ed automatically, here t = 0.0516subplot(231);imshow(f, );title('the original image');subplot(232);imshow(gv, );title('vertical edge with threshold determined automatically');gv1 = edge(f, 'sobel', 0.15, 'vertical'); % using a specified threshold.subplot(233);imshow(g

29、v1, );title('vertical edge with a specified threshold');gboth = edge(f, 'sobel', 0.15); % edge detection of two directionssubplot(234);imshow(gboth, );title('horizontal and vertical edge');% edge detection of 450 direction using imfilter functionw45 = -2 -1 0; -1 0 1; 0 1 2;g

30、45 =imfilter(double(f), w45, 'replicate');T = 0.3*max(abs(g45(:);g45 = g45 >= T;subplot(235);imshow(g45, );title('edge at 45 with imfilter');wm45 = 0 1 2; -1 0 1; -2 -1 0;g45 =imfilter(double(f), wm45, 'replicate');T = 0.3*max(abs(g45(:);g45 = g45 >= T;subplot(236);imsh

31、ow(g45, );title('edge at -45 with imfilter');另一个实验：为比较三种检测方法的相对性能：Sobel, LoG 和 Canny edge detectors，和为了改善检测效果所需使用的技巧。% using the default thresholdf = imread('Fig1006(a)(building).tif');gs_default, ts = edge(f, 'sobel'); % ts = 0.074gl_default, tl = edge(f, 'log'); % t

32、l = 0.002 and the default sigma = 0.2gc_default, tc = edge(f, 'canny'); % tc = 0.0189 0.047 and the default sigma = 0.1% using the optimal threshold acquired by manual testgs_best = edge(f, 'sobel', 0.05);gl_best = edge(f, 'log', 0.003, 2.25);gc_best = edge(f, 'canny'

33、, 0.04 0.1, 1.5); The left column in above figure shows the edge images obtained using the default syntax for the sobel, log and canny operator respectively, whereas the right column are the results using optimal threshold and sigma values obtained by try.2Hough变换用于线检测从而增强边缘的连续性2.1 Hough transform f

34、or peak detectionPeak detection is the first step in using Hough transform for line detection and linking. However, finding a meaningful set of distinct peaks in a Hough transform can be challenging. Because of the quantization in space of the digital image and in parameter space of the Hough transf

35、orm, as well as the fact that edges in typical images are not perfectly straight, Hough transform peaks tend to lie in more than one Hough transform cell. One strategy to overcome this problem is following:(1) find the HT cell containing the highest value and record its location;(2) suppress (set to

36、 zero) HT cells in the immediate neighborhood of the maximum;(3) repeat until the desired number of peaks has been found, or until a specified threshold has been reached.function r, c, hnew = houghpeaks(h, numpeaks, threshold, nhood)%HOUGHPEAKS Detect peaks in Hough transform.% R, C, HNEW = HOUGHPEA

37、KS(H, NUMPEAKS, THRESHOLD, NHOOD) detects% peaks in the Hough transform matrix H. NUMPEAKS specifies the% maximum number of peak locations to look for. Values of H below% THRESHOLD will not be considered to be peaks. NHOOD is a% two-element vector specifying the size of the suppression% neighborhood

38、. This is the neighborhood around each peak that is% set to zero after the peak is identified. The elements of NHOOD% must be positive, odd integers. R and C are the row and column% coordinates of the identified peaks. HNEW is the Hough transform% with peak neighborhood suppressed. % If NHOOD is omi

39、tted, it defaults to the smallest odd values >=% size(H)/50. If THRESHOLD is omitted, it defaults to% 0.5*max(H(:). If NUMPEAKS is omitted, it defaults to 1. =2.2 HT for line detection and linkingFor each peak, the first step is to find the location of all nonzero pixels in the image that contrib

40、uted to that peak. This purpose can be implemented by the following function:function r, c = houghpixels(f, theta, rho, rbin, cbin)%HOUGHPIXELS Compute image pixels belonging to Hough transform bin.% R, C = HOUGHPIXELS(F, THETA, RHO, RBIN, CBIN) computes the% row-column indices (R, C) for nonzero pi

41、xels in image F that map% to a particular Hough transform bin, (RBIN, CBIN). RBIN and CBIN% are scalars indicating the row-column bin location in the Hough% transform matrix returned by function HOUGH. THETA and RHO are% the second and third output arguments from the HOUGH function. x, y, val = find

42、(f);x = x - 1; y = y - 1; theta_c = theta(cbin) * pi / 180;rho_xy = x*cos(theta_c) + y*sin(theta_c);nrho = length(rho);slope = (nrho - 1)/(rho(end) - rho(1);rho_bin_index = round(slope*(rho_xy - rho(1) + 1);idx = find(rho_bin_index = rbin);r = x(idx) + 1; c = y(idx) + 1;The pixels associated with th

43、e locations found using houghpixles must be grouped into line segments, which is programmed into the following function:function lines = houghlines(f,theta,rho,rr,cc,fillgap,minlength)%HOUGHLINES Extract line segments based on the Hough transform.% LINES = HOUGHLINES(F, THETA, RHO, RR, CC, FILLGAP,

44、MINLENGTH)% extracts line segments in the image F associated with particular% bins in a Hough transform. THETA and RHO are vectors returned by% function HOUGH. Vectors RR and CC specify the rows and columns% of the Hough transform bins to use in searching for line% segments. If HOUGHLINES finds two

45、line segments associated with% the same Hough transform bin that are separated by less than% FILLGAP pixels, HOUGHLINES merges them into a single line% segment. FILLGAP defaults to 20 if omitted. Merged line% segments less than MINLENGTH pixels long are discarded.% MINLENGTH defaults to 40 if omitted. % LINES is a structure array whose length equals the number of% merged line segments found. Each elem