Python計算機視覺——Harris角點檢測器

Harris角點檢測器

1 知識點

角點定義：

a. 區域視窗沿各方向移動，均產生明顯變化的點

b. 影像區域曲率突變的點

角點圖釋：
常見的角點：

相關名詞：

兩個像素塊（大小相同）： I 1 ( x ) I_1(x) I1?(x), I 2 ( x ) I_2(x) I2?(x)
視窗函式： Σ x , y w ( x , y ) \Sigma_{x,y}w(x,y) Σx,y?w(x,y)
灰度變化： E ( u , v ) E(u,v) E(u,v)
角點回應函式： R R R
權重矩陣： W W W(通常為高斯濾波器 G σ G_\sigma Gσ?)

相關計算：

1. 將視窗平移 [ u , v ] [u,v] [u,v]產生灰度變化 E ( u , v ) E(u,v) E(u,v)

   E ( u , v ) = Σ x , y w ( x , y ) [ I ( x + u , y + v ) ? I ( x , y ) ] E(u,v)=\Sigma_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)] E(u,v)=Σx,y?w(x,y)[I(x+u,y+v)?I(x,y)]

   利用一級泰勒展開式，可得，

   I ( x + u , y + v ) = I ( x , y ) + I x u + I y v + O ( u 2 , v 2 ) I(x+u,y+v)=I(x,y)+I_xu+I_yv+O(u^2,v^2) I(x+u,y+v)=I(x,y)+Ix?u+Iy?v+O(u2,v2)

   E ( u , v ) ≈ Σ x , y w ( x , y ) [ I x u + I y v ] 2 ? [ u v ] M [ u v ] E(u,v)≈\Sigma_{x,y}w(x,y)[I_xu+I_yv]^2\cong \begin{bmatrix}u \ v\end{bmatrix}M \begin{bmatrix}u \\ v\end{bmatrix} E(u,v)≈Σx,y?w(x,y)[Ix?u+Iy?v]2?[u v?]M[uv?]
   其中，

   M = Σ x , y w ( x , y ) [ I x 2 I x I y I x I y I y 2 ] M=\Sigma_{x,y}w(x,y)\begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \\ \end{bmatrix} M=Σx,y?w(x,y)[Ix2?Ix?Iy??Ix?Iy?Iy2??]

2. 點 x x x的對稱半正定矩陣:
   M I = M I ( x ) = ? I ? I T = [ I x I y ] [ I x I y ] = [ I x 2 I x I y I x I y I y 2 ] M_I=M_I(x)=\nabla I \nabla I^T=\begin{bmatrix} I_x \\ I_y \\ \end{bmatrix} \begin{bmatrix} I_x \ I_y\\ \end{bmatrix} =\begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \\ \end{bmatrix} MI?=MI?(x)=?I?IT=[Ix?Iy??][Ix? Iy??]=[Ix2?Ix?Iy??Ix?Iy?Iy2??]
   M I M_I MI?周圍像素上的區域平均為：
   M I ˉ = W ? M I \bar{M_I}=W*M_I MI?ˉ?=W?MI?
   通過式 λ E ? M I = 0 \lambda E-M_I=0 λE?MI?=0，得到特征值 λ 1 = ∣ ? I ∣ 2 \lambda_1=|\nabla I|^2 λ1?=∣?I∣2， λ 2 = 0 \lambda_2=0 λ2?=0
   特征值會依賴于區域影像變化，具體表現為：影像梯度變化 ?? ? ?? \implies ? M I ˉ \bar{M_I} MI?ˉ?特征值變化，

不同的特征值對應不同的像素情況：
a. λ 1 \lambda_1 λ1?、 λ 2 \lambda_2 λ2?都是很大的正數 ?? ? ?? \implies ? x x x為角點
b. λ 1 \lambda_1 λ1?很大， λ 2 = 0 \lambda_2=0 λ2?=0 ?? ? ?? \implies ?區域存在一個邊， M I ˉ \bar{M_I} MI?ˉ?變化不大
c. λ 1 = λ 2 = 0 \lambda_1= \lambda_2=0 λ1?=λ2?=0 ?? ? ?? \implies ?區域為空

