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 教授主页上也给出了编程实现。感兴趣的同学可以直接下载试用。