matlab-霍夫變換詳解（判斷正方形長方形）

霍夫變換

霍夫變換是1972年提出來的，最開始就是用來在圖像中過檢測直線，后來擴展能檢測圓、曲線等。

直線的霍夫變換就是 把xy空間的直線 換成成 另一空間的點。就是直線和點的互換。

 

我們在初中數學中了解到，一條直線可以用如下的方程來表示：y=kx+b，k是直線的斜率，b是截距。

我們轉換下變成：b=-kx+y。我們是不是也可以把（k，b）看作另外一個空間中的點？這就是k-b參數空間。 這樣，我們就把一條x-y直線用一個（k，b）的點表示出來了。

我們看到，在x-y圖像空間中的一個點，變成了k-b參數空間中的一條直線，而x-y圖像空間中的2點連成的直線，變成了k-b參數空間中的一個交點。

如果x-y圖像空間中有很多點在k-b空間中相交於一點，那么這個交點就是我們要檢測的直線。這就是霍夫變換檢測直線的基本原理。

當然，有一個問題需要注意，圖像空間中如果一條直線是垂直的，那么斜率k是沒有定義的（或者說無窮大）。為了避免這個問題，霍夫變換采用了另一個參數空間：距離-角度參數空間。也就是極坐標。

我們在中學中學過，平面上的一個點也可以用距離-角度來定義，也就是極坐標。那么在圖像中，每一個點都可以用距離和角度來表達：

 

但是，使用距離-角度后，點（x，y）與距離，角度的關系變成了：

ρ=xcosθ+ysinθ

於是，在新的距離-角度參數空間中，圖像中的一個點變成了一個正弦曲線，而不是k-b參數空間中的直線了。這些正弦曲線的交點就是圖像空間中我們要檢測的直線了。

 

Matlab霍夫變換的函數詳解

[H, theta, rho] = hough(BW,ParameterName, ParameterValue)

BW：二值圖

ParameterName：'RhoResolution'或'Theta'

RhoResolution-指定在累計數組中（檢測極值）的檢測間隔？默認為1

Theta-指定檢測的角度范圍（不超過-90~90度）以及間隔，例如-90:0.5:89.5，默認-90:1:89

H：累計數組

Theta：H對應的θ，實際上H的大小就是Rho×Theta

Rho：H對應的ρ

 這兩個參數值的注意，RhoResolution太大覆蓋不到極值點，檢測到一些不對應直線的次極值，

 

峰值提取

peaks = houghpeaks(H,numpeaks)

peaks = houghpeaks(...,param1, val1, param2, val2)

H：累計數組；

Numpeaks：指定需要檢測的峰值個數；

Param1：可以是'Threshold'或'NHoodSize'

'Threshold'-指定峰值的域值，默認是0.5*max(H(:))

'NHoodSize'-是個二維向量[m,n]，檢測到一個峰值后，將峰值周圍[m,n]內元素置零。

 

畫直線段

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

lines = houghlines(...,param1, val1, param2, val2)

 

BW：二值圖

Theta、rho、peaks：分別來自函數hough和houghpeaks

Lines：結構數組，大小等於檢測到的直線段數，每個單元包含

Point1、point2：線段的端點

Theta、rho：線段的theta和rho

 

下面實例演示

用霍夫變換判斷矩形

所用圖形是我用畫圖工具制作的，圖片很干凈，所以不需要濾波等措施。以下是圖案。

img  = imread('a.png');
img_gray = rgb2gray(img);

%背景是黑的！！！
threshold =graythresh(img_gray);%取閾值
bw=im2bw(img_gray,threshold);
if length(bw(bw==1))>length(bw(bw==0))
    bw=~bw;
end
%填充
bw = imfill(bw,'holes');

figure(), imshow(img_gray), title('image');



[B,L] = bwboundaries(bw,'noholes');
stats = regionprops(L,'Area','Centroid','image');
boundary = B{1};
delta_sq = diff(boundary).^2;
perimeter = sum(sqrt(sum(delta_sq,2)));
% obtain the area calculation corresponding 
area = stats(1).Area;
% compute the roundness metric
metric = 4*pi*area/perimeter^2;
if metric>0.9
    disp('yuan')
else
    
    % the canny edge of image
    BW = edge(bw,'canny');
    % Iedge=edge(Ihsv,'sobel');    %邊沿檢測
    %BW = imdilate(BW,ones(3));%圖像膨脹
    figure(), imshow(BW), title('image edge');
    % the theta and rho of transformed space
    [H,Theta,Rho] = hough(BW,'RhoResolution',1.2);
    figure(), imshow(H,[],'XData',Theta,'YData',Rho,'InitialMagnification','fit'),...
        title('rho\_theta space and peaks');
    xlabel('\theta'), ylabel('\rho');
    axis on, axis normal, hold on;
    % label the top 5 intersections
    P  = houghpeaks(H,4,'threshold',ceil(0.1*max(H(:))),'NHoodSize',[35,11]);
    x = Theta(P(:,2)); 
    y = Rho(P(:,1));
    plot(x,y,'*','color','r');


    % find lines and plot them
    lines = houghlines(BW,Theta,Rho,P);
    figure(), imshow(img),title('final image');
    hold on
    b_len = ones(1,length(lines));
    for k = 1:length(lines)
     xy = [lines(k).point1; lines(k).point2];
     b_len(1,k)=sqrt((xy(1,1)-xy(2,1))^2+(xy(1,2)-xy(2,2))^2);
     plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','r');
     hold on
    end
    if var(b_len)<10
        disp('正方形')
    else
        disp('長方形')
    end