3. 角點回應
   R = d e t M I ? k ( t r a c e M I ) 2 R = det M_I-k(trace M_I)^2 R=detMI??k(traceMI?)2或
   R = d e t M I ( t r a c e M I ) 2 R = \frac{det M_I}{(trace M_I)^2} R=(traceMI?)2detMI??
   其中，
   d e t M I = λ 1 λ 2 det M_I = \lambda_1 \lambda_2 detMI?=λ1?λ2?
   t r a c e M I = λ 1 + λ 2 trace M_I = \lambda_1+\lambda_2 traceMI?=λ1?+λ2?

R 只與M的特征值有關
角點：R 為大數值正數
邊緣：R為大數值負數
平坦區：R為小數值

4. 角點進行挑選方法：
   a. 設定閾值，選取像素值高于閾值的所有影像點
   b. 角點之間間隔必須大于設定的最小距離

5. 在影像中尋找對應點：
   Harris角點檢測器僅能夠檢測出影像中的興趣點，但沒有給出通過比較影像間的興趣點來尋找匹配角點的方法，
   Harris角點的描述子通常是由周圍影像像素塊的灰度值以及用于比較的歸一化互相關矩陣構成的，
   兩個像素塊 I 1 ( x ) I_1(x) I1?(x), I 2 ( x ) I_2(x) I2?(x)的相關矩陣為
   c ( I 1 , I 2 ) = Σ x f ( I 1 ( x ) , I 2 ( x ) ) , f 指 的 是 相 關 方 法 c(I_1,I_2)=\Sigma_{x}f(I_1(x),I_2(x)),f指的是相關方法 c(I1?,I2?)=Σx?f(I1?(x),I2?(x)),f指的是相關方法
   互相關矩陣為
   f ( I 1 , I 2 ) = I 1 I 2 f(I_1,I_2)=I_1I_2 f(I1?,I2?)=I1?I2?
   ∴ c ( I 1 , I 2 ) = I 1 ? I 2 \therefore c(I_1,I_2)=I_1·I_2 ∴c(I1?,I2?)=I1??I2?
   ? · ? 指向量乘法， c ( I 1 , I 2 ) c(I_1,I_2) c(I1?,I2?)值越大， I 1 , I 2 I_1,I_2 I1?,I2?相似度越高，
   歸一化互相關矩陣為：
   n c c ( I 1 , I 2 ) = 1 n ? 1 Σ x ( I 1 ( x ) ? μ 1 ) σ 1 ? ( I 2 ( x ) ? μ 2 ) σ 2 ncc(I_1,I_2)=\frac{1}{n-1}\Sigma_{x}\frac{(I_1(x)-\mu_1)}{\sigma_1}·\frac{(I_2(x)-\mu_2)}{\sigma_2} ncc(I1?,I2?)=n?11?Σx?σ1?(I1?(x)?μ1?)??σ2?(I2?(x)?μ2?)?
   n n n指像素塊中像素的數目，像素塊即以該像素點為中心的周圍矩形部分影像
   μ \mu μ指平均像素強度
   σ \sigma σ指像素塊的標準差

2 實驗部分

1. 角點檢測
   a. 回傳角點回應值 R R R

def compute_harris_response(im,sigma=3):
    """ Compute the Harris corner detector response function 
        for each pixel in a graylevel image. """
    
    # derivatives
    imx = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    
    # compute components of the Harris matrix
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)
    
    # determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy
    
    return Wdet / Wtr

b. 回傳角點并標記

