Python计算机视觉第2章-局部图像描述子

2.1 Harris角点检测器

Harris角点检测算法（也称Harris & Stephens角点检测器）是一个极为简单的角点检测算法。该算法的主要思想是，如果像素周围显示存在多于一个方向的边，我们认为该点为兴趣点。该点就称为角点。

Harris角点检测算法

Harris角点检测算法基于图像灰度值的变化来识别角点。角点通常出现在图像中灰度值变化剧烈的区域，特别是梯度在两个方向上都有较大变化的地方。

具体来说，算法的步骤如下：

计算图像梯度：
- 计算图像在x方向和y方向上的梯度，通常使用Sobel算子或其他类似的滤波器来计算。
构造矩阵 $M$ ：
- 对于图像中的每个像素，计算自相关矩阵 $M$ ：
  
  $\begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix}$
  
  其中， $I_x$ 和 $I_y$ 是图像在 x 方向和 y 方向的梯度，矩阵中的各项可以通过对梯度值进行加权求和得到。
响应函数 $R$ ：
- 通过矩阵 $M$ 的特征值分析，计算响应函数 $R$ ：
  
  $\text{det}(M) - k \cdot (\text{trace}(M))^2$
  
  其中， $\text{det}(M)$ 是 $M$ 的行列式， $\text{trace}(M)$ 是 $M$ 的迹， $k$ 是一个经验常数，通常取值为 0.04 到 0.06 之间。
角点检测：
- 对于每个像素，计算对应的 $R$ 值。如果 $R$ 值大于某个阈值并且是局部最大值，则该像素被认为是一个角点。

以下是算法的示例代码：

from matplotlib import pyplot as plt
from PIL import Image
import numpy as np
from scipy.ndimage import gaussian_filter


def compute_harris_response(img, sigma=3):
    # 计算图像的x和y方向的梯度
    imgx = np.zeros(img.shape)
    gaussian_filter(img, (sigma, sigma), (0, 1), imgx)
    imgy = np.zeros(img.shape)
    gaussian_filter(img, (sigma, sigma), (1, 0), imgy)

    # 计算Harris矩阵的分量
    Wxx = gaussian_filter(imgx * imgx, sigma)
    Wxy = gaussian_filter(imgx * imgy, sigma)
    Wyy = gaussian_filter(imgy * imgy, sigma)

    # 计算Harris响应，并处理可能的除零情况
    Wdet = Wxx * Wyy - Wxy * Wxy
    Wtran = Wxx + Wyy
    Wtran[Wtran == 0] = 1e-10  # 避免除零错误
    return Wdet / Wtran


def get_harris_points(harrisimg, min_dist=10, threshold=0.1):
    # 找到角点
    corners = harrisimg.max() * threshold
    harrisimg_t = (harrisimg > corners) * 1
    coords = np.array(harrisimg_t.nonzero()).T
    candidate_values = [harrisimg[c[0], c[1]] for c in coords]
    index = np.argsort(candidate_values)

    # 将坐标保存在允许的位置
    allow_locations = np.zeros(harrisimg.shape)
    allow_locations[min_dist:-min_dist, min_dist:-min_dist] = 1
    filtered_coords = []
    for idx in index:
        if allow_locations[coords[idx][0], coords[idx][1]] == 1:
            filtered_coords.append(coords[idx])
            allow_locations[coords[idx][0] - min_dist:coords[idx][0] + min_dist,
            coords[idx][1] - min_dist:coords[idx][1] + min_dist] = 0
    return filtered_coords


# 加载图像并进行Harris角点检测
img = Image.open(r'img.png').convert('L')
img = np.array(img)
harrisimg = compute_harris_response(img)
filtered_coords1 = get_harris_points(harrisimg, 6, threshold=0.1)
filtered_coords2 = get_harris_points(harrisimg, 6, threshold=0.05)

# 显示结果并保存图像
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(141), plt.imshow(img, cmap='gray'), plt.title('原始图像'), plt.axis('off')
plt.subplot(142), plt.imshow(harrisimg, cmap='gray'), plt.title('角点图像'), plt.axis('off')
plt.subplot(143), plt.imshow(img, cmap='gray'), plt.plot([p[1] for p in filtered_coords1],
                                                         [p[0] for p in filtered_coords1], '.'), plt.title(
    '阈值为0.1'), plt.axis('off')
plt.subplot(144), plt.imshow(img, cmap='gray'), plt.plot([p[1] for p in filtered_coords2],
                                                         [p[0] for p in filtered_coords2], '.'), plt.title(
    '阈值为0.05'), plt.axis('off')