end

　然后得到結果

變換參數域圖像（很完美，4個點）

結果圖像

基本包含在原圖形上，測試成功！

 

參考文獻：

[1] http://blog.sina.com.cn/s/blog_ac7218750101giyf.html

 

[2] https://blog.csdn.net/saltriver/article/details/80547245

 

[3] https://ww2.mathworks.cn/help/images/ref/hough.html

 

盤外篇

自己寫個Hough變換？ hough的原理很簡單，可以自己嘗試以下。

那么我們檢測一下 直線的斜率。

下面上代碼

clear all;

img  = imread('demo1.png');
img_gray = rgb2gray(img);

%二值化
threshold =graythresh(img_gray);%取閾值
bw=im2bw(img_gray,threshold);
if length(bw(bw==1))>length(bw(bw==0))
    bw=~bw;
end

figure(), imshow(bw), title('image');
% %填充
% bw = imfill(bw,'holes');
% %邊沿檢測
%the canny edge of image
%bw = edge(bw,'canny');

bw = imdilate(bw,ones(8)); %圖像膨脹(因為直線太細會使誤差增大)
%figure(), imshow(bw), title('image ');
%the theta and rho of transformed space
[H,Theta,Rho] = hough1(bw);
figure(), imshow(H,[],'XData',Theta,'YData',Rho,'InitialMagnification','fit'),...
    title('rho\_theta space and peaks');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
%label the top 1 intersections
%P = houghpeaks(H,1);
[xx,yy]=find(H==max(max(H)),1);
P=[xx,yy];
x = Theta(P(:,2));
y = Rho(P(:,1));
plot(x,y,'*','color','r');

[m,n]=size(bw);
label=zeros(m,n);
xy=[];
for ii=1:length(bw)
    jj=ceil((y-ii*cos(x/180*pi))/sin(x/180*pi));
    if jj>=1&& jj<=n
        label(ii,jj)=1;
        if bw(ii,jj)==1
            xy=[xy;jj,ii];  %因為plot和imshow的xy軸是鏡像關系！
        end            
    end
end
    
    
%find lines and plot them
%lines = houghlines(bw,Theta,Rho,P);
figure(), imshow(img),title('final image');
hold on
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','r');

k=-sin(x/180*pi)/cos(x/180*pi);
b=y/cos(x/180*pi);
disp(['斜率是',num2str(k)])
disp(['y軸截距是',num2str(b)])

%ezplot('x*cos(-72/180*pi)+y*sin(-72/180*pi)=-31',[0,256])

　　下面就是自己寫的hough變換函數了

function [h, theta, rho] = hough1(f, dtheta, drho)

if  nargin < 3
   drho = 1;
end
if  nargin < 2
   dtheta = 1;
end
 
f = double(f);
[M,N] = size(f);
theta = linspace(-90, 0, ceil(90/dtheta) + 1);
theta = [theta -fliplr(theta(2:end - 1))];
ntheta = length(theta);
 
D = sqrt((M - 1)^2 + (N - 1)^2);
q = ceil(D/drho);
nrho = 2*q - 1;
rho = linspace(-q*drho, q*drho, nrho);
 
[x, y, val] = find(f);
x = x - 1;  y = y - 1;
 
h = zeros(nrho, length(theta));
 
for  k = 1:ceil(length(val)/1000)
   first = (k - 1)*1000 + 1;
   last  = min(first+999, length(x));
   
   x_matrix     = repmat(x(first:last), 1, ntheta);
   y_matrix     = repmat(y(first:last), 1, ntheta);
   val_matrix   = repmat(val(first:last), 1, ntheta);
   theta_matrix = repmat(theta, size(x_matrix, 1), 1)*pi/180;
   
   rho_matrix = x_matrix.*cos(theta_matrix) + ...
       y_matrix.*sin(theta_matrix);
   slope = (nrho - 1)/(rho(end) - rho(1));
   rho_bin_index = round(slope*(rho_matrix - rho(1)) + 1);
   
   theta_bin_index = repmat(1:ntheta, size(x_matrix, 1), 1);
   
   h = h + full(sparse(rho_bin_index(:), theta_bin_index(:), ...
                       val_matrix(:), nrho, ntheta));
end

　　結果還不錯