Harris角点检测

Harris角点检测算法（也称Harris检测器），主要思想：若该像素周围存在多个方向的边，则该像素点即为角点

1. 原理分析

什么是角点

局部窗口沿各方向移动，均产生明显变化的点
图像局部曲率突变的点

Harris角点检测基本思想

平坦区域
任意方向移动均无灰度变化
边缘
沿着边缘方向移动无灰度变化
角点
沿着任意方向均有明显灰度变化

数学表达

数学定义

将图像窗口平移[u,v]产生灰度变化E(u,v)

对公式进行展开化简:
延伸知识点：一元函数泰勒展开
延伸知识点：二元函数泰勒展开

2. 示例

Harris 角点检测步骤

读入图像转灰度图

# 读入图像ort compute_harris_response
im = array(Image.open('Image/revolve-1.jpg').convert('L'))

计算响应函数

harrisim = compute_harris_response(im, 3)


# 在灰度图像中，对每个像素计算Harris角点响应函数
def compute_harris_response(im, sigma=3):
    # 计算导数
    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)

    # 矩阵分量
    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

# Harris响应函数
harrisim1 = 255 - harrisim

基于响应值选择角点


# 从响应图中返回角点
def get_harris_points(harrisim, min_dist=10, threshold=0.5):
    # 挑选高于阈值
    corner_threshold = threshold
    # coner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    # 得到坐标
    coords = array(harrisim_t.nonzero()).T
    # 响应值
    candidate_values = [harrisim[c[0], c[1]]for c in coords]

    # 按响应值排序
    index = argsort(candidate_values)[::-1]

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

    # 选择最佳角点
    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

在原图中画出检测的角点

# 画出Harris响应图
subplot(221)
imshow(im)
axis('off')
axis('equal')

threshold = [0.1, 0.2, 0.3]
for i, thres in enumerate(threshold):
    filtered_coords = get_harris_points(harrisim, 50, thres)
    subplot(2, 2, i+2)
    title(str(threshold[i]))
    imshow(im)
    plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
    axis('off')

show()

def plot_harris_points(image, filtered_coords):
    figure(dpi=120)
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],
         [p[0] for p in filtered_coords], '*')
    axis('off')
    show()

结果显示
在这里插入图片描述

SIFT特征

1.SIFT算法实现特征匹配主要有三个流程

提取关键点；
对关键点附加详细的信息（局部特征），即描述符
通过特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，建立景物间的对应关系

2. 哪些点是SIFT中要查找的关键点（特征点）？

十分突出的点,不会因光照、尺度、旋转等因素的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点。

3. 尺度空间

什么是尺度空间

主要思想是通过对原始图像进行尺度变换，获得图像多尺度下的空间如下表示。从而实
现边缘、角点检测和不同分辨率上的特征提取，以满足特征点的尺度不变性。
在这里插入图片描述

4. 特征匹配算法

提取图像特征点
提取特征点周围图像块，构造特征描述符
通过特征描述符对比，实现特征匹配


imname = 'Image/revolve-1.jpg'
im = array(Image.open(imname).convert('L'))
# 得到SIFT特征
process_image(imname, 'revolve.sift')
# 取出坐标及特征描述
l1, d1 = read_features_from_file('revolve.sift')
figure(dpi=180)
gray()
subplot(221)
plt.axis('off')
imshow(Image.open(imname))
title('src')
subplot(222)
# 可视化特征点
plot_features(im, l1, circle=False)
title('SIFT')
subplot(223)
plot_features(im, l1, circle=True)
title('scale')
# 检测harris角点，对比
harrisim = compute_harris_response(im)
subplot(224)
filtered_coords = get_harris_points(harrisim, 60, 0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
title('Harris')
show()



# 检测兴趣点
def process_image(imagename, resultname, params = "--edge-thresh 10 --peak-thresh 5"):

    if imagename[-3:] != 'pgm':
        im = Image.open(imagename).convert('L')
        im.save('tmp.pgm')
        imagename = 'tmp.pgm'
    cmmd = str(r'E:\computerstudy\python\SIFT\ch02\sift.exe '+imagename+" --output="+resultname +
               " "+params)
    os.system(cmmd)
    print('processed', imagename, 'to', resultname)


def read_features_from_file(filename):
    f = loadtxt(filename)
    return f[:, :4], f[:, 4:]


def plot_features(im, locs, circle=False):
    def draw_circle(c, r):
        t = arange(0, 1.01, .01)*2*pi
        x = r*cos(t) + c[0]
        y = r*sin(t) + c[1]
        plot(x, y, 'b', linewidth=2)
    imshow(im)
    if circle:
        for p in locs:
            draw_circle(p[:2], p[2])
    else:
        plot(locs[:, 0], locs[:, 1], 'ob')
    axis('off')
    return