实验四 用Hough变换提取曲线的参数

实验目的

（1）了解边缘检测算子的原理，并利用边缘算子对图像进行检测；

（2）掌握Hough变换的基本原理。

实验内容

（1）分别将原始图像及加高斯噪声、椒盐噪声后的图像中圆形边缘检测出来；

（2）用Hough变换对边缘进行参数提取。

实验要求

（1）实验用图像文件：原始图像（houghorg.bmp）、加高斯噪声后图像（houghgau.bmp）和加椒盐噪声后图像（houghsalt.bmp）；

img

图6.3 原始图像

（2）在含有噪声的背景下，先对图像中值滤波，再进行边缘检测；

（3）将目标的边界提取出来。边缘检测算子可利用matlab自带函数实现，使用Robert、Sobel和Laplacian算子；

（4）利用Hough变换提取的参数绘制曲线，并叠加在噪声图像上。

4．实验要点：

（1）利用算子进行边缘检测：可先将加噪以后的图像进行平滑滤波，如采用9*9的掩膜模板进行中值滤波；为了对图像中图形边缘进行线性提取，可通过设置阈值将图像变为二值图像，再利用三种不同的算子（Robert、Sobel和Laplacian）来完成边缘的检测；

（2）Hough变换进行曲线参数提取：在使用三种算子对加噪后图像进行边缘检测以后，使用Hough变换对检测后图像进行参数提取，并在提取成功以后，使用提取获得的参数进行图像的重建，最后将重建图像叠加到加噪图像中。注意在进行Hough变换时，对比观察获得图像与使用算子进行边缘检测获得图像之间的区别

实验原理

Hough变换：Hough变换是1962年由Hough提出来的，用于检测图像中直线、圆、抛物线、椭圆等形状能够用一定函数关系描述的曲线，它在影像分析，模式识别等很多领域中得到了成功的应用。

Hough变换的基本原理是将影像空间中的曲线(包括直线)变换到参数空间中，通过检测参数空间中的极值点，确定出该曲线的描述参数，从而提取影像中的规则曲线。

实验结果与分析

实验结果：

（一）加高斯噪声后图像：

img

加椒盐噪声后图像：

img

（二）在含有噪声的背景下，先对图像中值滤波，再进行边缘检测

img

（三）将目标的边界提取出来。边缘检测算子可利用matlab自带函数实现，使用Robert、Sobel和Laplacian算子；

（四）利用Hough变换提取的参数绘制曲线，并叠加在噪声图像上。

Robert：

image-20220606150118140

Sobel：

image-20220606150125451

Laplacian：

image-20220606150137145

表 1 Hough变换边界提取结果

高斯噪声 椒盐噪声
圆心 半径 圆心 半径
Robert 242 176 97 241 170 98
Sobel 243 179 99 242 178 98
Laplacian 242 177 97 242 177 97

实验分析：

观察实验结果，三种算子都能有效提取边界，但使用不同的算子提取出来的边缘质量有差别。Laplacian算子在Hough变换后重建的边缘与原图像有一定偏离，这可能是由于Laplacian算子对噪声点有一定的放大作用。

由于Hough变换提取的是曲线边界，也就将之前算子检测到的图像边缘的边界滤除了，最终效果很好。程序运行用时较长，是Hough变换算法的遍历过程所致。

主要代码

主程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
clc
clear all
houghorg = imread('houghorg.bmp');
houghorg = rgb2gray(houghorg);

% 产生高斯噪声和椒盐噪声 
houghgau = imnoise(houghorg,'gaussian',0,0.01);
houghsalt = imnoise(houghorg,'salt & pepper',0.01);
imwrite(houghgau,'houghgau.bmp')
imwrite(houghsalt,'houghsalt.bmp')

% 中值滤波
z = 9;
houghgau_z = medfilt2(houghgau,[z z]); 
houghsalt_z = medfilt2(houghsalt,[z z]); 
imwrite(houghgau_z,'houghgau_z.bmp')
imwrite(houghsalt_z,'houghsalt_z.bmp')

% 二值图像，利用类间方差阈值算法对滤波处理后图像进行分割处理
t = third1(houghgau_z);
houghgau_er = houghgau_z;
houghgau_er(houghgau_er<=t) = 0;
houghgau_er(houghgau_er>=t) = 255;
imwrite(houghgau_er,'houghgau_er.bmp')

