最小二乘圆

日期2011-03-22（最后修改），2011-03-19（创建）

简介

此页面收集了用于查找拟合一组 2D 点 (x,y) 的最小二乘圆的不同方法。

此分析的完整代码可在此处获取：least_squares_circle_v1d.py。

查找最小二乘圆对应于查找圆心 (xc, yc) 及其半径 Rc，它们使下面定义的残差函数最小化

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

这是一个非线性问题。我们将看到三种解决问题的方法，并比较它们的结果以及速度。

使用代数近似

如此文档中所述，如果我们将要最小化的函数定义如下，则此问题可以用线性问题来近似

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

这导致了以下方法，使用 linalg.solve

在 [ ]
#! python
# == METHOD 1 ==
method_1 = 'algebraic'

# coordinates of the barycenter
x_m = mean(x)
y_m = mean(y)

# calculation of the reduced coordinates
u = x - x_m
v = y - y_m

# linear system defining the center (uc, vc) in reduced coordinates:
#    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)

# Solving the linear system
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

# Calcul des distances au centre (xc_1, yc_1)
Ri_1     = sqrt((x-xc_1)**2 + (y-yc_1)**2)
R_1      = mean(Ri_1)
residu_1 = sum((Ri_1-R_1)**2)

使用 scipy.optimize.leastsq

Scipy 附带了一些用于解决上述非线性问题的工具。其中，scipy.optimize.leastsq 在这种情况下非常易于使用。

实际上，一旦定义了圆心，半径可以直接计算出来，并且等于 mean(Ri)。因此，只剩下两个参数：xc 和 yc。

基本用法

在 [ ]
#! python
#  == METHOD 2 ==
from scipy      import optimize

method_2 = "leastsq"

def calc_R(xc, yc):
    """ calculate the distance of each 2D points from the center (xc, yc) """
    return sqrt((x-xc)**2 + (y-yc)**2)

def f_2(c):
    """ calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """
    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)

高级用法，使用雅可比函数

为了提高速度，可以通过添加 Dfun 参数来告诉 optimize.leastsq 如何计算函数的雅可比矩阵

在 [ ]
#! python
# == METHOD 2b ==
method_2b  = "leastsq with jacobian"

def calc_R(xc, yc):
    """ calculate the distance of each data points from the center (xc, yc) """
    return sqrt((x-xc)**2 + (y-yc)**2)

def f_2b(c):
    """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
    Ri = calc_R(*c)
    return Ri - Ri.mean()

def Df_2b(c):
    """ Jacobian of f_2b
    The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
    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)

使用 scipy.odr

Scipy 有一个专门用于处理正交距离回归的包，即 scipy.odr。这个包可以处理显式和隐式函数定义，在本例中我们将使用第二种形式。

以下是圆的隐式定义

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

基本用法

这将导致以下代码

在 [ ]
#! python
# == METHOD 3 ==
from scipy      import  odr

method_3 = "odr"

def f_3(beta, x):
    """ implicit definition of the circle """
    return (x[0]-beta[0])**2 + (x[1]-beta[1])**2 -beta[2]**2

# initial guess for parameters
R_m = calc_R(x_m, y_m).mean()
beta0 = [ x_m, y_m, R_m]

# for implicit function :
#       data.x contains both coordinates of the points (data.x = [x, y])
#       data.y is the dimensionality of the response
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)

高级用法，使用雅可比函数

隐式函数定义的优点之一是它的导数非常容易计算。

这可以用来完成模型

在 [ ]
#! python
# == METHOD 3b ==
method_3b  = "odr with jacobian"

def f_3b(beta, x):
    """ implicit definition of the circle """
    return (x[0]-beta[0])**2 + (x[1]-beta[1])**2 -beta[2]**2

def jacb(beta, x):
    """ Jacobian function with respect to the parameters beta.
    return df_3b/dbeta
    """
    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):
    """ Jacobian function with respect to the input x.
    return 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):
    """ Return a first estimation on the parameter from the 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

# for implicit function :
#       data.x contains both coordinates of the points
#       data.y is the dimensionality of the response
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)    # beta0 has been replaced by an estimate function
lsc_odr.set_job(deriv=3)                    # use user derivatives function without checking
lsc_odr.set_iprint(iter=1, iter_step=1)     # print details for each iteration
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)

三种方法的比较

我们将比较这三种方法在两种情况下的结果

* 当 2D 点 都在 圆 周围 时

* 当 2D 点 在一个 小 弧 上 时

数据点都在圆周围

以下是一个数据点都在圆周围的示例

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

这将给出

||||||||||||||||<tablewidth="100%">摘要|| ||方法|| Xc || Yc || Rc ||nb_calls || std(Ri)|| residu || residu2 || ||代数 || 10.55231 || 9.70590|| 23.33482|| 1|| 1.135135|| 7.731195|| 16236.34|| ||leastsq || 10.50009 || 9.65995|| 23.33353|| 15|| 1.133715|| 7.711852|| 16276.89|| ||leastsq with jacobian || 10.50009 || 9.65995|| 23.33353|| 7|| 1.133715|| 7.711852|| 16276.89|| ||odr || 10.50009 || 9.65995|| 23.33353|| 82|| 1.133715|| 7.711852|| 16276.89|| ||odr with jacobian || 10.50009 || 9.65995|| 23.33353|| 16|| 1.133715|| 7.711852|| 16276.89||

注意

* `nb_calls` 对应于 要 最小化 的 函数 的 调用 次数， 不 考虑 对 导数 函数 的 调用 次数。 这 与 迭代 次数 不同， 因为 ODR 可以在 一次 迭代 中 进行 多次 调用。

* 如 以下 图形 所示， 两个 函数 `residu` 和 `residu_2` 并不 等效， 但 它们的 最小值 在这种 情况下 很 接近。

数据点在弧周围

以下是一个数据点形成弧的示例

在 [ ]
#! python
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
leastsq 9.88760 3.68753 27.35040 24 0.820825 4.042522 12001.98
leastsq with jacobian 9.88759 3.68752 27.35041 10 0.820825 4.042522 12001.98
odr 9.88757 3.68750 27.35044 472 0.820825 4.042522 12002.01
odr with jacobian 9.88757 3.68750 27.35044 109 0.820825 4.042522 12002.01

arc_v5.png arc_residu2_v6.png

结论

ODR 和 leastsq 给出了相同的结果。

Optimize.leastsq 是最有效的方法，至少在函数调用次数方面，它可以比 ODR 快两到十倍。

添加一个计算雅可比矩阵的函数可以将函数调用次数减少两到五倍。

在这种情况下，使用 ODR 似乎有点过分，但对于更复杂的用例（如省略号）来说，它非常方便。

当所有点都围绕圆形时，代数近似法可以得到良好的结果，但当只有弧线要拟合时，它就受到限制。

事实上，当数据点不完全位于圆上时，要最小化的两个误差函数并不等效。代数方法在大多数情况下会导致比最小二乘圆更小的半径，因为它的误差函数基于平方距离，而不是距离本身。

章节作者：Elby

附件