# 保存图像为文件
plt.savefig('harris_corner_detection.png', dpi=300, bbox_inches='tight')

# 显示图像
plt.show()

在这里插入图片描述

实验图1 Harris角点检测算法处理结果

在图像间寻找对应点

Harris 角点检测器仅仅能够检测出图像中的兴趣点，但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点上加入描述子信息，并给出一个比较这些描述子的方法。
兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。
Harris 角点的描述子通常是由周围图像像素块的灰度值，以及用于比较的归一化互相关矩阵构成的。图像的像素块由以该像素点为中心的周围矩形部分图像构成。

通常，两个（相同大小）像素块 $I_1(x)$ 和 $I_2(x)$ 的相关矩阵定义为：
$c(\boldsymbol{I}_1,\boldsymbol{I}_2)=\sum_{x}f(\boldsymbol{I}_1(\mathbf{x}),\boldsymbol{I}_2(\mathbf{x}))$

其中，函数 $f$ 随着相关方法的变化而变化。上式取像素块中所有像素位置 $x$ 的和。对于互相关矩阵，函数 $f(I_1 ,I_2 )=I_1 I_2$ , 因此， $c(I_1,I_2 )=I_1 · I_2$ , 其中· 表示向量乘法（按照行或者列堆积的像素）。 $c(I_1,I_2)$ 的值越大，像素块 $I_1$ 和 $I_2$ 的相似度越高。

归一化的互相关矩阵是互相关矩阵的一种变形，可以定义为：
$ncc(I_1,I_2)=\frac{1}{n-1}\sum_{x}\frac{(I_1(\mathbf{x})-\mu_1)}{\sigma_1}\cdot\frac{(I_2(\mathbf{x})-\mu_2)}{\sigma_2}$

其中， $n$ 为像素块中像素的数目， $μ_1$ 和 $μ_2$ 表示每个像素块中的平均像素值强度， $σ_1$ 和 $σ_2$ 分别表示每个像素块中的标准差。通过减去均值和除以标准差，该方法对图像亮度变化具有稳健性。

2.2　SIFT（尺度不变特征变换）

David Lowe 在文献[17]中提出的SIFT（Scale-Invariant Feature Transform，尺度不变特征变换）是过去十年中最成功的图像局部描述子之一。SIFT特征后来在文献中得到精炼并详述，经受住了时间的考验。SIFT特征包括兴趣点检测器和描述子。SIFT描述子具有非常强的稳健性，这在很大程度上也是SIFT特征能够成功和流行的主要原因。自从SIFT特征的出现，许多其他本质上使用相同描述子的方法也相继出现。现在，SIFT描述符经常和许多不同的兴趣点检测器相结合使用（有些情况下是区域检测器），有时甚至在整幅图像上密集地使用。SIFT特征对于尺度、旋转和亮度都具有不变性，因此，它可以用于三维视角和噪声的可靠匹配。

2.2.1　兴趣点

SIFT 特征使用高斯差分函数来定位兴趣点：
$D(\mathbf{x},\sigma)=[G_{\kappa\sigma}(\mathbf{x})-G_{\sigma}(\mathbf{x})]*I(\mathbf{x})=[G_{\kappa\sigma}-G_{\sigma}]*I=I_{\kappa\sigma}-I_{\sigma}$

其中， $G_σ$ 是上一章中介绍的二维高斯核， $I_σ$ 是使用 $G_σ$ 模糊的灰度图像， $κ$ 是决定相差尺度的常数。兴趣点是在图像位置和尺度变化下 $D (x, σ)$ 的最大值和最小值点。这些候选位置点通过滤波去除不稳定点。基于一些准则，比如认为低对比度和位于边上的点不是兴趣点，我们可以去除一些候选兴趣点。

2.2.2　描述子

上面讨论的兴趣点（关键点）位置描述子给出了兴趣点的位置和尺度信息。为了实现旋转不变性，基于每个点周围图像梯度的方向和大小，SIFT描述子又引入了参考方向。SIFT描述子使用主方向描述参考方向。主方向使用方向直方图（以大小为权重）来度量。

2.2.3　检测兴趣点

以下是实验代码：

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

# 加载图像
img = cv2.imread('img.png', cv2.IMREAD_GRAYSCALE)
img_color = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB)  # 转换为 RGB 图像

# 计算 Harris 角点
dst = cv2.cornerHarris(img, 2, 3, 0.04)

# 结果膨胀以标记角点
dst = cv2.dilate(dst, None)

# 设置阈值以标记角点
img_color[dst > 0.01 * dst.max()] = [0, 0, 255]  # 使用红色标记角点

# 显示结果
plt.imshow(img_color)
plt.title('Harris Corners')
plt.axis('off')