t = third1(houghsalt_z);
houghsalt_er = houghsalt_z;
houghsalt_er(houghsalt_er<=t) = 0;
houghsalt_er(houghsalt_er>=t) = 255;
imwrite(houghsalt_er,'houghsalt_er.bmp')

% % Robert完成边缘的检测 houghgau_er
% houghgau_er = double(houghgau_er);
% [height,width,channels] = size(houghgau_er);
% G = zeros(height,width);
% for i=1:height-1
%     for j=1:width-1
%        G(i,j)=sqrt((houghgau_er(i,j)-houghgau_er(i+1,j+1))^2+(houghgau_er(i+1,j)-houghgau_er(i,j+1))^2);
%     end
% end
% houghgau_er_Robert = uint8(G);
% imwrite(houghgau_er_Robert,'houghgau_er_Robert.bmp')
% 
% % Robert完成边缘的检测 houghsalt_er
% houghsalt_er = double(houghsalt_er);
% [height,width,channels] = size(houghsalt_er);
% G = zeros(height,width);
% for i=1:height-1
%     for j=1:width-1
%        G(i,j)=sqrt((houghsalt_er(i,j)-houghsalt_er(i+1,j+1))^2+(houghsalt_er(i+1,j)-houghsalt_er(i,j+1))^2);
%     end
% end
% houghsalt_er_Robert = uint8(G);
% imwrite(houghsalt_er_Robert,'houghsalt_er_Robert.bmp')


% Robert完成边缘的检测 houghgau_er
houghgau_er = double(houghgau_er);
G = edge(houghgau_er,'Roberts'); % sobel算子滤波锐化
houghgau_er_Robert = G;
imwrite(G,'houghgau_er_Robert.bmp')

% Robert完成边缘的检测 houghsalt_er
houghsalt_er = double(houghsalt_er);
G = edge(houghsalt_er,'Roberts'); % sobel算子滤波锐化
houghsalt_er_Robert = G;
imwrite(G,'houghsalt_er_Robert.bmp')

% Sobel完成边缘的检测 houghgau_er
% m = fspecial('sobel'); % 应用sobel算子图像边缘检测
% G = filter2(m,houghgau_er); % sobel算子滤波锐化
G = edge(houghgau_er,'sobel'); % sobel算子滤波锐化
houghgau_er_Sobel = G;
imwrite(G,'houghgau_er_Sobel.bmp')

% Sobel完成边缘的检测  houghsalt_er
% m = fspecial('sobel'); % 应用sobel算子图像边缘检测
% G = filter2(m, houghsalt_er); % sobel算子滤波锐化
G = edge(houghsalt_er, 'sobel'); % sobel算子滤波锐化
houghsalt_er_Sobel = G;
imwrite( G,'houghsalt_er_Sobel.bmp')


% Laplacian完成边缘的检测 houghgau_er
% m = fspecial('Laplacian'); % 应用sobel算子图像边缘检测
% G = filter2(m,houghgau_er); % sobel算子滤波锐化
G = edge(houghgau_er,'log'); % sobel算子滤波锐化
houghgau_er_Laplacian = G;
imwrite(houghgau_er_Laplacian,'houghgau_er_Laplacian.bmp')

% Laplacian完成边缘的检测 houghsalt_er
% m = fspecial('Laplacian'); % 应用sobel算子图像边缘检测
% G = filter2(m,houghsalt_er); % sobel算子滤波锐化
G = edge(houghsalt_er,'log'); % sobel算子滤波锐化
houghsalt_er_Laplacian = G;
imwrite(houghsalt_er_Laplacian,'houghsalt_er_Laplacian.bmp')

[circle_xy,circle_r] = Hough(houghgau_er_Robert)
imwrite(ht(circle_xy,circle_r), 're_houghgau_er_Robert.bmp')
re_houghgau_er_Robert = imread('re_houghgau_er_Robert.bmp');
% he_houghgau_er_Robert = re_houghgau_er_Robert;
% he_houghgau_er_Robert(houghgau==255)=255;
he_houghgau_er_Robert = houghgau;
he_houghgau_er_Robert(re_houghgau_er_Robert==255)=255;
imwrite(he_houghgau_er_Robert,'he_houghgau_er_Robert.bmp');


