Python计算机视觉——Harris角点检测器

Harris角点检测器

1 知识点

角点定义：

a. 局部窗口沿各方向移动，均产生明显变化的点

b. 图像局部曲率突变的点

角点图释：
常见的角点：

相关名词：

两个像素块（大小相同）：$I_1(x)$,$I_2(x)$
窗口函数：${\Sigma }_{x,y}w(x,y)$
灰度变化：$E(u,v)$
角点响应函数：$R$
权重矩阵：$W$(通常为高斯滤波器$G_{\sigma }$)

相关计算：

1. 将窗口平移$[u,v]$产生灰度变化$E(u,v)$

   $$E(u,v)={\Sigma }_{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)$$

   $$E(u,v)\approx {\Sigma }_{x,y}w(x,y)[I_{x}u+I_{y}v{]}^2\cong \left[\begin{array}{c}uv\end{array}\right]M\left[\begin{array}{c}u\\ v\end{array}\right]$$
   其中，

   $$M={\Sigma }_{x,y}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$$

2. 点$x$的对称半正定矩阵:
   $$M_I=M_I(x)=\nabla I\nabla I^T=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]\left[\begin{array}{c}{I}_x{I}_y\end{array}\right]=\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$$
   $M_I$周围像素上的局部平均为：
   $$\overline{M_I}=W\ast M_I$$
   通过式 $\lambda E-M_I=0$，得到特征值 ${\lambda }_1=|\nabla I{|}^2$，${\lambda }_2=0$
   特征值会依赖于局部图像变化，具体表现为：图像梯度变化$⟹$$\overline{M_I}$特征值变化。

不同的特征值对应不同的像素情况：
a. ${\lambda }_1$、${\lambda }_2$都是很大的正数$⟹$$x$为角点
b. ${\lambda }_1$很大，${\lambda }_2=0$$⟹$区域存在一个边，$\overline{M_I}$变化不大
c. ${\lambda }_1={\lambda }_2=0$$⟹$区域为空

3. 角点响应
   $$R=det{M}_I-k(trace{M}_I{)}^2$$或
   $$R=\frac{det{M}_I}{(trace{M}_I{)}^2}$$
   其中，
   $$det{M}_I={\lambda }_1{\lambda }_2$$
   $$trace{M}_I={\lambda }_1+{\lambda }_2$$

R 只与M的特征值有关
角点：R 为大数值正数
边缘：R为大数值负数
平坦区：R为小数值

4. 角点进行挑选方法：
   a. 设置阈值。选取像素值高于阈值的所有图像点
   b. 角点之间间隔必须大于设定的最小距离

5. 在图像中寻找对应点：
   Harris角点检测器仅能够检测出图像中的兴趣点，但没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。
   Harris角点的描述子通常是由周围图像像素块的灰度值以及用于比较的归一化互相关矩阵构成的。
   两个像素块$I_1(x)$,$I_2(x)$的相关矩阵为
   $$c(I_1,I_2)={\Sigma }_{x}f(I_1(x),I_2(x)),f指的是相关方法$$
   互相关矩阵为
   $$f(I_1,I_2)=I_1{I}_2$$
   $$\therefore c(I_1,I_2)=I_1\cdot I_2$$
   $\cdot$ 指向量乘法，$c(I_1,I_2)$值越大，$I_1,I_2$相似度越高。
   归一化互相关矩阵为：
   $$ncc(I_1,I_2)=\frac{1}{n-1}{\Sigma }_x\frac{(I_1(x)-{\mu }_1)}{{\sigma }_1}\cdot \frac{(I_2(x)-{\mu }_2)}{{\sigma }_2}$$
   $n$指像素块中像素的数目，像素块即以该像素点为中心的周围矩形部分图像
   $\mu$指平均像素强度
   $\sigma$指像素块的标准差

2 实验部分

1. 角点检测
   a. 返回角点响应值$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不具有尺度不变性和旋转不变性，算法中的像素块的大小也会影响对应匹配的结果。