阿齐兹的PyCV学习笔记——Harris角点检测器

Harris角点检测器

Harris角点检测算法：如果像素周围显示存在多于一个方向的边，就认为该点是兴趣点，该点就称为角点

对图像中的每一个像素，计算Harris矩阵，Harris矩阵的特征值有三种情况：

• 都是很大的正数，则该点是角点
• 一个很大，一个约等于0，则该区域存在一个边
• 都约等于0，该区域为空

在不需要计算特征值的情况下，引入指示函数区分重要的情况

定义角点响应函数，参数sigma定义了高斯滤波器的尺度大小，返回像素值为Harris响应函数的一幅图像。需要使用scipy.ndimage.filters模块中的高斯导数滤波器计算导数，因为在角点检测中需要抑制噪声

from scipy.ndimage import filters
def compute_harris_response(im,sigma=3):
    """在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数"""

    #计算导数
    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)

    #计算Harris矩阵分量
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    #计算特征值和迹
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    return Wdet / Wtr

定义函数get_harris_points()，选取像素值高于阈值的所有图像点，再加上额外的限制：角点之间的间隔必须大于设定的最小距离，返回以角点响应值递减的顺序排序的角点

def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    """从一幅Harris响应图像中返回角点，min_dist为分割角点和图像边界的最少像素数目"""

    #寻找高于阈值的候选角点
    corner_threshold = harrisim.max()*threshold
    harrisim_t = (harrisim>corner_threshold)*1

    #得到候选点的坐标
    coords = array(harrisim_t.nonzero()).T

    #以及他们的Harris响应值
    candidate_values = [harrisim[c[0],c[1]] for c in coords]

    #对候选点按照Harris响应值进行排序
    index = argsort(candidate_values)

    #将可行点的位置保存到数组中
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1

    #按照min_distance原则，选择最佳Harris点
    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

使用Matplotlib模块绘制角点

def plot_harris_points(image,filtered_coords):
    """绘制图像中检测到的角点"""

    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()

检测角点的测试代码：

im = array(Image.open('C:/Users/0AQZ0/Documents/exercisecode/Python/PyCV/Images/001.jpg').convert('L'))
harrisim = compute_harris_response(im)
filtered_coords = get_harris_points(harrisim,6)
plot_harris_points(im,filtered_coords)

效果图：

兴趣点描述子是分配给兴趣点的一个向量，描述该点附近图像的表观信息。描述子越好，寻找到的对应点越好。对应点是指描述相同物体和场景在不同图像上形成的像素点

Harris角点的描述子通常是周围图像像素块的灰度值，以及用于比较的归一化互相关矩阵构成的。图像的像素块由以该像素点为中心的周围矩阵部分图像构成

get_descriptors()方法可以将图像块像素值压平成一个向量，然后添加到描述子列表中，第一个参数是奇数大小长度的方形灰度图像块，图像块的中心是处理的像素点

def get_descriptors(image,filtered_coords,wid=5):
    """对每个返回的点，返回点周围2*wid+1个像素的值（假设选取点的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

match()方法使用归一化的互相关矩阵，将每个描述子匹配到另一个图像中的最优的候选点

def match(desc1,desc2,threshold=0.5):
    """对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配焦点"""

    n=len(desc1[0])

    #点对的距离
    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

match_twosided()方法可以从第二幅图像向第一幅图像匹配，过滤掉在两种方法中不都是最好的匹配

def match_twosided(desc1,desc2,threshold=0.5):
    """两边对称版本的match()"""

    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)

    ndx_12 = where(matches_12 >= 0)[0]

    #去除非对称的匹配
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1
    return matches_12

appendimages()方法、plot_matches()方法可以实现匹配点的可视化，绘制图像使用线段连接匹配的像素点

def appendimages(im1,im2):
    """返回将两幅图像并排拼接成的一幅新图像"""

    #选取具有最少行数的图像，然后填充足够的空行
    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)
    #如果这些情况都没有，那么它们的行数相同，不需要进行填充

    return concatenate((im1,im2),axis=1)

def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
    """显示一幅带有连接匹配之间连线的图片
       输入：im1，im2（数组图像），locs1，locs2（特征位置），matchscores（match（）的输出），
       show_below（如果图像应该显示在匹配的下方）"""

    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')

使用归一化互相关矩阵寻找对应点的实例：

im1 = array(Image.open('C:/Users/0AQZ0/Desktop/中控杯/视觉识别代码/images/01.jpg').convert('L'))
im2 = array(Image.open('C:/Users/0AQZ0/Desktop/中控杯/视觉识别代码/images/02.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 'start matching'
matches = match_twosided(d1,d2)

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

效果图：

缺点：算法的结果存在一些不正确匹配，因为图像像素块的互相关矩阵具有较弱的描述性，而且描述符不具有尺度不变性和旋转不变性，算法中像素块的大小也会影响匹配的结果