局部图像描述子------Harris角点检测器

一.Harris角点检测

1.实验原理

1.1什么是角点

角点就是极值点，即在某方面属性特别突出的点，是在某些属性上强度最大或者最小的孤立点，线段的终点。而对于图像而言，即是图像的角点，其是物体轮廓线的连接点。因此在角点时，窗口任意方向的移动都导致图像灰度的明显变化。

数学表达式

假设图像像素点(x,y)的灰度为I(x,y),以像素为中心的窗口沿x和y方向分别移动u和v的灰度强度变化的表达式为：

其中 E(u,v)是灰度变化，w(x,y) 是窗口函数，一般是高斯函数，所以可以把w(x,y)看做是高斯滤波器。I(x,y)是图像灰度, I(x+u,y+v)是平移后的图像灰度。在这里我们可以将 I(x+u,y+v)函数在（x,y）处泰勒展开：
将[ Ixu+Iyv ]展开后整理可以用矩阵表达为：
最后我们可以近似得到E(x,y)的表达式,将其化为二次型后得到：
其中M是一个2X2的矩阵，称为像素点的自相关矩阵，可以由图像的导数求得。M=窗口函数*偏导矩阵，表达式为：
因为u,v是局部微小的移动变量，所以我们对M进行讨论，M是一个2X2的矩阵，M的表达式中与点的位置（x,y）具体强相关性，记M得特征值为λ1,λ2。我们可以简单理解为该点的灰度值变化速度，那么a1和a2可以分别看做是x方向和y方向的灰度变化速率，就可以用a1,a2两者的大小关系来进行分类。
（1）当两个特征值λ1和λ2都偏小的时候，表示窗口沿任意方向移动都会使灰度变化很细微，该点处于图像的平坦区域。
（2）当λ1>>λ2或者λ1<<λ2时，说明该点向水平（垂直）方向移动时变化会很明显，而向垂直（水平）方向则变化不明显，该点处于图像的边缘区。
（3）当两个特征值λ1和λ2都很大的时候，表示窗口沿任意方向移动都会使灰度变化很明显，该点位置就是图像角点的位置。

然而在实际中，经常使用的是角点响应函数CRF这一概念，以此更加准确的计算所需角点,方法如下：
当|R|很小时，R<T , 认为该点处于图像的平坦区域。
当R<0时，R<T , 认为该点处于图像的边缘区。
当R>0时，R>T, 认为该点位置就是图像角点。

2.实验代码

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
im = cv2.imread(r'D:\software\pycharm\PycharmProjects\computer-version\two\images\IMG_20230412_093735.jpg', cv2.IMREAD_GRAYSCALE)

# 计算 Harris 响应函数
def compute_harris_response(im, sigma=3):
    # 计算 Sobel 导数
    dx = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=25)
    dy = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=25)

    # 计算卷积核
    ksize = int(6 * sigma + 1)
    w = cv2.getGaussianKernel(ksize, sigma)
    w = w * w.T

    # 计算 Harris 响应函数
    w_dx2 = cv2.filter2D(dx*dx, -1, w)
    w_dy2 = cv2.filter2D(dy*dy, -1, w)
    w_dxdy = cv2.filter2D(dx*dy, -1, w)

    det = w_dx2 * w_dy2 - w_dxdy**2
    trace = w_dx2 + w_dy2
    response = det - 0.05 * trace**2

    return response

# 不同阈值下的 Harris 角点检测
thresholds = [0.01, 0.05, 0.1] # 阈值列表
min_dist = 10 # 最小距离

# 非极大值抑制
def get_harris_points(harrisim, threshold=0.01, min_dist=10):
    # 找到高于阈值的像素
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

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

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

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

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

    # 按照最小距离原则选择最佳 Harris 点
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i, 0], coords[i, 1]] == np.array([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

for threshold in thresholds:
    harrisim = compute_harris_response(im)
    filtered_coords = get_harris_points(harrisim, threshold, min_dist)

    # 显示结果
    plt.figure()
    plt.gray()
    plt.imshow(im)
    plt.plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], 'r.')
    plt.axis('off')
    plt.title('Threshold = {}'.format(threshold))
    plt.show()

这段代码实现了 Harris 角点检测算法。具体来说，它首先读取一张灰度图像，然后计算 Harris 响应函数。接着，它使用非极大值抑制方法找到图像中的角点，并将它们保存在 filtered_coords 数组中。最后，它显示了原始图像，并在角点处绘制了红色的点。其中，compute_harris_response 函数计算了 Harris 响应函数，它使用了 Sobel 算子计算图像的梯度，然后使用高斯卷积核计算了每个像素的权重。最后，它使用这些权重计算了 Harris 响应函数。get_harris_points 函数实现了非极大值抑制方法，它首先找到高于阈值的像素，然后对这些像素进行排序，按照最小距离原则选择最佳的 Harris 点，并将它们保存在 filtered_coords 数组中。最后，代码使用 Matplotlib 库显示了原始图像，并在角点处绘制了红色的点。并且展示了不同阈值下的 Harris 角点检测。

3. 结果展示

分析：上面分别展示了阈值分别为：0.01，0.05，0.1的Harris角点检测结果。阈值越大，检测出的角点越少。

二.Harris特征匹配

1.实验原理

Harris角点检测器可以给出图像中检测到兴趣点，但是它并没有提供在图像间对兴趣点进行比较的方法，因此我们需要在每个角点添加描述子，以及对这些描述子进行比较。兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。

2.实验代码

# 导入必要的库
import cv2
import numpy as np

# 读取两张图像
img1 = cv2.imread(r'D:\software\pycharm\PycharmProjects\computer-version\two\images\IMG_20230413_105726.jpg')
img2 = cv2.imread(r'D:\software\pycharm\PycharmProjects\computer-version\two\images\IMG_20230413_105730.jpg')

# 转换为灰度图像
gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)

# 使用SIFT算法检测关键点和描述符
sift = cv2.xfeatures2d.SIFT_create()
kp1, des1 = sift.detectAndCompute(gray1, None)
kp2, des2 = sift.detectAndCompute(gray2, None)

# 使用FLANN匹配器进行特征点匹配
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des1, des2, k=2)

# 选择最佳匹配
good_matches = []
for m, n in matches:
    if m.distance < 0.7 * n.distance:
        good_matches.append(m)

# 绘制匹配结果
img_matches = cv2.drawMatches(img1, kp1, img2, kp2, good_matches, None, flags=2)

# 显示匹配结果
cv2.imshow('Matches', img_matches)
cv2.waitKey(0)
cv2.destroyAllWindows()

3.结果展示

分析：该算法的结果存在很多不正确匹配。这是因为，与现代的一些方法（下面将会提到）相比，图像像素块的互相关矩阵具有较弱的描述性。实际运用中，我们通常使用更稳健的方法来处理这些对应匹配。这些描述符还有一个问题，它们不具有尺度不变性和旋转不变性，而算法中像素块的大小也会影响对应匹配的结果。