代数逼近法、最小二乘法、正交距离回归法来拟合圆及其结果对比（Python）

0 引言

在进行动态跟踪时，有时可能会关注轨迹的运动状态，例如获取沿圆弧轨迹运动物体的运动半径大小。本文介绍了几种算法对点集（xi，yi）进行圆拟合的方法：代数逼近法、最小二乘法和正交距离回归法。
其中，最常用的是最小二乘法，求最小二乘法的圆就是求圆心（xc，yc）和其半径Rc，使残余函数最小，残余函数定义如下：

#! python
Ri = sqrt( (x - xc)**2 + (y - yc)**2)
residu = sum( (Ri - Rc)**2)

这是一个非线性问题，本文将采用几种方法来解决，并对比其运算结果和速度。
特别指出的是，在进行圆拟合时，分两种情况进行对比分析，即数据点呈圆分布、数据点呈圆弧分布（圆的一小部分），两者有不同，注意区分。
本文着重介绍算法的Python程序实现过程与结果，对于理论部分，读者可参考相关链接或其他资料。

1 拟合圆的算法

1.1 代数逼近法

如果作如下定义的残余函数，代数逼近算法可以被看作是线性问题，详细过程参看此文档：

#! python
residu_2 = sum( (Ri**2 - Rc**2)**2)

使用linalg.solve：

#! python
# == 方法 1 ==
method_1 = 'algebraic'

# 质心坐标
x_m = mean(x)
y_m = mean(y)

# 相对坐标
u = x - x_m
v = y - y_m

# 相对坐标下定义中心(uc, vc)的线性系统：
#    Suu * uc +  Suv * vc = (Suuu + Suvv)/2
#    Suv * uc +  Svv * vc = (Suuv + Svvv)/2
Suv  = sum(u*v)
Suu  = sum(u**2)
Svv  = sum(v**2)
Suuv = sum(u**2 * v)
Suvv = sum(u * v**2)
Suuu = sum(u**3)
Svvv = sum(v**3)

# 求线性系统
A = array([ [ Suu, Suv ], [Suv, Svv]])
B = array([ Suuu + Suvv, Svvv + Suuv ])/2.0
uc, vc = linalg.solve(A, B)

xc_1 = x_m + uc
yc_1 = y_m + vc

#半径残余函数
Ri_1     = sqrt((x-xc_1)**2 + (y-yc_1)**2)
R_1      = mean(Ri_1)
residu_1 = sum((Ri_1-R_1)**2)

程序中代码比较直观，读者可根据代码写出相应的计算公式，这样便于理解。

1.2 最小二乘法

Scipy库提供了一些用于解决非线性问题的工具，其中，本文使用的scipy.optimize.leastsq便是非常简单的一个，可解决最小二乘法的问题。
实际上，一但圆心确定，半径便可直接计算，其等于半径的均值，即mean(Ri)。因此，这只考虑两个参数，即圆心坐标：xc、yc。

1.2.1 基本用法

#! python
#  == 最小二乘法 ==
from scipy import optimize

method_2 = "leastsq"

def calc_R(xc, yc):
    """ 计算s数据点与圆心(xc, yc)的距离 """
    return sqrt((x-xc)**2 + (y-yc)**2)

def f_2(c):
    """ 计算半径残余"""
    Ri = calc_R(*c)
    return Ri - Ri.mean()

center_estimate = x_m, y_m
center_2, ier = optimize.leastsq(f_2, center_estimate)

xc_2, yc_2 = center_2
Ri_2       = calc_R(*center_2)
R_2        = Ri_2.mean()
residu_2   = sum((Ri_2 - R_2)**2)

1.2.2 基于Jacobian函数的最小二乘法

为获取计算速率，通过增加一个参数，便可采用optimize.leastsq计算雅可比矩阵，如下：

#! python
# == 基于Jacobian函数的最小二乘法 ==
method_2b  = "leastsq with jacobian"

