第八單元 Hough轉換

陳慶瀚

2004-11-15

 

Hough transform (HT) 於1962年由P.Hough (P.V.C. Hough. Method and Means for Recognizing Complex Patterns, US Patent 3,069,654, December 1962)提出， 繼而於1972年首度以論文形式發表(R.O.Duda, R.E.Hart. Use of the Hough Transform to Detect Lines and Curves in Pictures, CACM(15). No. 1, January 1972, pp. 11-15.)。

HT用在二值化影像的形狀偵測(shape detection)，主要原理是利用影像中分散的點位置找出特定形狀(例如直線或圓)的參數值，每一個點藉由一對多的映射(由影像空間映射到參數空間)產生參數的所有可能值，再累計全部點所產生的參數值，最後在得以在參數空間決定表現最明顯的形狀參數。

P.V.C. Hough. Method and Means for Recognizing Complex Patterns, US Patent 3,069,654, December 1962.

R.O.Duda, R.E.Hart. Use of the Hough Transform to Detect Lines and Curves in Pictures, CACM(15). No. 1, January 1972, pp. 11-15.

 

1.直線偵測(Straight Lines )的Hough轉換

對於二值化影像上的任一點(x,y)，通過這一點的直線方程式可以表示為

f(x,y) = y-ax-b = 0

a和b分別是直線的斜率和截距。

我們可以把上式視為相互限制條件(mutual constraint)的映射關係，由影像點(x,y)映射至多重的參數(a,b)，或由參數(a,b)映射至多重的影像點(x,y)。也就是說，影像空間上的一個點(x,y)，可以定義參數空間(a,b)上一條直線的多重(理論上無限多)點；同樣情形，參數空間的一點(a,b)，也可以定義出影像空間上一條直線的眾多點(x,y)，如下圖。


影像空間  x-y              參數空間a-b

 

使用累增器(Accumulator)

 

由於Hough transform把每一個影像點(x,y)映射至多重的參數點(a,b)，我們使用一個累增器來紀錄每一組(a,b)出現的次數，出現頻率最高的一組(a,b)，就是影像空間上最具代表性的一條直線。例如

Hough transform演算法：

1. 在影像空間上決定所有可能的特徵點，通常是邊緣點或骨架點。.

2. 尋找影像空間的一個特徵點(x,y)，計算

  2.1 對於中每一個a，計算通過(x,y)的所有直線的(a’,b’)

2.2 在accumulator的(a’,b’)位置上累加一次

尋找下一個特徵點，重複步驟2.1和2.2，直到所有特徵點都通過計算。

3. 找到accumulator的所有區域最大值(maxima) 。

4. 將每一個maxima映射回影像空間的每一條代表直線。

 

使用極座標(Polar system)

 

從實作的角度，使用f(x,y) = y-ax-b = 0並不可行，因為當影像空間上的直線是一垂直線時，則斜率a=∞，而無法納入有限的accumulator。

 

所以，實用的Hough transform都採用極座標(ρ, θ)來取代(a, b)，極座標與直角座標的轉換公式為：

x cos(θ) + y sin(θ) = r

影像空間上的每一條直線被映射到(ρ, θ)參數空間的一個點，而影像空間上每一個點則被映射到(ρ, θ)參數空間的一條曲線，如下圖：

如果有影像空間同一直線上的三個點，則(ρ, θ)參數空間的三條曲線將交於一點，也就是accumulator在此一位置的累加值為3。

2. Hough轉換程式設計

#include <fstream.h>

#include <math.h>

#include "array.h"

#define feature_point 0

#define background 255

 

class HoughLineDetector

{

   float sin_tab[180],cos_tab[180];

public:

    HoughLineDetector()// Generate sin and cos look-up table

    {

        float theta;

        for(int angle=0;angle<180;angle++)

        {

             theta=angle*3.14159265358979/180.0;

             cos_tab[angle]=cos(theta);

             sin_tab[angle]=sin(theta);

        }   

    }

   void Transform(uc2D &im,uc2D &hough);

};

 

void  HoughLineDetector::Transform(uc2D &im, uc2D &hough_ima)

{

   float rmax=sqrt((float)im.nr*(float)im.nr+

(float)im.nc*(float)im.nc);

   hough_ima.Initialize(2*rmax,180);

   for(int j=0;j<hough_ima.nr;j++)

for(int i=0;i<hough_ima.nc;i++)

hough_ima.m[j][i]=0;

   float r;

   int nr,angle;

   for(int i=0;i<im.nr;i++)for(int j=0;j<im.nc;j++)

   {

      if(im.m[j][i]==feature_point)

      {

           for(angle=0;angle<180;angle++)  // for each theta

           {

               r=i*sin_tab[angle]+j*cos_tab[angle];// compute r

               nr=int(r+rmax);                     //shift r to positive value

               hough_ima.m[nr][angle]++;           // accumulation

           }

      }

   }

}

 

程式執行範例如下：

 

 

 

 