[circle_xy,circle_r] = Hough(houghsalt_er_Robert)
imwrite(ht(circle_xy,circle_r), 're_houghsalt_er_Robert.bmp')
re_houghsalt_er_Robert = imread('re_houghsalt_er_Robert.bmp');
% he_houghsalt_er_Robert = re_houghsalt_er_Robert;
% he_houghsalt_er_Robert(houghsalt==255)=255;
he_houghsalt_er_Robert = houghsalt;
he_houghsalt_er_Robert(re_houghsalt_er_Robert==255)=255;
imwrite(he_houghsalt_er_Robert,'he_houghsalt_er_Robert.bmp');


[circle_xy,circle_r] = Hough(houghgau_er_Sobel)
imwrite(ht(circle_xy,circle_r), 're_houghgau_er_Sobel.bmp')
re_houghgau_er_Sobel = imread('re_houghgau_er_Sobel.bmp');
% he_houghgau_er_Sobel = re_houghgau_er_Sobel;
% he_houghgau_er_Sobel(houghgau==255)=255;
he_houghgau_er_Sobel = houghgau;
he_houghgau_er_Sobel(re_houghgau_er_Sobel==255)=255;
imwrite(he_houghgau_er_Sobel,'he_houghgau_er_Sobel.bmp');
  

[circle_xy,circle_r] = Hough(houghsalt_er_Sobel)
imwrite(ht(circle_xy,circle_r), 're_houghsalt_er_Sobel.bmp')
re_houghsalt_er_Sobel = imread('re_houghsalt_er_Sobel.bmp');
% he_houghsalt_er_Sobel = re_houghsalt_er_Sobel;
% he_houghsalt_er_Sobel(houghsalt==255)=255;
he_houghsalt_er_Sobel = houghsalt;
he_houghsalt_er_Sobel(re_houghsalt_er_Sobel==255)=255;
imwrite(he_houghsalt_er_Sobel,'he_houghsalt_er_Sobel.bmp');


[circle_xy,circle_r] = Hough(houghgau_er_Laplacian)
imwrite(ht(circle_xy,circle_r), 're_houghgau_er_Laplacian.bmp')
re_houghgau_er_Laplacian = imread('re_houghgau_er_Laplacian.bmp');
% he_houghgau_er_Laplacian = re_houghgau_er_Laplacian;
% he_houghgau_er_Laplacian(houghgau==255)=255;
he_houghgau_er_Laplacian = houghgau;
he_houghgau_er_Laplacian(re_houghgau_er_Laplacian==255)=255;
imwrite(he_houghgau_er_Laplacian,'he_houghgau_er_Laplacian.bmp');


[circle_xy,circle_r] = Hough(houghsalt_er_Laplacian)
imwrite(ht(circle_xy,circle_r), 're_houghsalt_er_Laplacian.bmp')
re_houghsalt_er_Laplacian = imread('re_houghsalt_er_Laplacian.bmp');
% he_houghsalt_er_Laplacian = re_houghsalt_er_Laplacian;
% he_houghsalt_er_Laplacian(houghsalt==255)=255;
he_houghsalt_er_Laplacian = houghsalt;
he_houghsalt_er_Laplacian(re_houghsalt_er_Laplacian==255)=255;
imwrite(he_houghsalt_er_Laplacian,'he_houghsalt_er_Laplacian.bmp');

Hough：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
function [ par1, par3 ] = Hough( BW )

r_max=100;
r_min=40;
step_r=1;
step_angle=pi/20;
p=0.5;
[m,n] = size(BW);
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
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;
    par(:,k) = [par1,par2,par3];
end

par1=[par2 par1];

end

ht:

1
2
3
4
5
6
7
8
9
10
11
12
13
function hb = ht( par1, par3 )
par1_y = par1(1);
par1_x = par1(2);
hb = zeros(325, 440);
k = 1;
for i = 1:325
     for j = 1:440
         if ((abs(i-par1_x)^2 + abs(j-par1_y)^2) >= (par3-k)^2) && ((abs(i-par1_x)^2 + abs(j-par1_y)^2) <= (par3+k)^2)
             hb(i, j) = 255;
         end
     end
end
hb = uint8(hb);