def calc_R(xc, yc):
    """ 计算数据点据圆心(xc, yc)的距离 """
    return sqrt((x-xc)**2 + (y-yc)**2)

def f_2b(c):
    """ 计算数据点与以c=(xc, yc)为圆心圆的半径差值 """
    Ri = calc_R(*c)
    return Ri - Ri.mean()

def Df_2b(c):
    """ 雅可比矩阵"""
    xc, yc     = c
    df2b_dc    = empty((len(c), x.size))

    Ri = calc_R(xc, yc)
    df2b_dc[0] = (xc - x)/Ri                   # dR/dxc
    df2b_dc[1] = (yc - y)/Ri                   # dR/dyc
    df2b_dc    = df2b_dc - df2b_dc.mean(axis=1)[:, newaxis]

    return df2b_dc

center_estimate = x_m, y_m
center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True)

xc_2b, yc_2b = center_2b
Ri_2b        = calc_R(*center_2b)
R_2b         = Ri_2b.mean()
residu_2b    = sum((Ri_2b - R_2b)**2)

1.3 正交距离回归法

Scipy提供了一个专用包scipy.odr来解决正交距离回归问题。该软件包可以处理显式函数定义和隐式函数定义，本文采用后者。
圆的隐式函数定义如下：

#! python
(x - xc)**2 + (y-yc)**2 - Rc**2 = 0

1.3.1 基本用法

代码如下：

#! python
# == 正交距离回归法 ==
from scipy import odr

method_3 = "odr"

def f_3(beta, x):
    """ 圆的隐式定义 """
    return (x[0]-beta[0])**2 + (x[1]-beta[1])**2 -beta[2]**2

"""参数初始化"""
R_m = calc_R(x_m, y_m).mean()
beta0 = [ x_m, y_m, R_m]

lsc_data  = odr.Data(row_stack([x, y]), y=1)
lsc_model = odr.Model(f_3, implicit=True)
lsc_odr   = odr.ODR(lsc_data, lsc_model, beta0)
lsc_out   = lsc_odr.run()

xc_3, yc_3, R_3 = lsc_out.beta
Ri_3 = calc_R([xc_3, yc_3])
residu_3 = sum((Ri_3 - R_3)**2)

1.3.2 基于Jacobian函数的正交距离回归法

隐式函数定义便于求解倒数。
计算模型如下：

#! python
# == 基于Jacobian函数的正交距离回归法 ==
method_3b  = "odr with jacobian"

def f_3b(beta, x):
    """ 圆定义 """
    return (x[0]-beta[0])**2 + (x[1]-beta[1])**2 -beta[2]**2

def jacb(beta, x):
    """ 计算关于参数β的雅可比函数，返回df_3b/dβ。"""
    xc, yc, r = beta
    xi, yi    = x

    df_db    = empty((beta.size, x.shape[1]))
    df_db[0] =  2*(xc-xi)                     # d_f/dxc
    df_db[1] =  2*(yc-yi)                     # d_f/dyc
    df_db[2] = -2*r                           # d_f/dr

    return df_db

def jacd(beta, x):
    """ 计算关于输入x的雅可比函数，返回df_3b/dx。"""
    xc, yc, r = beta
    xi, yi    = x

    df_dx    = empty_like(x)
    df_dx[0] =  2*(xi-xc)                     # d_f/dxi
    df_dx[1] =  2*(yi-yc)                     # d_f/dyi

    return df_dx

def calc_estimate(data):
    """ 从数据中返回对参数的第一次估计。 """
    xc0, yc0 = data.x.mean(axis=1)
    r0 = sqrt((data.x[0]-xc0)**2 +(data.x[1] -yc0)**2).mean()
    return xc0, yc0, r0

lsc_data  = odr.Data(row_stack([x, y]), y=1)
lsc_model = odr.Model(f_3b, implicit=True, estimate=calc_estimate, fjacd=jacd, fjacb=jacb)
lsc_odr   = odr.ODR(lsc_data, lsc_model)
lsc_odr.set_job(deriv=3)
lsc_odr.set_iprint(iter=1, iter_step=1)
lsc_out   = lsc_odr.run()