def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    """ Return corners from a Harris response image
        min_dist is the minimum number of pixels separating 
        corners and image boundary. """
    
    # find top corner candidates above a threshold
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    
    # get coordinates of candidates
    coords = array(harrisim_t.nonzero()).T
    
    # ...and their values
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    
    # sort candidates (reverse to get descending order)
    index = argsort(candidate_values)[::-1]
    
    # store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    
    # select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist), 
                        (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    
    return filtered_coords
def plot_harris_points(image,filtered_coords):
    """ Plots corners found in image. """
    
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],
                [p[0] for p in filtered_coords],'ro')
    axis('off')
    show()

c. 運行主函式

# 讀入影像
im = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert('L'))
# 檢測harris角點
harrisim = compute_harris_response(im)
# Harris回應函式
harrisim1 = 255 - harrisim

figure()
gray()

#畫出Harris回應圖
subplot(141)
imshow(harrisim1)
axis('off')
axis('equal')

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    filtered_coords = get_harris_points(harrisim, 6, thres)
    subplot(1, 4, i+2)
    imshow(im)
    plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '.')
    axis('off')
show()

d. 運行結果圖

域值越大，得到的角點越稀疏

也可以直接呼叫函式得到角點，結果圖為：

代碼為：

im = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert('L'))
harrisim = compute_harris_response(im)
filter_coords = get_harris_points(harrisim,6)
plot_harris_points(im, filtered_coords)

2. 影像中找對應點
   a. 選用奇數像素的灰度影像塊

def get_descriptors(image,filtered_coords,wid=5):
    """ For each point return pixel values around the point
        using a neighbourhood of width 2*wid+1. (Assume points are 
        extracted with min_distance > wid). """
    
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1,
                            coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)
    
    return desc

b. 歸一化互相關

def match(desc1,desc2,threshold=0.5):
    """ For each corner point descriptor in the first image, 
        select its match to second image using
        normalized cross correlation. """
    
    n = len(desc1[0])
    
    # pair-wise distances
    d = -ones((len(desc1),len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
            d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1) 
            if ncc_value > threshold:
                d[i,j] = ncc_value
            
    ndx = argsort(-d)
    matchscores = ndx[:,0]
    
    return matchscores

c. 兩邊對稱匹配

def match_twosided(desc1,desc2,threshold=0.5):
    """ Two-sided symmetric version of match(). """
    
    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)
    
    ndx_12 = where(matches_12 >= 0)[0]
    
    # remove matches that are not symmetric
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1
    
    return matches_12

d. 兩幅影像拼接

def appendimages(im1,im2):
    """ Return a new image that appends the two images side-by-side. """
    
    # select the image with the fewest rows and fill in enough empty rows
    rows1 = im1.shape[0]    
    rows2 = im2.shape[0]
    
    if rows1 < rows2:
        im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
    elif rows1 > rows2:
        im2 = concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)
    # if none of these cases they are equal, no filling needed.
    
    return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
    """ Show a figure with lines joining the accepted matches 
        input: im1,im2 (images as arrays), locs1,locs2 (feature locations), 
        matchscores (as output from 'match()'), 
        show_below (if images should be shown below matches). """
    
    im3 = appendimages(im1,im2)
    if show_below:
        im3 = vstack((im3,im3))
    
    imshow(im3)
    
    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0:
            plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    axis('off')

e. 運行結果：

im1 = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert("L"))
im2 = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert("L"))

wid = 5
harrisim = compute_harris_response(im1, 5)
filtered_coords1 = get_harris_points(harrisim, wid+1)
d1 = get_descriptors(im1, filtered_coords1, wid)

harrisim = compute_harris_response(im2, 5)
filtered_coords2 = get_harris_points(harrisim, wid+1)
d2 = get_descriptors(im2, filtered_coords2, wid)

print('starting matching')
matches = match_twosided(d1, d2)

figure()
gray() 
plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
show()

圖形像素塊的互相關矩陣有較弱的描述性，Harris不具有尺度不變性和旋轉不變性，演算法中的像素塊的大小也會影響對應匹配的結果，