由於Hough轉換是一個運算耗時、且記憶體需求龐大的計算方法，1986年Lavin和Le Master提出一個快速Hough轉換(FHT)方法 [H. Li, M.A. Lavin, and R.J. Le Master. Fast hough transform: A hierarchical approach. Computer Vision, Graphics and Image Processing, 36:139-161, 1986.]

 

習題1、請在HoughLineDetector物件類別加入兩個成員函式：

void FindLocalMaximum(uc2D &hough)

功能：找出hough space的所有local maximum，並紀錄在兩個rmax, thetamax陣列(宣告為物件類別的i1D的成員變數)中。

void backProjection(uc2D &ima)

功能：使用rmax, thetamax的資訊，繪出對應x-y image sapce的直線。將miat400x200.raw影像所抽取的骨架影像進行Hough轉換，輸出Hough image，再求出所有local maximum，最後以backProjection()將所偵測出來的直線繪在x-y影像上。

3. 圓形偵測的Hough Transform

一個圓可以用下面三個參數來描述：圓心座標(a,b)和半徑r：

仿照直線偵測的方法，可以將影像空間上的一點(x,y)映射至三維的參數空間(a,b,r)上的圓錐(circular cone)：

 

如果影像空間上的兩個點位於同一個圓的軌跡，則參數空間上的兩個圓錐將交於一點(a’ ,b’, r’)。

 

累增器(Accumulator)的使用如同直線偵測的Hough轉換，以下程式片段即為圓形偵測的Hough Transform演算法：

 

for(i=0; i<nr; i++)for (j=0; j<nc; j++)

{

   if (edge_image[i] > 0) // This is an edge pixel

   {

      x = j;

      y = i;

   // populate the accumulator array

      for(a= 0;a<max_a;a++)for(b=0;b<max_b;b++)

//For each circle center,execute:

      {

          r = (int) sqrt((x-a)*(x-a) + (y-b)*(y-b));

          if ((r>0)&&(r<max_r))accumulator[r][a][b]++;

      }

   }

}

 

習題2、請仿照HoughLineDetector程式寫一個HoughCircleDetector物件類別，並應用在multicoins240x320.raw，偵測出影像中的四枚錢幣。

 

4. 泛用型Hough轉換(Generalized Hough Transform)

為了得以偵測無法被參數化(沒有解析模型)的不規則形狀，Ballard首先提出泛化Hough轉換(Generalized Hough Transform，GHT)[D.H. Ballard, “Generalized the Hough Transform to Detect Arbitrary Shapes,” Pattern Recognition, Vol. 3, No. 2, pp. 110-122, 1981.]

 

定義任意形狀(generalized shape)的表示參數:

a = {y, s, θ}

其中， y=(xr, yr)表示一個shape在影像中的參考原點， θ是shape方向(orientation)， s=(sx,sy)是兩個正交(x軸和y軸方向)的比率參數(scale factors)。

 

上面的表示式共有5個參數。有時為了簡化計算複雜度，我們將s參數視為純量(scalar), 使得參數空間剩下4維。如果考慮一個半徑r0的圓形，由於半徑固定，同時其方向性參數θ也是固定，因此此一表示式仍符合圓形偵測的Hough表示式。

 

對於任一形狀的邊界點x位置的梯度(gradient)方向與x-軸的夾角為Φ，

 

泛用型Hough轉換演算法分成兩個階段，分別是樣版(template)建立階段和形狀偵測階段。

I.樣版建立階段演算法：建立R-Table

Step 1、決定位於樣版形狀內部的一個參考點─例如重心點(xc, yc)；

Step 2、建立一個空的R-Table，其索引值為Φi,i=1,2,K,其角度增量為p/K，亦即Φ從0變化到180度，角度增量表示梯度方向的解析度；

Step 3、針對每一個邊緣點(x,y)，計算(r,α)參數:

r = sqrt((x-xc)2+(y-yc)2)

及

α = tan-1((y-yc)/(x-xc))

Step 4、計算f (垂直梯度方向)，並將(r,α)加入與f值最接近的Φi；

Step 5、重複Step 4和5，直到所有的邊緣點都完成測試，得到R-table:

II. 形狀偵測階段

Step 1、建立一個2D Hough table H(xc, yc)，將其初值全設為0；

Step2、針對每一個邊緣點(x,y)，計算垂直梯度方向的夾角f’；

Step 3、在R-Table中，找出最接近f’的fi，對於fi內存放的所有(ri,αi)，計算

Step 4、將H(xc, yc)累加1，重複步驟2和3，直到所有的邊緣點都完成測試。

Step 5、找出H(xc, yc)中的區域最大值，其(xc, yc)即為偵測形狀的位置。

5. 考慮旋轉和縮放的任意形狀偵測

原理與上一節相同，但是我們加上旋轉參數q和縮放參數S, Hough table因此變成H(xc, yc,q,S)。至於偵測演算法的step 3變成

xc = x + rScos(α+q)

yc = y + rSsin(α+q)

其餘步驟皆同，但計算複雜度則大幅增加。

 

習題3(option)、請以下圖(F-16戰機)為對象，設計一個4D的GHT，可偵測旋轉的、外型大小可變的F-16形狀。