xc_3b, yc_3b, R_3b = lsc_out.beta
Ri_3b       = calc_R(xc_3b, yc_3b)
residu_3b   = sum((Ri_3b - R_3b)**2)

2 算法对比

本文分两种情况来对比算法，即：
（1）数据点在圆上的情况，即绕圆的数据点；
（2）数据点在圆的一段圆弧上的情况，即绕圆弧的数据点。

2.1 绕圆的数据点

有以下绕圆的数据点：

#! python
# Coordinates of the 2D points
x = r_[  9,  35, -13,  10,  23,   0]
y = r_[ 34,  10,   6, -14,  27, -10]

可得到

方法	xc	yc	Rc	nb_calls	std(Ri)	residu	residu2
代数逼近法	10.55231	9.70590	23.33482	1	1.135135	7.731195	16236.34
最小二乘法	10.50009	9.65995	23.33353	15	1.133715	7.711852	16276.89
基于Jacobian的最小二乘法	10.50009	9.65995	23.33353	7	1.133715	7.711852	16276.89
正交距离回归法	10.50009	9.65995	23.33353	82	1.133715	7.711852	16276.89
基于Jacobian的正交距离回归法	10.50009	9.65995	23.33353	16	1.133715	7.711852	16276.89

注：
*“nb_calls”对应于要最小化的函数调用数，而不考虑对导数函数的调用次数。这与迭代次数不同，因为ODR可以在迭代期间进行多次调用。
*如下图所示，这两个函数“residu”和“residu2”不是等价的，但它们的最小值在这种情况下是非常接近的。
运行结果

2.2 绕圆弧的数据点

有以下绕圆弧的数据点：

#! python
# Coordinates of the 2D points
x = r_[36, 36, 19, 18, 33, 26]
y = r_[14, 10, 28, 31, 18, 26]

可得到

方法	xc	yc	Rc	nb_calls	std(Ri)	residu	residu2
代数逼近法	15.05503	8.83615	20.82995	1	0.930508	5.195076	9153.40
最小二乘法	9.88760	3.68753	27.35040	24	0.820825	4.042522	12001.98
基于Jacobian的最小二乘法	9.88759	3.68752	27.35041	10	0.820825	4.042522	12001.98
正交距离回归法	9.88757	3.68750	27.35044	472	0.820825	4.042522	12002.01
基于Jacobian的正交距离回归法	9.88757	3.68750	27.35044	109	0.820825	4.042522	12002.01

运行结果

3 结论

正交距离回归算法（ODR）和最小二乘法可得到相同的结果。最小二乘法（Optimize.leastsq）是最有效的方法，运算速度是ODR算法的2~10倍，其与函数调用的次数也是相关的。
添加一个计算雅可比的函数可以使函数调用的次数减少2到5倍。
此情况下，正交距离回归算法（ODR）算法看似不必要，但对于更复杂的情况（如椭圆）来说，是非常方便的。
当数据点绕圆分布时，代数逼近算法可得到较好的结果；但当数据点只在一段圆弧分布时，此算法便有较大的局限性。
实际上，当数据点并非都呈圆分布时，用于求最小值的两个误差函数是不一样的。在大部分情况下，代数逼近算法比最小二乘算法得到的圆半径更小，因为代数逼近算法的误差是基于距离的平方值来计算，而不是直接用距离计算。

4参考资料

[1] https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html#Advanced-usage,-with-jacobian-functions

声明

特别感谢Elby。
为更直观理解，下文介绍了圆的拟合过程。<该文给出了完整的代码案例，若无法打开，可留言、私信或邮箱获取！>
本文程序只做学习研究，详情请参考链接内容！
如有疑问，可留言反馈或通过邮件交流（hardenqiu@foxmail.com）。