椭球面拟合方法及一般多项式函数拟合拓展

基于对一般二次曲面拟合实际应用效果的不满，特地整理这一篇文章。不加任何限制的一般二次曲面拟合在机器视觉实际应用时会出现很多意外的情况。比如文章《匹配位姿拟合求精方法 - 兜尼完 - 博客园 (cnblogs.com)》和《9点拟合梯度边缘亚像素方法 - 兜尼完 - 博客园 (cnblogs.com)》，这两种方法都来自于硕士论文，在实际应用时都会出现亚像素量超出理论范围[-0.5,0.5]的意外情况。主要原因是二次曲面的种类太多了，包含很多不同的形状，如果不施加任何限制那么拟合结果可能不符合预期。比如在梯度边缘亚像素时，你可能期望拟合一个抛物面，但结果却可能是一个马鞍面，从而使亚像素量超范围（这句话是猜测没有实际验证过）。

一、椭球面拟合

这里首先叙述椭球面拟合方法。它可以很方便地扩展成三元二次函数的拟合从而应用于匹配位姿拟合的求精中。下面是关于拟合椭球面的几篇论文：

这里叙述上方第四个文章的方法。设待拟合椭球面方程为：

$${ ax^{2}+by^{2}+cz^{2}+2fyz+2gxz+2hxy+2px+2qy+2rz+d=0 }$$

我们有N组数据。用矩阵D表示样本矩阵，矩阵v表示待定系数矩阵。有：

$${ \mathbf D = \begin{pmatrix} \mathbf X_{1} & \mathbf X_{2} & \cdots & \mathbf X_{N} \end{pmatrix},其中\mathbf X_{i}= \begin{pmatrix} x^{2}_{i} & y^{2}_{i} &  z^{2}_{i} & 2y_{i}z_{i} & 2x_{i}z_{i} & 2x_{i}y_{i} & 2x_{i} & 2y_{i} & 2z_{i} & 1 \end{pmatrix}^{T} }$$

$${ \mathbf v= \begin{pmatrix} a & b & c & f & g & h & p & q & r & d \end{pmatrix}^{T} }$$

定义误差函数如下：

$${ e = \mathbf v^{T} \mathbf D \mathbf D^{T} \mathbf v,需要\mathbf v^{T} \mathbf C \mathbf v = 1 }$$

上式中${ \mathbf C = \begin{pmatrix} \mathbf C_{1} & \mathbf 0_{6×4} \\ \mathbf 0_{4×6} & \mathbf 0_{4×4} \end{pmatrix},其中\mathbf C_{1}=\begin{pmatrix} -1 & 0.5k-1 & 0.5k-1 & 0 & 0 & 0 \\ 0.5k-1 & -1 & 0.5k-1 & 0 & 0 & 0 \\ 0.5k-1 & 0.5k-1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -k & 0 & 0 \\ 0 & 0 & 0 & 0 & -k & 0 \\ 0 & 0 & 0 & 0 & 0 & -k \end{pmatrix} }$。对于一般的椭圆而言可取${ k=4 }$，它可以确保拟合结果是一个椭球体。需要了解的是不是所有的椭圆都满足${ k=4 }$时的条件，详情见论文原文。对误差函数加上拉格朗日乘子变为：

$${ e = \mathbf v^{T} \mathbf D \mathbf D^{T} \mathbf v - \lambda \left( \mathbf v^{T} \mathbf C \mathbf v - 1 \right) }$$

对矩阵v求偏导数可得下式。对矩阵求导的含义是对矩阵中的每一个元素分别求导。对于二次型${ \mathbf v^{T} \mathbf M \mathbf v }$而言，它关于矩阵v的导数是${ 2 \mathbf M \mathbf v }$。可以把矩阵表达式展开进行验证：

$${ \frac{\partial e}{\partial \mathbf v}= 2\mathbf D \mathbf D^{T} \mathbf v - 2\lambda \mathbf C \mathbf v = \mathbf 0,需要\mathbf v^{T} \mathbf C \mathbf v=1 }$$

由于矩阵C是没有逆矩阵的，所以可以做如下变换使上式变成矩阵特征值和特征向量的形式：

$${ \frac{1}{\lambda} \mathbf v =  \left( \mathbf D \mathbf D^{T} \right)^{-1} \mathbf C \mathbf v }$$

上式的解${ \mathbf v }$就是右侧系数矩阵${ \left( \mathbf D \mathbf D^{T} \right)^{-1} \mathbf C }$的最大正特征值对应的特征向量。注意这里是最大特征值不是最小特征值哦。另外矩阵${ \mathbf D \mathbf D^{T} }$也不一定有逆矩阵。可以对矩阵D按照矩阵C的形式分为4个分块矩阵进行更细一点的处理以提高可求逆的概率。详情可参考论文原文。

二、一般多项式函数拟合的拓展

对于一般的多项式函数有多种拟合方法。设待拟合的多项式函数如下：

$${ f(x,y,\cdots)=c_{1}x^{p} + c_{2}y^{p}+ \cdots + c_{M} = 0 }$$

我们有N组数据。用矩阵D表示采集数据的样本矩阵，矩阵x表示待定系数${c_{i} }$组成的列向量。

$${ \mathbf D=\begin{pmatrix} x_{1}^{p} & y_{1}^{p} & \cdots & 1 \\ x_{2}^{p} & y_{2}^{p} & \cdots & 1 \\ \vdots & \vdots & \ddots & 1 \\ x_{N}^{p} & y_{N}^{p} & \cdots & 1 \end{pmatrix},\mathbf x=\begin{pmatrix} c_{1} \\ c_{2} \\ \vdots \\ c_{M} \end{pmatrix} }$$

此时，我们要求解的方程组可以写成矩阵形式：

$${ \mathbf D \mathbf x = \mathbf 0 }$$

它是齐次线性方程组，不能直接使用最小二乘法。但可以增加限定条件将其转换成最小二乘问题，一般有如下两种方法求最优解。

1）代数距离法

这种求解方法一般被称为代数距离（algebraic distance）法。其误差函数如下：

$${ e=\left\| \mathbf D \mathbf x \right\|^{2} = \mathbf x^{T} \mathbf D^{T} \mathbf D \mathbf x,需要\mathbf x^{T} \mathbf x=1 }$$

可以看出来上式就是一个二次型求最小值的问题。它的解就是矩阵${ \mathbf D^{T} \mathbf D }$的最小特征值对应的特征向量。需要注意的是矩阵${ \mathbf D^{T} \mathbf D }$可能没有逆矩阵，此时需要用广义逆矩阵代替。在OpenCV中可以使用SVD分解求广义逆矩阵。

2）几何距离法

这种求解方法一般被称为几何距离（geometric distance）法。其误差函数如下：

$${ e=\frac{\left\| \mathbf D \mathbf x \right\|^{2} }{ \left\| \nabla f(x,y,\cdots) \right\|^{2} } }$$

这个误差函数直接求解会比较困难。所以为了简化问题，一般会限定分母等于1。此时问题转换成等式限制的最小二乘问题，求解方式和代数距离法相似。具体例子可以参考论文《37401.37420 (acm.org)》，此文第149页讲述了拟合圆的方法，在限定圆的一般式的梯度模长等于1的情况下得到等式条件${ D^{2}+ E^{2}-4AF=1 }$。

posted @ 2023-07-27 20:11  兜尼完  阅读(730)  评论(0)  编辑  收藏  举报