数字图像处理期末大作业剖析

数字图像处理期末大作业

一、问题描述

实现第十章中采用Hough变换来检测图像中圆的过程。，通过包括平滑（把细节去除），边缘检测（得到轮廓）以及Hough变换得到的圆，并把结果叠加到原来的灰度图像上。给出具体的过程，中间结果，最后结果，实现的代码，并写出报告。

二、图片的获取以及预处理

针对老师提供的一副硬币图片，要求检测出其中的hough圆，并叠加到原图像上以便增强图像。在检测hough圆之前，首先要对图像进行平滑处理，进行拉普拉斯变换，然后检测垂直方向，水平方向，+45度和-45度方向的边缘，将四个方向的边缘叠加起来，得到总的边缘，对该图像进行二值化，然后对得到的图像检测其hough圆，得到圆形边缘，将该图像叠加到原图像上，就实现了图像边缘增强的目的。

三、图像处理算法的基本原理以及处理结果

本实验流程图如下：

开始

读取图像

图像预处理

边缘检测

二值化

Hough变换

Hough变换结果与原图像相加

结束

1.读取图像

图像处理的第一步就是对所采集的图像进行读入，本次实验的输入图像是一幅灰度图像，不需要将图像转换成为灰度图像，直接利用函数imread （）完成。

原图像如下所示：

原图像

2.图像预处理

在图像预处理中，我们完成了两步工作，首先使用方差为1的高斯噪声对图

