代数逼近法、最小二乘法、正交距离回归法来拟合圆及其结果对比（Python） |
您所在的位置：网站首页 › 函数拟合方法有哪几种 › 代数逼近法、最小二乘法、正交距离回归法来拟合圆及其结果对比（Python） |
代数逼近法、最小二乘法、正交距离回归法来拟合圆及其结果对比（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]可得到 方法xcycRcnb_callsstd(Ri)residuresidu2代数逼近法10.552319.7059023.3348211.1351357.73119516236.34最小二乘法10.500099.6599523.33353151.1337157.71185216276.89基于Jacobian的最小二乘法10.500099.6599523.3335371.1337157.71185216276.89正交距离回归法10.500099.6599523.33353821.1337157.71185216276.89基于Jacobian的正交距离回归法10.500099.6599523.33353161.1337157.71185216276.89注： *“nb_calls”对应于要最小化的函数调用数，而不考虑对导数函数的调用次数。这与迭代次数不同，因为ODR可以在迭代期间进行多次调用。 *如下图所示，这两个函数“residu”和“residu2”不是等价的，但它们的最小值在这种情况下是非常接近的。 运行结果 有以下绕圆弧的数据点： #! python # Coordinates of the 2D points x = r_[36, 36, 19, 18, 33, 26] y = r_[14, 10, 28, 31, 18, 26]可得到 方法xcycRcnb_callsstd(Ri)residuresidu2代数逼近法15.055038.8361520.8299510.9305085.1950769153.40最小二乘法9.887603.6875327.35040240.8208254.04252212001.98基于Jacobian的最小二乘法9.887593.6875227.35041100.8208254.04252212001.98正交距离回归法9.887573.6875027.350444720.8208254.04252212002.01基于Jacobian的正交距离回归法9.887573.6875027.350441090.8208254.04252212002.01运行结果 正交距离回归算法（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。 为更直观理解，下文介绍了圆的拟合过程。 本文程序只做学习研究，详情请参考链接内容！ 如有疑问，可留言反馈或通过邮件交流（[email protected]）。 |
【本文地址】