亚像素角点检测

本文最后更新于：2024年1月14日晚上

cv.goodFeaturesToTrack() 提取到的角点只能达到像素级别，在很多情况下并不能满足实际的需求，这时，我们则需要使用 cv.cornerSubPix() 对检测到的角点作进一步的优化计算，可使角点的精度达到亚像素级别。

检测

前面已经提及 goodFeaturesToTrack() 提取到的角点只能达到像素级别，获取的角点坐标是整数，但是通常情况下，角点的真实位置并不一定在整数像素位置，因此为了获取更为精确的角点位置坐标，需要角点坐标达到亚像素（subPixel）精度。这时，我们则需要使用cv::cornerSubPix()对检测到的角点作进一步的优化计算，可使角点的精度达到亚像素级别。

原理解析

在亚像素级精度的角点检测算法中，一种方法是从亚像素角点到周围像素点的矢量应垂直于图像的灰度梯度这个观察事实得到的，通过最小化误差函数的迭代方法来获得亚像素级精度的坐标值。

对于一个角点 $ q $ ，考虑在 $ q $ 附近的窗口 win 内的任意点 $ p $ 。 $ p $ 的图像梯度与向量 $ q-p $ 的点积总是为零。

原因是如果 $ p $ 在平坦区域，则 $ p $ 的梯度为零，如果 $ p $ 在边缘处， $ p $ 的梯度总是与 $ q-p $ 垂直。

据此可以列出公式：

$$ \begin{array}{c} I_{x}(p)\left(x_{q}-x_{p}\right)+I_{y}(p)\left(y_{q}-y_{p}\right)=0 \\ I_{x}(p) x_{q}+I_{y}(p) y_{q}=I_{x}(p) x_{p}+I_{y}(p) y_{p} \end{array} $$

联立窗口 win内的所有点 $ p_{i} $ 得到方程组，写成矩阵形式如下:

$$ \left(\begin{array}{cc}I_{x}\left(p_{1}\right) & I_{y}\left(p_{1}\right) \\ \vdots & \vdots \\ I_{x}\left(p_{n}\right) & I_{y}\left(p_{n}\right)\end{array}\right)\left(\begin{array}{c}x_{q} \\ y_{q}\end{array}\right)=\left(\begin{array}{c}I_{x}\left(p_{1}\right) x_{p_{1}}+I_{y}\left(p_{1}\right) y_{p_{1}} \\ \vdots \\ I_{x}\left(p_{n}\right) x_{p_{n}}+I_{y}\left(p_{n}\right) y_{p_{n}}\end{array}\right) $$

这是一个超约束(over-constrained)线性方程组，可以用最小二乘法求最佳解。

$$ A=\left(\begin{array}{cc}I_{x}\left(p_{1}\right) & I_{y}\left(p_{1}\right) \\ \vdots & \vdots \\ I_{x}\left(p_{n}\right) & I_{y}\left(p_{n}\right)\end{array}\right), b=\left(\begin{array}{c}I_{x}\left(p_{1}\right) x_{p_{1}}+I_{y}\left(p_{1}\right) y_{p_{1}} \\ \vdots \\ I_{x}\left(p_{n}\right) x_{p_{n}}+I_{y}\left(p_{n}\right) y_{p_{n}}\end{array}\right) $$ $$ \begin{array}{c} A^{T} A\left(\begin{array}{l}x_{q} \\ y_{q}\end{array}\right)=A^{T} b \\ A^{T} A=\left(\begin{array}{cc}\sum I_{x}^{2} & \sum I_{x} I_{y} \\ \sum I_{x} I_{y} & \sum I_{y}^{2}\end{array}\right) \\ A^{T} b=\left(\begin{array}{l}\sum I_{x}^{2} x+\sum I_{x} I_{y} y \\ \sum I_{x} I_{y} x+\sum I_{y}^{2} y\end{array}\right) \\ \left(\begin{array}{l}x_{q} \\ y_{q}\end{array}\right)=\left(A^{T} A\right)^{-1} A^{T} b \end{array} $$

或直接沿用伪逆求解最小二乘解的结论：

$$ \left(\begin{array}{l}x_{q} \\ y_{q}\end{array}\right)= A^{+} b $$

然后以新的 $ \left(x_{q}, y_{q}\right) $ 为初始角点，重新执行以上优化过程，反复迭代，直至 $ \left(x_{q}, y_{q}\right) $ 收敛。

OpenCV 函数

函数定义：