像进行平滑，然后进行拉普拉斯变换，即)],(*)([2

y x f r h ?，22

2r 2

e 21)(σσ

-=πr h 为方差为2σ的高斯噪声，本实验中12=σ。又),(*)]([)],(*)([2

2y x f r h y x f r h ?=?，

其中2

224

2

22]2[

)(σσσ

r e

r r h -

-=?，将)(2r h ?和),(y x f 分别进行傅里叶别换，将其逐

点相乘，再进行傅里叶反变换，就得到了预处理后的图像。 3.边缘检测

对水平，垂直，+45度，-45度方向进行边缘检测，本实验中我们采用了Prewitt 梯度算子，它用于检测水平方向，垂直方向，+45度方向和-45度方向的掩膜分别如下：

水平掩膜 垂直掩膜 +45度掩膜 -45度掩膜

使用这四个掩膜分别对上一步得到的图像逐点进行处理，就可以得到四个方向的边缘了（本实验中边缘的一个像素都不处理），再将它们加起来，就得到了总的边缘，实验结果如下：

水平边缘

垂直边缘

-45度边缘

总的边缘如下图所示：

4.二值化

对上图得到的图像进行二值化，这里我采用的是循环方式确定图像全局阈值，即首先以图像的平均值作为阈值，将图像分成两部分，分别求两部分的平均值，新的阈值为这两个平均值的均值，重复上述过程，直到两次阈值之差小于特定的值时停止，并以最后一次得到的阈值对图像进行二值化，本实验中我要求两次阈值之差小于0.5时停止，最后得到的全局阈值为 -102.1332，二值化后的图像如下所示：

二值化后的图像

5.Hough变换检测圆形边界

Hough 变换的原理就是利用图像全局特征将边缘像素连接起来组成区域封闭边界，它将图像空间转换到参数空间，在参数空间对点进行描述，达到检测图像边缘的目的。该方法把所有可能落在边缘上的点进行统计计算，根据对数据的统计结果确定属于边缘的程度。其实质就是对图像进行坐标变换，把平面坐标变换为参数坐标，使变换的结果更易识别和检测。

Hough变换的基本原理在于，利用点与线的对偶性，将图像空间的线条变为参数空间的聚集点，从而检测给定图像是否存在给定性质的曲线。圆的方程为：(x-a)^2+(y-b)^2=r^2，通过Hough变换，将图像空间对应到参数空间。把X-Y 平面上的圆转换到a-b-r参数空间，则图像空间中过(x, y)点圆对应参数空间中，高度r变化下的一个三维锥面，如下

同理，过图像空间中任意一点的圆对应于参数空间中的一个三维锥面。因此，过图像空间上同一圆上的点，对应的参数空间中的三维锥面，在r高度必然相交于一点(a, b, r)。这样通过检测这一点可以得到圆的参数，相应的圆也可求得了。图像平面的方程转化为参数平面上的示意图如图所示：

Hough变换是基于通过提取分布于目标圆周上的参数及点的特征值的来检测圆或圆弧的。为了检测目标外形,对图像上的每点定义一个参数空间的映射。Hough变换通过在参数空间找寻特征(峰值或最大值点)得到位于图像空间中的特征(目标形状)来转换问题。

对于已知半径的圆Hough变换可以检测任意已知表达形式的曲线，关键在于选择合适的参数空间。我们可以根据曲线的表达形式决定其参数空间。当检测某一已知半径的圆时，可以使用与原图像空间相同的空间作为其参数空间。则原图像空间中的一个圆对应参数空间中的一个点，参数空间的一个点对应图像空间中一个圆，原图像空间中在同一圆上的点，它们的参数相同即a，b相同，它们在参数空间对应的圆就会过同一点（a，b），因此，将原图像中的所有点变换到参数空间之后，依据参数空间中点的聚集度就可判断出原图像空间中有无近似于圆的图形。

对于未知半径的圆，在一个xy平面图像中确定一个圆至少需要三个元素，即圆心的x轴和y轴坐标，圆的半径，因此Hough变换检测圆的目的就是检测出图像中各个圆的圆心坐标以及圆的半径。其基本思想是将原图像空间中边缘点映射至参数空间中，再将参数空间中得到的全部坐标点元素所对应的累加值进行统计，并根据此累加值来判断圆的大小和圆心的位置。例如，在xy平面上的方程为(x-a)2+(y-b)2=r2 其中点（a，b）为圆心坐标，r为圆半径，点（x，y）为圆周上的一点，将其转换为参数坐标系（a，b，r），方程为(a-x)2+(b-y)2=r2,可看出次方程为圆锥面，对于原图像中任意确定的一个点在参数空间都有一个三维锥面与其对应。

在Hough变换检测圆时，可以利用梯度信息在很大程度上加快圆检测的速度。对圆周而言，其梯度方向只有背离圆心或者指向圆心，当梯度指向圆心，圆心就在梯度的延长线上，而当梯度背离圆心，圆心则在梯度的反向延长线上。所以，边缘梯度信息的加入可以预估圆心的位置，这样可以使算法的运算量明显减少，并且可以有效抑制虚假局部最大值。圆心位置可以用极坐标方程形式表达:a= x - r·co s (θ(x , y ) ), a=y - r·sin (θ(x , y ) )或者a= x+ r·co s (θ(x , y ) ), a=y+r·sin (θ(x , y ) )，边缘像素（x，y）处的梯度方向为θ(x , y )。前一组公式是梯度方向背离圆心，后一组为梯度方向指向圆心。以前一种情况为例，边缘图像中每个边缘像素点（x，y）都可以算出其相应的梯度方向θ(x , y )。通过前一个公式可算出圆心坐标（a0，b0），对于参数空间可能的半径r0，其相应的参数空间累加器单元加一，最后找到累加器的局部最大值，就得到一个圆。

本实验中，我们检测的圆半径范围为10-100，半径步长为1，角度步长为π/18,p=0.7。得到的结果如下：

hough变化后的图像

6.Hough变换结果与原图像相加

上一步得到的图像是一副逻辑图像，这一步我们首先判断上一步得到的图像每一个像素是否为1，如果是则相加后该点的像素值为255，否则该点的像素值为原图像在该点的像素值，这样就得到了hough变换后图像与原图像相加的图像，结果如下所示：

最终的结果

四、源程序

本实验采用matlab来实现，原程序如下所示：

clear all;

%1.读入图像并显示

S=imread('E:\Sweden-coins.bmp');

figure;imshow(S);title('原图像');

%2.引入均值为零，方差为a的高斯噪声对图像进行平滑，然后进行拉普拉斯变换

[m,n]=size(S);

a=1;

h=zeros(m,n);

h=double(h);

for i=1:m

for j=1:n

t=i^2+j^2;

h(i,j)=(t-2*a)*exp(-(t/(2*a)))/a^2;

end

end

H=fftshift(fft2(h));

G=fftshift(fft2(S));

W=H.*G;

w=ifft2(ifftshift(W));

figure;imshow(w);title('平滑后进行拉普拉斯变换后的图像');

%3.对水平，垂直，+45度，-45度方向进行边缘检测，采用Prewitt梯度算子

IMG=double(w);

H1=zeros(m,n);

H2=zeros(m,n);

G1=zeros(m,n);

G2=zeros(m,n);

for i=2:m-1

for j=2:n-1 H1(i,j)=-IMG(i-1,j-1)-IMG(i-1,j)-IMG(i-1,j+1)+IMG(i+1,j-1)+IMG(i+

1,j)+IMG(i+1,j+1);

H2(i,j)=-IMG(i-1,j-1)-IMG(i,j-1)-IMG(i+1,j-1)+IMG(i-1,j+1)+IMG(i,

j+1)+IMG(i+1,j+1);

G1(i,j)=IMG(i-1,j)+IMG(i-1,j+1)+IMG(i,j+1)-IMG(i,j-1)-IMG(i+1,j-1

)-IMG(i+1,j);

G2(i,j)=-IMG(i-1,j-1)-IMG(i-1,j)-IMG(i,j-1)+IMG(i+1,j)+IMG(i,j+1)

+IMG(i+1,j+1);

end

end

figure;imshow(H1);title('水平边缘');

figure;imshow(H2);title('垂直边缘');

figure;imshow(G1);title('-45度边缘');

figure;imshow(G2);title('+45度边缘');

Z=zeros(m,n);

for i=1:m

for j=1:n

Z(i,j)=H1(i,j)+H2(i,j)+G1(i,j)+G2(i,j);

end

end

figure;imshow(Z);title('提取的边缘');

%4.对提取的边缘进行二值化

T=sum(Z(:))/(m*n);

T2=20;

while T2>0.5

Y1=0;

Y2=0;

N1=0;

N2=0;

N1=double(N1);

N2=double(N2);

Y1=double(Y1);

Y2=double(Y2);

for i=1:m

for j=1:n

if(Z(i,j)>T)

Y1=Y1+double(Z(i,j));

N1=N1+1;

else

Y2=Y2+double(Z(i,j));

N2=N2+1;

end

end

end

t1=Y1/N1;

t2=Y2/N2;

T1=(t1+t2)/2;

T2=abs(T-T1);

T=T1

end

GS=zeros(m,n);

for i=1:m

for j=1:n

if Z(i,j)>T

GS(i,j)=0;

else

GS(i,j)=255;

end

end

end

figure;imshow(GS);title('二值化后的图像');

%5.进行hough变换,检测圆形边界

BW=double(GS);

r_max=100;

r_min=10;step_r=1;step_angle=pi/18;p=0.7;

size_r = round((r_max-r_min)/step_r)+1;

size_angle = round(2*pi/step_angle);

hough_space = zeros(m,n,size_r);

[rows,cols] = find(BW);

ecount = size(rows);

for i=1:ecount

for r=1:size_r

for k=1:size_angle

a = round(rows(i)-(r_min+(r-1)*step_r)*cos(k*step_angle));

b = round(cols(i)-(r_min+(r-1)*step_r)*sin(k*step_angle));

if(a>0&&a<=m&&b>0&&b<=n)

hough_space(a,b,r) = hough_space(a,b,r)+1;

end

end

end

end

max_para = max(max(max(hough_space)));

index = find(hough_space>=max_para*p);

length = size(index);

hough_circle = false(m,n);

for i=1:ecount

for k=1:length

par3 = floor(index(k)/(m*n))+1;

par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;

par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;

if((rows(i)-par1)^2+(cols(i)-par2)^2<(r_min+(par3-1)*step_r)^2+5& &...

(rows(i)-par1)^2+(cols(i)-par2)^2>(r_min+(par3-1)*step_r)^2-5)

hough_circle(rows(i),cols(i)) = true;

end

end

for k=1:length

par3 = floor(index(k)/(m*n))+1;

par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;

par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;

par3 = r_min+(par3-1)*step_r;

fprintf(1,'Center %d %d radius %d\n',par1,par2,par3);

para(:,k) = [par1,par2,par3];

end

figure;imshow(hough_circle);title('hough变化后的图像');

%6.将hough变化后提取的图像与原图像相加，增强其边缘

LST=zeros(m,n);

for i=1:m

for j=1:n

if hough_circle(i,j)==1;

LST(i,j)=255;

else

LST(i,j)=S(i,j);

end

end

end

figure;imshow(LST,[]);title('最终的结果');

五、心得体会

我们选择用MATLAB来编写，因为感觉自己对C等其他语言掌握的不是最熟练的。刚拿到题目的时候我们三个感觉还挺简单的，平滑、边缘检测以及霍夫变换都是上课时候讲过的，可是到真正着手做的时候感觉理论与实践搭不上，很多在理论想清楚，可是一旦到编程实践，有些地方就会无从下手，我们又仔细研究了一下算法，首先使用高斯函数对图像进行平滑，然后进行拉普拉斯变换，对于得到的图像，为了进行边缘检测，经过讨论之后我们选择了Prewitt梯度算子，检测水平，垂直，+45度和-45度边缘，经过反复修改参数，终于得到了比较理想的边缘，边缘粗但清晰。

在求Hough圆的时候，刚开始对算法不是太清晰，导致编程很难顺利进行下去，后来，我们又仔细看了课本，上网查了一些相关的资料，最后结果顺利地显示了。但是犹豫循环计算多，导致识别圆的时间比较长，同时耗费了大量的存储空间，这一点对实际应用存在较大的影响，需要我们不断改进。总之，这次实验使我们都收获很多，不仅熟悉了算法，对图像处理的理解又多了一分，通过编程实践，一定程度上提高了我们相关方面的知识，以后我们会多多增加实践，将理论与实践结合，不断提高我们对知识的深入理解和编程实践能力。