通过函数拟合实现亚像素模板匹配

本文最后更新于:2023年4月15日 晚上

OpenCV 的模板匹配仅能实现像素级的模板匹配,如果想要在其基础上实现亚像素精度可以使用函数拟合的方法。

像素级模板匹配

像素级模板匹配可以用 OpenCV 的 matchTemplate 函数实现,算法运用快速傅里叶变换求图像相关度,以此为基础实现像素级模板匹配。

受到方法的限制,矩阵仅能按照整数值位置平移,因此无法直接实现亚像素精度。

上采样亚像素匹配

一种直接的调整思路是将图像进行差值上采样,在上采样后的图像上进行像素级模板匹配,得到结果后再按照采样倍率计算得到亚像素精度的匹配位置;

该方法在理论上是可行而且应该是有效的,但是这个方法比较大的问题是会大量增加运算量,例如想要得到0.1像素精度需要图像尺寸增大10倍,那么图像的面积就会增加100倍,性价比是很低的,所以不推荐这种方案。

《机器视觉算法与应用》中提到了可以用函数拟合的方法提升模板匹配精度,我尝试做了实现。

函数拟合亚像素匹配

拟合亚像素匹配的思想是在像素级模板匹配已经找到了像素级的最优匹配位置之后,拟合最优点周围的 Loss 值构成的三维流形函数,之后找到函数的最小值作为亚像素模板匹配的结果。

实现步骤

  1. 完成像素级模板匹配,并记录在匹配最优位置周围的 Loss 数据

    1
    2
    3
    4
    5
     [	[11350720.,  5760784.,  5126464.,  9760384., 17425104.]	,
    [ 7811536., 2045104., 1459056., 6331616., 14327520.] ,
    [ 7130720., 1393984., 902112., 5880512., 13967680.] ,
    [ 9504192., 3992352., 3592784., 8491984., 16391968.] ,
    [14032592., 8856816., 8525248., 13176016., 20679584.] ]
  2. 设计模型拟合上述数据,多项式当然是最常用的了,这里我用比较复杂的 4 次多项式,其实 2 次、3 次的结果也差不多,模型:

    $$ L(x,y)=p_1x^4+p_2x^3y+p_3x^2y^2+p_4xy^3+p_5y^4+p_6x^3+p_7x^2y+p_8xy^2+p_9y^3+p_{10}x^2+p_{11}xy+p_{12}y^2+p_{13}x+p_{14}y+p_{15} $$

    这里是一个二元多项式形式的模型,可以用伪逆矩阵求解最小二乘解

    以 ${\bf p} \in R^{15}$ 为未知数的模型,在最优点附近每一对整数 $(x,y)$ 都计算出了 $Loss$,那么此时就能写出形如
    $$
    Ap=b
    $$
    的等式关系,可以据此求出 $\bf p$ 的最小二乘解:
    $$
    p=A^+b
    $$

  3. 求解出模型后,还需要获取三维流形最小值所在的位置,可以使用最优化的方法直接限制取值范围求局部最优解,将该值作为亚像素模板匹配的结果。

示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import cv2
import mtutils as mt
import numpy as np
from scipy.optimize import minimize


def img_subpixel_match(img, golden_img):
H, W = golden_img.shape[:2]

# 像素级模板匹配,获取 Loss Map
res = cv2.matchTemplate(img, golden_img, cv2.TM_SQDIFF)
# 最优位置
_, _, (x_pos, y_pos), _ = cv2.minMaxLoc(res)

# 构造数据 A b
data_list = list()
score_list = list()
for dx in range(-4, 5, 1):
for dy in range(-4, 5, 1):
X = x_pos + dx
Y = y_pos + dy
# 模型
data = [X**4, Y**4, Y**3 * X, X**3 * Y,Y**2 * X**2, X**3, Y**3, Y**2 * X, X**2 * Y, X*Y, Y*Y, X*X, X, Y, 1]
data_list.append(data)
score_list.append(res[Y,X])

A = np.array(data_list)
B = np.array([score_list]).T

# 伪逆矩阵求解模型参数
w = np.dot(np.linalg.pinv(A), B)

# 绘图用数据
loss_list_list = list()
for dx in range(-7, 8, 1):
loss_list = list()
for dy in range(-7, 8, 1):
X = x_pos + dx
Y = y_pos + dy
data = [X**4, Y**4, Y**3 * X, X**3 * Y,Y**2 * X**2, X**3, Y**3, Y**2 * X, X**2 * Y, X*Y, Y*Y, X*X, X, Y, 1]
d = np.array([data]).T
loss = np.dot(w.T, d)
loss_list.append(loss)
loss_list_list.append(loss_list)
T = np.array(loss_list_list).squeeze()

temp_score = res[y_pos-7:y_pos+8, x_pos-7:x_pos+8]

L = temp_score.T
Xs, Ys = np.meshgrid(np.arange(15), np.arange(15))

# 绘图
plt = mt.plt
fig = plt.figure(figsize=(12, 6)) # 创建一个大小为12*6的绘图空间
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(Xs, Ys, L)
ax.plot_surface(Xs, Ys, T)
col1, col2 = cm.jet(np.array([0.15, 0.8]))
col1_patch = mpatches.Patch(color=col1, label='Origin Data')
col2_patch = mpatches.Patch(color=col2, label='Model Output')
plt.legend(handles=[col1_patch, col2_patch])
plt.show()
pass

# 最优化求解模型最小值位置
def Loss(x):
X, Y = x
data = [X**4, Y**4, Y**3 * X, X**3 * Y,Y**2 * X**2, X**3, Y**3, Y**2 * X, X**2 * Y, X*Y, Y*Y, X*X, X, Y, 1]
d = np.array([data]).T
loss = np.dot(w.T, d)
return loss[0][0]

x0 = np.asarray((x_pos, y_pos))
opt_res = minimize(Loss, x0, method='Powell')

sub_pix_x = opt_res.x[0]
sub_pix_y = opt_res.x[1]
return sub_pix_x, sub_pix_y

将 Halcon 得到的亚像素位置作为 GroundTruth 的话,该方法做亚像素模板匹配的平均误差在 0.05 个像素左右。

拟合效果

对于我手头的一批数据,展示不同阶数模型的三维流形示例与平均精度:

四阶模型

in pixel X Y
平均误差 0.02734 0.05427

三阶模型

in pixel X Y
平均误差 0.02314 0.05105

二阶模型

in pixel X Y
平均误差 0.02722 0.03232

参考资料



文章链接:
https://www.zywvvd.com/notes/study/image-processing/template-match/subpix-tm-fitting/subpix-tm-fitting/


通过函数拟合实现亚像素模板匹配
https://www.zywvvd.com/notes/study/image-processing/template-match/subpix-tm-fitting/subpix-tm-fitting/
作者
Yiwei Zhang
发布于
2023年4月8日
许可协议