void cv::cornerSubPix(
	cv::InputArray image, // 输入图像
	cv::InputOutputArray corners, // 角点（既作为输入也作为输出）
	cv::Size winSize, // 区域大小为 NXN; N=(winSize*2+1)
	cv::Size zeroZone, // 类似于winSize，但是总具有较小的范围，Size(-1,-1)表示忽略
	cv::TermCriteria criteria // 停止优化的标准
);

参数详解：

参数	含义
image	输入图像，和 `cv::goodFeaturesToTrack()` 中的输入图像是同一个图像。
corners	检测到的角点，即是输入也是输出。
winSize	计算亚像素角点时考虑的区域的大小，大小为NXN; N=(winSize*2+1)。
zeroZone	作用类似于winSize，但是总是具有较小的范围，通常忽略（即Size(-1, -1)）。
criteria	表示计算亚像素时停止迭代的标准，可选的值有`cv::TermCriteria::MAX_ITER` 、`cv::TermCriteria::EPS`（可以是两者其一，或两者均选），前者表示迭代次数达到了最大次数时停止，后者表示角点位置变化的最小值已经达到最小时停止迭代。二者均使用 `cv::TermCriteria()`构造函数进行指定。

cornerSubPix 的实现中：

随着角点位置的细化，每次迭代都要重新计算窗口的像素值。由于中心坐标并非整数，因此整个窗口的像素坐标也不是整数，需要用插值算法来计算每个点的像素值
使用加权最小二乘法优化结果，用高斯核让算法给离中心近的点更高的权重

Python 实现

示例图像：

示例代码：

import numpy as np
import cv2 as cv
import mtutils as mt

# termination criteria
criteria = (cv.TERM_CRITERIA_EPS + cv.TERM_CRITERIA_MAX_ITER, 30, 0.001)

W_num = 7
H_num = 7

img = mt.cv_rgb_imread('undistort.png', 1)

ret, corners = cv.findChessboardCorners(img, (W_num,H_num), None)
corners2 = cv.cornerSubPix(img,corners, (11,11), (-1,-1), criteria)
cv.drawChessboardCorners(img, (W_num, H_num), corners2, ret)
print(np.squeeze(corners2).tolist())
mt.PIS(img)
pass

输出亚像素角点结果：

[[37.20782470703125, 28.19313621520996], [61.632537841796875, 25.13271713256836], [88.52682495117188, 22.974666595458984], [116.57200622558594, 22.350011825561523], [144.58685302734375, 22.965282440185547], [171.49815368652344, 25.12959861755371], [196.05784606933594, 28.172744750976562], [33.98270034790039, 53.4093017578125], [59.12421417236328, 50.916385650634766], [87.07998657226562, 49.2585563659668], [116.57009887695312, 48.41100311279297], [146.17893981933594, 49.26211166381836], [174.11178588867188, 50.913177490234375], [199.19203186035156, 53.42321014404297], [31.714290618896484, 81.21449279785156], [57.28245162963867, 79.64262390136719], [85.70196533203125, 78.4729995727539], [116.56983947753906, 78.20845031738281], [147.414794921875, 78.47453308105469], [175.7760772705078, 79.6445083618164], [201.4080810546875, 81.2107925415039], [31.063425064086914, 110.45411682128906], [56.46329116821289, 110.45072174072266], [85.46632385253906, 110.4486312866211], [116.5677490234375, 110.44801330566406], [147.5690155029297, 110.44960021972656], [176.6145477294922, 110.4505844116211], [202.2876739501953, 110.44986724853516], [31.70306396484375, 139.6763916015625], [57.28034973144531, 141.36480712890625], [85.69633483886719, 142.5038299560547], [116.56790924072266, 142.6646270751953], [147.42013549804688, 142.50144958496094], [175.77560424804688, 141.36170959472656], [201.41563415527344, 139.68089294433594], [33.96387481689453, 167.4934844970703], [59.08919906616211, 169.99136352539062], [87.06123352050781, 171.72630310058594], [116.5658187866211, 172.5388946533203], [146.1722412109375, 171.72296142578125], [174.13735961914062, 169.99850463867188], [199.19720458984375, 167.4800262451172], [37.18704605102539, 192.6940460205078], [61.63749313354492, 195.68472290039062], [88.49982452392578, 197.91934204101562], [116.56771850585938, 198.51622009277344], [144.60362243652344, 197.91807556152344], [171.47967529296875, 195.6854705810547], [196.0402069091797, 192.72105407714844]]