最小二乘圆拟合: Pratt 方法

kasa方法圆拟合作为最常见的圆拟合方法，虽然计算方法简单，效率快，但是拟合结果存在较大偏差(Heavy bias)。通过将\(D=\sqrt{(x-a)^2+(y-b)^2}\)与半径R的差值\(D-R\)转换为\(D^2-R^2\)，将非线性问题转换为线性问题。但是因为\(D^2-R^2 = d(d+2R) (令d=D-R)\),当偏离值d较大时，kasa方法导致R明显变小。Pratt通过将kasa方法的目标除以\((2R)^2\)的方式，取得更准确的结果。

$$f=  \frac{ \sum_{i=1}^{n}((x_i-a)^2+(y_i-b)^2-R^2)^2}{(2R)^2}$$类似kasa method，pratt方法将圆方程描述为 \(A(x^2+y^2) + Bx + Cy+D=0\)。这样圆心\((a,b) = (-\frac{B}{2A},-\frac{C}{2A})\)半径\(R = \frac{\sqrt{B^2+C^2-4AD}}{2A}\)。注意，A有可能为0，此时圆拟合方法得到的是一条直线，或者说一个曲率为0的圆。但是在这里，因为A,B,C,D的值乘以一个标量也不会改变圆的方程，可以让A=1,所以目标方程可以改写为：

$$f= (\frac{\sum A(x_i^2+y_i^2) + Bx_i + Cy_i+D}{B^2+C^2-4AD})^2 $$

Pratt将这个方程转化为求$$g=(\sum A(x_i^2+y_i^2) + Bx_i + Cy_i+D)^2$$的最小值，并且在约束条件$$B^2+C^2-4AD=1$$的条件下。

惯例将上式写成矩阵形式:$$g= A^TMA  (\vec A=[A,B,C,D]^T ; M=\begin{bmatrix}\sum z^2 & \sum xz & \sum yz & \sum z \\ \sum xz & \sum x^2 & \sum xy & \sum x \\ \sum yz & \sum xy &\sum y^2 & \sum y \\ \sum z & \sum x & \sum y & n \end{bmatrix}; z_i = x_i^2 + y_i^2)$$

$$A^TBA = \begin{bmatrix}A&B&C&D\end{bmatrix}\begin{bmatrix}0&0&0&-2\\0&1&0&0\\0&0&1&0\\-2&0&0&0 \end{bmatrix}\begin{bmatrix}A\\B\\C\\D\end{bmatrix}=1 $$

通过引入拉格朗日乘子解决最小值问题：

$$g(A,\eta) = A^TMA – \eta(A^TBA-1)$$

对\(\vec A\)求偏导: \(MA – \eta BA=0\)。因为B可逆，可以变换为 \(B^{-1})MA = \eta A\)。可知 \(\eta\)为矩阵\(B^{-1}M\)的特征值，A为其特征向量。

由\(MA – \eta BA=0\)可得\(A^TMA – \eta A^TBA=ATMA-\eta=0\)即\(A^TMA\)的极值为\(\eta\)，那么矩阵\(B^{-1}M\)有四个特征值，最小的特征值就是\(A^TMA\)的最小值？不然。如果\(\eta<0\) ,那么拟合出的圆将没有意义，因为\(A^TMA \ge 0\)。（Nikolai Chernov 教授经过复杂的证明，\(B^{-1}M\) 存在一个负特征值。）因此，A的值应该是对应矩阵\(B^{-1}M\)最小非负特征值的特征向量。

最后，通过采用Eigen库的EigenSolver可以计算特征值和特征向量。

---

算法分析：

pratt 方法有两个明显优势，一个是消除了半径的影响，使得曲率半径非常大、噪音大的情况下算法也能正确运行。

一个是在方程式上融合了直线方程，当系数A=0时，求得的曲线是直线，在实际应用中，不清楚待拟合曲线为直线还是曲线的情况下非常有效。

还有一点就是pratt方法的效率也比较高，因此是一种非常实用的方法。

---

本文参考材料 Nikolai Chernov 教授的 < Circular and Linear Regression > Chernov 教授主页上也给出了编程实现。感兴趣的同学可以直接下载试用。