# 保存结果图像
plt.savefig('harris_corners.png', dpi=300, bbox_inches='tight')

plt.show()

在这里插入图片描述

实验图2 处理结果

2.2.4　匹配描述子

对于将一幅图像中的特征匹配到另一幅图像的特征，一种稳健的准则（同样是由Lowe 提出的）是使用这两个特征距离和两个最匹配特征距离的比率。相比于图像中的其他特征，该准则保证能够找到足够相似的唯一特征。使用该方法可以使错误的匹配数降低。
以下为实验代码：

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

# 加载图像
img1 = cv2.imread('img.png', cv2.IMREAD_GRAYSCALE)  # 查询图像
img2 = cv2.imread('img.png', cv2.IMREAD_GRAYSCALE)  # 训练图像

# 检查图像是否成功加载
if img1 is None:
    raise ValueError("无法读取图像文件 'img1.png'。请检查文件路径和文件名。")
if img2 is None:
    raise ValueError("无法读取图像文件 'img2.png'。请检查文件路径和文件名。")

# 创建 SIFT 检测器
sift = cv2.SIFT_create()

# 检测关键点和计算描述符
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)

# 创建 BFMatcher 对象
bf = cv2.BFMatcher(cv2.NORM_L2, crossCheck=True)

# 进行描述符匹配
matches = bf.match(des1, des2)

# 按照距离排序
matches = sorted(matches, key=lambda x: x.distance)

# 绘制匹配结果
img_matches = cv2.drawMatches(img1, kp1, img2, kp2, matches[:50], None, flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)

# 显示匹配结果
plt.imshow(img_matches)
plt.title('Feature Matches')
plt.axis('off')

# 保存匹配结果图像
plt.savefig('feature_matches.png', dpi=300, bbox_inches='tight')

plt.show()

在这里插入图片描述

实验图3 处理结果

2.3　匹配地理标记图像

我们首先通过图像间是否具有匹配的局部描述子来定义图像间的连接，然后可视化这些连接情况。为了完成可视化，我们可以在图中显示这些图像，图的边代表连接。

2.3.1　使用局部描述子匹配

在这种情况下，我们将使用前面部分讲述的SIFT特征描述子。我们假设已经对这些图像使用SIFT特征提取代码进行了处理，并且将特征保存在和图像同名（但文件名后缀是.sift，而不是.jpg）的文件中。假设imlist和featlist列表中包含这些文件名。我们可以对所有组合图像对进行逐个匹配，以下为示例关键代码：

import sift
import numpy as np

# 图像列表和特征文件列表
imlist = [...]  # 替换为实际图像文件列表
featlist = [...]  # 替换为实际特征文件列表

# 初始化
nbr_images = len(imlist)
matchscores = np.zeros((nbr_images, nbr_images))

# 计算每对图像的匹配得分
for i in range(nbr_images):
    for j in range(i, nbr_images):  # 仅计算上三角
        print('Comparing', imlist[i], imlist[j])
        
        # 读取特征
        l1, d1 = sift.read_features_from_file(featlist[i])
        l2, d2 = sift.read_features_from_file(featlist[j])
        
        # 匹配描述符
        matches = sift.match_twosided(d1, d2)
        nbr_matches = np.sum(matches > 0)
        
        print('Number of matches =', nbr_matches)
        
        # 填充匹配得分矩阵
        matchscores[i, j] = nbr_matches

# 复制值以填充对称矩阵
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):  # 不需要复制对角线
        matchscores[j, i] = matchscores[i, j]

2.3.3　可视化连接的图像

我们首先通过图像间是否具有匹配的局部描述子来定义图像间的连接，然后可视化这些连接情况。为了完成可视化，我们可以在图中显示这些图像，图的边代表连接。以下为示例关键代码：

import pydot

# 创建一个无向图
g = pydot.Dot(graph_type='graph')

# 添加节点0，设置字体颜色为透明
g.add_node(pydot.Node(str(0), fontcolor='transparent'))

# 添加其他节点及边
for i in range(5):
    # 添加节点 i+1
    g.add_node(pydot.Node(str(i + 1)))
    # 添加边连接节点 0 和 i+1
    g.add_edge(pydot.Edge(str(0), str(i + 1)))
    
    for j in range(5):
        # 添加节点 j+1-i+1
        g.add_node(pydot.Node(str(j + 1) + '-' + str(i + 1)))
        # 添加边连接节点 j+1-i+1 和 j+1
        g.add_edge(pydot.Edge(str(j + 1) + '-' + str(i + 1), str(j + 1)))

# 将图写入文件
g.write_png('graph.jpg', prog='neato')