圆与椭圆拟合推导
本文整理圆拟合与椭圆拟合从代数方程到几何参数的核心推导。
1. 圆拟合
圆心为 \((a,b)\)、半径为 \(r\) 的圆满足
\[ (x-a)^2+(y-b)^2=r^2. \]
展开并移项：
\[ x^2+y^2=2ax+2by+c, \qquad c=r^2-a^2-b^2. \]
因此，方程对参数 \((a,b,c)\) 是线性的。对样本点 \((x_i,y_i)\)，可写成最小二乘问题
\[ \begin{bmatrix} 2x_1&2y_1&1\\ 2x_2&2y_2&1\\ \vdots&\vdots&\vdots\\ 2x_n&2y_n&1 \end{bmatrix} \begin{bmatrix}a\\b\\c\end{bmatrix} \approx \begin{bmatrix} x_1^2+y_1^2\\ x_2^2+y_2^2\\ \vdots\\ x_n^2+y_n^2 \end{bmatrix}. \]
求出 \((a,b,c)\) 后，半径为
\[ r=\sqrt{c+a^2+b^2}. \]
设计矩阵必须满列秩，即 \(\operatorname{rank}=3\)，否则不能唯一确定圆。共线或近似共线的点通常会导致秩不足或数值不稳定；平移和缩放只能改善数值条件，不能消除这种几何退化。
2. 椭圆拟合
2.1 拟合一般二次曲线
先将椭圆写成一般二次曲线
\[ Ax^2+Bxy+Cy^2+Dx+Ey+F=0. \]
对每个样本点 \((x_i,y_i)\)，有
\[ [x_i^2,\ x_iy_i,\ y_i^2,\ x_i,\ y_i,\ 1]q=0, \qquad q=[A,B,C,D,E,F]^\mathsf{T}. \]
堆叠所有样本后得到齐次系统
\[ \mathcal Dq\approx0. \]
由于 \(q=0\) 是无意义的平凡解，需要固定参数尺度：
\[ \min_{\lVert q\rVert=1}\lVert\mathcal Dq\rVert^2. \]
若 \(\mathcal D=U\Sigma V^\mathsf{T}\)，解就是最小奇异值对应的右奇异向量，即最接近 \(\mathcal D\) 零空间的方向。
2.2 求椭圆中心
令 \(z=[x,y]^\mathsf{T}\)，将二次曲线写成
\[ f(z)=z^\mathsf{T}Mz+l^\mathsf{T}z+c, \]
其中
\[ M=\begin{bmatrix}A&B/2\\B/2&C\end{bmatrix}, \qquad l=\begin{bmatrix}D\\E\end{bmatrix}, \qquad c=F. \]
中心 \(z_0\) 是消去一次项的平移点。由
\[ \nabla f(z)=2Mz+l=0 \]
得到
\[ z_0=-\frac12M^{-1}l. \]
这里要求 \(M\) 可逆；否则二次曲线是退化的，或不具有有限中心。
2.3 平移到中心
引入以中心为原点的新坐标 \(u=[u_1,u_2]^\mathsf{T}\)：
\[ z=z_0+u. \]
将它代入原二次型：
\[ \begin{aligned} f(z_0+u) &=(z_0+u)^\mathsf{T}M(z_0+u) +l^\mathsf{T}(z_0+u)+c \\ &=z_0^\mathsf{T}Mz_0 +z_0^\mathsf{T}Mu +u^\mathsf{T}Mz_0 +u^\mathsf{T}Mu +l^\mathsf{T}z_0 +l^\mathsf{T}u+c. \end{aligned} \]
因为 \(M\) 对称，标量 \(z_0^\mathsf{T}Mu=u^\mathsf{T}Mz_0\)，所以含 \(u\) 的一次项可以合并为
\[ z_0^\mathsf{T}Mu+u^\mathsf{T}Mz_0+l^\mathsf{T}u =(2Mz_0+l)^\mathsf{T}u. \]
中心满足 \(2Mz_0+l=0\)，因此该项恰好消失，剩下
\[ f(z_0+u)=u^\mathsf{T}Mu+c', \]
其中新的常数项为
\[ c'=z_0^\mathsf{T}Mz_0+l^\mathsf{T}z_0+c. \]
再利用 \(l=-2Mz_0\)，可以得到两个等价形式：
\[ \begin{aligned} c' &=c-z_0^\mathsf{T}Mz_0 \\ &=c-\frac14l^\mathsf{T}M^{-1}l. \end{aligned} \]
因此，平移后的曲线方程为
\[ u^\mathsf{T}Mu+c'=0. \]
到这里，已经得到了以椭圆中心为原点的隐式方程。若只需要描述曲线，或判断某个点是否满足椭圆方程，可以停在这一步。但矩阵 \(M\) 中仍可能包含由原坐标轴方向引起的交叉项，因此半轴长度和主轴方向还没有被直接表示出来。
如果还需要得到半长轴、半短轴及旋转角，就要继续对 \(M\) 做特征分解，把坐标轴旋转到椭圆的主轴方向。平移只改变一次项和常数项，不改变二次项矩阵 \(M\)，所以这些几何信息仍由 \(M\) 决定。
2.4 求主轴方向与半轴长度
\(M\) 是实对称矩阵，可作特征分解
\[ M=Q\Lambda Q^\mathsf{T}, \qquad \Lambda=\operatorname{diag}(\lambda_1,\lambda_2). \]
\(Q=[q_1,q_2]\) 的两列是 \(M\) 的单位正交特征向量，因此 \(Q\) 是正交矩阵，满足
\[ Q^\mathsf{T}Q=I, \qquad Q^{-1}=Q^\mathsf{T}. \]
\(u\) 是平移到椭圆中心后的坐标，但仍沿用原坐标轴。为了消去二次型中的交叉项，将 \(u\) 改用特征向量 \(q_1,q_2\) 作为基底：
\[ u=v_1q_1+v_2q_2=Qv, \qquad v=Q^\mathsf{T}u. \]
这里的 \(v=[v_1,v_2]^\mathsf{T}\) 就是点在椭圆主轴坐标系中的坐标。这个变换只旋转坐标轴，不再进行平移。
将 \(u=Qv\) 代入平移后的二次型：
\[ \begin{aligned} u^\mathsf{T}Mu &=(Qv)^\mathsf{T}M(Qv) \\ &=v^\mathsf{T}Q^\mathsf{T}MQv \\ &=v^\mathsf{T}\Lambda v. \end{aligned} \]
因此，曲线方程变为
\[ \lambda_1v_1^2+\lambda_2v_2^2+c'=0, \]
即
\[ \frac{v_1^2}{-c'/\lambda_1} +\frac{v_2^2}{-c'/\lambda_2}=1. \]
因此，\(Q\) 的列向量给出主轴方向，对应半轴长度为
\[ a_i=\sqrt{-\frac{c'}{\lambda_i}}, \qquad i=1,2. \]
椭圆的旋转角可由选定主轴的特征向量计算。若约定 \(a\ge b\)，应先按半轴长度排序，再用长轴对应的特征向量确定角度。
3. 有效性条件
一般二次曲线的最小二乘解不一定是实椭圆。将整个方程乘以非零常数不会改变曲线，因此可先统一整体符号，使 \(M\) 为正定矩阵。随后必须满足
\[ \lambda_1>0, \qquad \lambda_2>0, \qquad c'<0. \]
这样 \(-c'/\lambda_i>0\)，两个半轴才是有限实数。若条件不成立，拟合结果可能是双曲线、抛物线、虚椭圆或退化二次曲线，不能转换为有效的几何椭圆参数。
4. 推导链条
- 圆：展开圆方程 → 线性最小二乘求圆心及辅助参数 → 恢复半径。
- 椭圆：齐次最小二乘求二次曲线 → 平移消去一次项 → 特征分解求方向和轴长 → 检查有效性。