Python图像处理笔记——基于Harris角点的图像匹配（skimage）

---

一、前言

基于Harris角点实现两幅（不同角度拍摄）的图像匹配，基于RANSAC算法进行鲁棒性估计。算法的流程如下：

1. 读取图像，进行仿射变换；
2. 计算两幅图像中的角点（可以类比为其他特征）；
3. 考虑点周围的小空间，使用误差平方加权和计算点之间的对应点；（该度量方法只在两幅图像非常相似的情况下可用）
4. 前一步中得到的对应点为一组坐标值（两幅不同图像），可以用来估计它们之间的几何变换；
5. 采用随机抽样一致性算法进行鲁棒性估计。

二、使用步骤

1.引入库

import numpy as np
from skimage.io import imread
from skimage.color import rgb2gray
from skimage.util import img_as_float
from skimage.exposure import rescale_intensity
from skimage.transform import AffineTransform, warp
from skimage.feature import corner_harris, corner_peaks, corner_subpix，plot_matches
from skimage.measure import ransac
from matplotlib import pylab

2.基于skimage的仿射变换

	temple = rgb2gray(img_as_float(imread(r'../../images/building.jpg'))) # 这样写可以读取两级相对路径
    img_orig = np.zeros(list(temple.shape) + [3])
    img_orig[..., 0] = temple
    grad_row, grand_col = (np.mgrid[0:img_orig.shape[0],
                           0:img_orig.shape[1]] / float(img_orig.shape[0]))
    img_orig[..., 1] = grad_row
    img_orig[..., 2] = grand_col
    img_orig = rescale_intensity(img_orig)
    img_orig_gray = rgb2gray(img_orig)
    affine_trans = AffineTransform(scale=(0.8, 0.9), rotation=0.1, translation=(120, -20))
    image_warped = warp(img_orig, affine_trans.inverse, output_shape=img_orig.shape)
    image_warped_gray = rgb2gray(image_warped)

3. 亚像素角点计算

    # 计算图像中的感兴趣点--------------Harris角点
    # harris corners
    coordinates_original = corner_harris(img_orig_gray)
    # threshold for an optimal value, depends on the image
    coordinates_original[coordinates_original > 0.01 * coordinates_original.max()] = 1
    corner_coordinates_original = corner_peaks(coordinates_original, threshold_rel=0.0001, min_distance=5)

    # harris corners
    coordinates_warped = corner_harris(image_warped_gray)
    # threshold for an optimal value, depends on the image
    coordinates_warped[coordinates_warped > 0.01 * coordinates_warped.max()] = 1
    corner_coordinates_warped = corner_peaks(coordinates_warped, threshold_rel=0.0001, min_distance=5)

    # 确定子像素角点的位置
    coordinates_subpix_original = corner_subpix(img_orig_gray, corner_coordinates_original, window_size=9)
    coordinates_subpix_warped = corner_subpix(image_warped_gray, corner_coordinates_warped, window_size=9)

4. 基于中心像素距离的权重计算

# 高斯权重：
def gaussian_weight(window_ext, sigma=1):
    y, x = np.mgrid[-window_ext:window_ext + 1, -window_ext:window_ext + 1]
    g_w = np.zeros(y.shape, dtype=np.double)
    g_w[:] = np.exp(-0.5 * (x ** 2 / sigma ** 2 + y ** 2 / sigma ** 2))
    g_w /= 2 * np.pi * sigma * sigma
    return g_w   

5. 角点匹配

权重像素取决于到中心像素的距离，计算变形后图像中所有角点的误差平方加权和，用最小误差平方加权和的角点作为对应点

def match_corner(coord, window_ext=3):
	row, col = np.round(coord).astype(np.intp)
	window_original = img_orig[row - window_ext:row + window_ext + 1, col - window_ext:col + window_ext + 1, :]
	# 基于中心像素距离的权重计算
	weights = gaussian_weight(window_ext, 3)
	weights = np.dstack((weights, weights, weights))
	
	# 计算所有角点的误差平方加权和
	SSDs = []
	for coord_row, coord_col in corner_coordinates_warped:
	window_warped = image_warped[coord_row - window_ext:coord_row + window_ext + 1,
	coord_col - window_ext:coord_col + window_ext + 1, :]
	if window_original.shape == window_warped.shape:
	SSD = np.sum(weights * (window_original - window_warped) ** 2)  # 加权求和
	SSDs.append(SSD)
	
	min_idx = np.argmin(SSDs) if len(SSDs) > 0 else -1
	return coordinates_subpix_warped[min_idx] if min_idx >= 0 else [None]

6. 确定匹配对应点

RANSAC算法进行鲁棒性估计。先将点分为内点和外点，然后在忽略外点的情况下将模型拟合到内点上，寻找与仿射变换一致的匹配。

    src = []
    dst = []
    for coord in coordinates_subpix_original:
        coord1 = match_corner(coord) # 匹配点
        if any(coord1) and len(coord1) > 0 and not all(np.isnan(coord1)):
            src.append(coord)
            dst.append(coord1)
    src = np.array(src)
    dst = np.array(dst)

    # 使用全部坐标估计仿射变换模型
    model = AffineTransform()
    model.estimate(src, dst)

    # 鲁棒性估计基于RANSAC
    model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3, residual_threshold=2, max_trials=100)
    outliers = inliers == False

    # 比较真实和估计的变换参数
    print(affine_trans.scale, affine_trans.translation, affine_trans.rotation)
    print(model.scale, model.translation, model.rotation)
    print(model_robust.scale, model_robust.translation, model_robust.rotation)

7. 可视化图像匹配结果

	fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))
    plt.gray()

    inlier_idxs = np.nonzero(inliers)[0]
    plot_matches(ax[0], img_orig_gray, image_warped_gray, src, dst,
                 np.column_stack((inlier_idxs, inlier_idxs)), matches_color='b')
    ax[0].axis('off')
    ax[0].set_title('Correct correspondences', size=20)

    outlier_idxs = np.nonzero(outliers)[0]
    plot_matches(ax[1], img_orig_gray, image_warped_gray, src, dst,
                 np.column_stack((outlier_idxs, outlier_idxs)), matches_color='r')
    ax[1].axis('off')
    ax[1].set_title('Faulty correspondences', size=20)

    fig.tight_layout()

    plt.show()

---