hough函数编写matlab,hough图像处理并非使用matlab自带函数（容易理解）

%%能够运行close all;

clc;

I=imread('C:\Users\Administrator.GHTT-20141012XI\Desktop\lena.jpg');

[x,y]=size(I);

BW1= dither(I);

BW=edge(BW1,'canny');

figure;imshow(I);title('原图')

%figure;imshow(BW);title('边缘检测图像')

rho_max=floor(sqrt(x^2+y^2))+1; %由原图数组坐标算出ρ最大值，并取整数部分加1

%此值作为ρ，θ坐标系ρ最大值

accarray=zeros(rho_max,180); %定义ρ，θ坐标系的数组，初值为0。

%θ的最大值，180度

Theta=[0:pi/180:pi]; %定义θ数组，确定θ取值范围

for n=1:x,

for m=1:y

if BW(n,m)==1

for k=1:180

%将θ值代入hough变换方程，求ρ值

rho=(n*cos(Theta(k)))+(m*sin(Theta(k)));

%将ρ值与ρ最大值的和的一半作为ρ的坐标值(数组坐标)，这样做是为了防止ρ值出现负数

rho_int=round(rho/2+rho_max/2);

%在ρθ坐标(数组)中标识点，即计数累加

accarray(rho_int,k)=accarray(rho_int,k)+1;

end

end

end

end

%figure;colormap gray;

%imagesc(accarray);title('hough变换后的图')

%xlabel('\theta'), ylabel('\rho');

% %=======利用hough变换提取直线======%

% %寻找100个像素以上的直线在hough变换后形成的点

K=1; %存储数组计数器

% case_accarray_n = zero(1000)

% case_accarray_m = zero(1000)

for rho_n=1:rho_max %在hough变换后的数组中搜索

for theta_m=1:180

if accarray(rho_n,theta_m)>100 %设定直线的最小值。

case_accarray_n(K)=rho_n; %存储搜索出的数组下标

case_accarray_m(K)=theta_m;

K=K+1;

end

end

end

% %把这些点构成的直线提取出来,输出图像数组为I_out

I_out=zeros(x,y);

for n=1:x,

for m=1:y

if BW(n,m)==1

for k=1:180

rho=(n*cos(Theta(k)))+(m*sin(Theta(k)));

rho_int=round(rho/2+rho_max/2);

%如果正在计算的点属于100像素以上点，则把它提取出来

for a=1:K-1

if rho_int==case_accarray_n(a)&k==case_accarray_m(a)%%%==gai==%%% k==case_accarray_m(a)&rho_int==case_accarray_n(a)

I_out(n,m)=BW(n,m);

end

end

end

end

end

end

figure;

imshow(I_out);title('利用经典hough变换提取的图像');

%%并不是hough检测直线

I=imread('C:\Users\Administrator.GHTT-20141012XI\Desktop\roadline.png');

I=rgb2gray(I);%彩色图象 to 二值图像

[m,n]=size(I); %获取维数

I=im2double(I); %便于计算 转成double

im=zeros(m,n);

%不懂这算子干嘛

for i=3:m-2

for j=3:n-2%处理领域较大，所以从图像(3,3)开始，在(m-2,n-2)结束

l(i,j)=-I(i-2,j)-I(i-1,j-1)-2*I(i-1,j)-I(i-1,j+1)-I(i,j-2)-2*I(i,j-1)+16*I(i,j)-2*I(i,j+1)-I(i,j+2)-I(i+1,j-1)-2*I(i+1,j)-I(i+1,j+1)-I(i+2,j);%LoG算子

end

end

[m,n]=size(l);

for i=2:m-1

for j=2:n-1

y(i,j)=l(i-1,j-1)+l(i-1,j)+l(i-1,j+1)+l(i,j-1)+l(i,j)+l(i,j+1)+l(i+1,j-1)+l(i+1,j)+l(i+1,j+1);

y(i,j)=y(i,j)/9; %LoG算子提取边缘后，对结果进行均值滤波以去除噪声，为下一步hough变换提取直线作准备

end

end

%换成整形

q=im2uint8(y);

[m,n]=size(q);

for i=1:m

for j=1:n

if q(i,j)>70; %设置二值化的阈值为70

q(i,j)=255; %对图像进行二值化处理，使图像边缘更加突出清晰

else

q(i,j)=0;

end

end

end

%Hough变换检测直线，使用(a，p)参数空间，a∈[0,180],p∈[0,2d]

a=180; %角度的值为0到180度

d=round(sqrt(m^2+n^2)); %图像对角线长度为p的最大值

s=zeros(a,2*d); %存储每个(a,p)个数

z=cell(a,2*d); %用元胞存储每个被检测的点的坐标

for i=1:m

for j=1:n%遍历图像每个点

if(q(i,j)==255)%只检测图像边缘的白点，其余点不检测

for k=1:a

p = round(i*cos(pi*k/180)+j*sin(pi*k/180));%对每个点从1到180度遍历一遍，取得经过该点的所有直线的p值(取整)

if(p > 0)%若p大于0，则将点存储在(d，2d)空间

s(k,d+p)=s(k,d+p)+1;%(a，p)相应的累加器单元加1

z{k,d+p}=[z{k,d+p} ,[i,j]];%存储点坐标

else

ap=abs(p)+1;%若p小于0，则将点存储在(0，d)空间

s(k,ap)=s(k,ap)+1;%(a，p)相应的累加器单元加1

z{k,ap}=[z{k,ap},[i,j]];%存储点坐标

end

end

end

end

end

for i=1:a

for j=1:d*2 %检查每个累加器单元中存储数量

if(s(i,j) >30) %将提取直线的阈值设为70

lp=z{i,j};%提取对应点坐标

for k=1:s(i,j)%对满足阈值条件的累加器单元中(a，p)对应的所有点进行操作

im(lp(1,k),lp(2,k))=255; %每个点R分量=255，G分量=0，B分量=0

end

end

end

end

subplot(1,2,1),imshow(I);title('原始图像');

subplot(1,2,2),imshow(im);title('hough变换边缘图像');