基于对一般二次曲面拟合效果的不满，特地整理这一篇文章。不加任何限制的一般二次曲面拟合在机器视觉实际应用时会出现很多意外的情况。比如文章《匹配位姿拟合求精方法 - 兜尼完 - 博客园 (cnblogs.com)》和《9点拟合梯度边缘亚像素方法 - 兜尼完 - 博客园 (cnblogs.com)》，这两种方法都来自于硕士论文，在实际应用时都会出现亚像素量超出理论范围[-0.5,0.5]的意外情况。主要原因是二次曲面的种类太多了，包含很多不同的形状，如果不施加任何限制那么拟合结果可能不符合预期。比如在梯度边缘亚像素时，你可能期望拟合一个抛物面，但结果却可能是一个马鞍面，从而使亚像素量超范围（这句话是猜测没有实际验证过）。
一、椭球面拟合
这里首先叙述椭球面拟合方法。它可以很方便地扩展成三元二次函数的拟合从而应用于匹配位姿拟合的求精中。下面是关于拟合椭球面的几篇论文：
- 最大似然估计法拟合椭球面：2110.13337.pdf (arxiv.org)，没看懂。
- 基于Cayley变换的拟合椭球面：2304.10630.pdf (arxiv.org)，没看懂。
- 基于ADMM方法的拟合椭球面：Fast multidimensional ellipsoid-specific fitting by alternating direction method of multipliers，这篇文章没搜到免费原文。
- 附加条件限制的最小二乘法拟合椭球面：(PDF) Least squares ellipsoid specific fitting (researchgate.net)，这个看懂了。
这里叙述上方第四个文章的方法。设待拟合椭球面方程为：
$${ ax^{2}+by^{2}+cz^{2}+2fyz+2gxz+2hxy+2px+2qy+2rz+d=0 }$$
我们有N组数据。用矩阵D表示样本矩阵，矩阵v表示待定系数矩阵。有：
$${ \mathbf D = \left( \begin{matrix} \mathbf X_{1} & \mathbf X_{2} & \cdots & \mathbf X_{N} \end{matrix} \right),其中\mathbf X_{i}= \left( \begin{matrix} 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{matrix} \right)^{T} }$$
$${ \mathbf v= \left( \begin{matrix} a & b & c & f & g & h & p & q & r & d \end{matrix} \right)^{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求偏导数可得：
$${ \frac{\partial e}{\partial \mathbf v}= 2\mathbf D \mathbf D^{T} \mathbf v - 2\lambda \mathbf C \mathbf v = 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 }$。
标签：right,mathbf,多项式,矩阵,椭球面,拟合,left From： https://www.cnblogs.com/mengxiangdu/p/17582190.html