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 = 01.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参考资料
声明
特别感谢Elby。
本文程序只做学习研究，详情请参考链接内容！
