9点拟合梯度边缘亚像素方法

本文方法取自硕士毕业论文《基于几何特征的快速模板匹配算法》，豆丁网链接：基于几何特征快速模板匹配算法 - 豆丁网 (docin.com)。仅供参考，由于拟合二次曲面耗时较大，工业应用一般不使用这种方法，而用拟合平面抛物线代替。

本文思路是先用整数边缘点${ \left ( x_{0},y_{0} \right ) }$周围9个点和对应的梯度模长${ z }$拟合一个二次曲面，然后计算这个二次曲面在整数边缘点梯度方向上的曲线的顶点作为拟合的亚像素点。整数边缘点${ \left ( x_{0},y_{0} \right ) }$可以使用Canny算子得到。现在我们先拟合二次曲面，定义待求二次曲面函数如下：

$${ z=ax^{2}+by^{2}+cxy+dx+ey+f }$$

将整数坐标${ \left ( x_{0},y_{0} \right ) }$周围9个点依次代入上式。为了计算方便将这9个点都减去${ \left ( x_{0},y_{0} \right ) }$平移到原点处，对应上图Z(0,0)位置。计算完成后再平移回来，不影响效果。可得如下方程：

$${ \begin{pmatrix} 1 & 1 & 1 & -1 & -1 & 1 \\ 0 & 1 & 0 & 0 & -1 & 1 \\ 1 & 1 & -1 & 1 & -1 & 1 \\ 1 & 0 & 0 & -1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 & 0 & 1 \\ 1 & 1 & -1 & -1 & 1 & 1 \\ 0 & 1 & 0 & 0 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} a\\ b\\ c\\ d\\ e\\ f \end{pmatrix}=\begin{pmatrix} z\left ( -1,-1 \right ) \\ z\left ( 0,-1 \right ) \\ z\left ( 1,-1 \right ) \\ z\left ( -1,0 \right ) \\ z\left ( 0,0 \right ) \\ z\left ( 1,0 \right ) \\ z\left ( -1,1 \right ) \\ z\left ( 0,1 \right ) \\ z\left ( 1,1 \right ) \end{pmatrix} }$$

可以看出方程左侧系数全是常量，因此可以提前算好以节省时间。利用求解公式${ \mathbf{K}=\left ( \mathbf{X^{T}}\mathbf{X} \right )^{-1} \mathbf{X^{T}}\mathbf{Y} }$可得：

$${ \begin{pmatrix} a\\ b\\ c\\ d\\ e\\ f \end{pmatrix}=\begin{pmatrix} 0.166667 & -0.333333 & 0.166667 & 0.166667 & -0.333333 & 0.166667 & 0.166667 & -0.333333 & 0.166667 \\ 0.166667 & 0.166667 & 0.166667 & -0.333333 & -0.333333 & -0.333333 & 0.166667 & 0.166667 & 0.166667 \\ 0.25 & 0 & -0.25 & 0 & 0 & 0 & -0.25 & 0 & 0.25 \\ -0.166667 & 0 & 0.166667 & -0.166667 & 0 & 0.166667 & -0.166667 & 0 & 0.166667 \\ -0.166667 & -0.166667 & -0.166667 & 0 & 0 & 0 & 0.166667 & 0.166667 & 0.166667 \\ -0.111111 & 0.222222 & -0.111111 & 0.222222 & 0.555556 & 0.222222 & -0.111111 & 0.222222 & -0.111111 \end{pmatrix} \begin{pmatrix} z\left ( -1,-1 \right ) \\ z\left ( 0,-1 \right ) \\ z\left ( 1,-1 \right ) \\ z\left ( -1,0 \right ) \\ z\left ( 0,0 \right ) \\ z\left ( 1,0 \right ) \\ z\left ( -1,1 \right ) \\ z\left ( 0,1 \right ) \\ z\left ( 1,1 \right ) \end{pmatrix} }$$

然后就是计算整数坐标${ \left ( x_{0},y_{0} \right ) }$的梯度。可以用Sobel算子的梯度${ \left ( G_{x},G_{y} \right ) }$，也可以用上述二次曲面函数在此处的偏导数${ \left ( \frac{\partial z}{\partial x}, \frac{\partial z}{\partial y} \right ) }$作为梯度（由于我们平移了坐标所以实际上就是${ z\left ( x,y \right ) }$在原点处的梯度）。这两种梯度经过对比效果差不多。最终就是求下式所表示的三维曲线的顶点。式中k是梯度直线的斜率${ \frac{G_{y}}{G_{x}} }$，当然在三维空间中它代表平面：

$${ \left\{\begin{matrix} z=ax^{2}+by^{2}+cxy+dx+ey+f \\ y=kx,需要k存在 \end{matrix}\right. }$$

消去y并对x求导数，使其为0，有：

$${ \frac{\mathrm{d} z}{\mathrm{d} x} = 2ax+2bk\left ( kx \right )+2c\left ( kx \right )+d+ek=0 }$$

解之可得亚像素量，这里使用${ {x}',{y}' }$表示亚像素结果：

$${ \left\{\begin{matrix} {x}'=-\frac{ek+d}{2\left ( bk^{2}+ck+a \right )}\\ {y}'=k{x}' \end{matrix}\right. }$$

当k不存在时：

$${ \left\{\begin{matrix} {x}'=\lim_{k \rightarrow \infty }-\frac{ek+d}{2\left ( bk^{2}+ck+a \right )}=0 \\ {y}'=\lim_{k \rightarrow \infty }k{x}'=-\frac{e}{2b} \end{matrix}\right. }$$

最终亚像素边缘为：

$${ \left\{\begin{matrix} x_{sub}=x_{0}+{x}' \\ y_{sub}=y_{0}+{y}' \end{matrix}\right. }$$

理论如上。实际撸代码时候可能不是区分k存不存在分别计算亚像素边缘，而是k<1和k>=1分别计算结果。