上一篇文章介绍了关于OPENCV椭圆拟合的方法 (Python opencv 椭圆拟合)。
在实际工程中，有些场景是需要精确的正圆测量和拟合，使用通用椭圆拟合方式容易引入较大的误差。这一篇文章将介绍关于正圆拟合的方法以及与椭圆拟合算法的对比。

实验数据

我们用sklearn里面的数据集来生成待拟合圆的数据。

    from sklearn.datasets import make_circles
    import numpy as np
    
    def get_unit_circle(shuffle=False, noise=0.02, factor=0.05):
        data = make_circles(shuffle=shuffle, noise=noise, factor=factor)
        idx = np.argwhere(data[1] == 0)
        x = data[0][idx, 0]
        y = data[0][idx, 1]
        return x, y

make_circles会生成内圆和外圆两组数据，对应的标签y分别1，0.
圆心为（0，0）半径为1.
可视化生成的数据

import matplotlib.pyplot as plt

x, y = get_unit_circle()
plt.scatter(x, y)
plt.axis('equal')
plt.show()

圆拟合

在这一步，我们会用SciPy中的optimize.leastsq去实现圆拟合，也就是调用的最小二乘法。

• Cost Function

  from scipy import optimize
  from math import pi
  
  def r(x, y, xc, yc):
  	'''
  	计算每一个点到圆心的距离
  	'''
      return np.sqrt((x-xc)**2 + (y-yc)**2)
  
  
  def f(c, x, y):
  	'''
  	Cost Function	
  	'''
      Ri = r(x, y, *c)
      return np.square(Ri - Ri.mean())
  

• 最小二乘法拟合圆

  def least_squares_circle(coords):
      """
      Circle fit using least-squares solver.
      Inputs:
  
          - coords, list or numpy array with len>2 of the form:
          [
      [x_coord, y_coord],
      ...,
      [x_coord, y_coord]
      ]
          or numpy array of shape (n, 2)
  
      Outputs:
  
          - xc : x-coordinate of solution center (float)
          - yc : y-coordinate of solution center (float)
          - R : Radius of solution (float)
          - residu : MSE of solution against training data (float)
      """
  
      x, y = None, None
      if isinstance(coords, np.ndarray):
          x = coords[:, 0]
          y = coords[:, 1]
      elif isinstance(coords, list):
          x = np.array([point[0] for point in coords])
          y = np.array([point[1] for point in coords])
      else:
          raise Exception("Parameter 'coords' is an unsupported type: " + str(type(coords)))
  
      # coordinates of the barycenter
      x_m = np.mean(x)
      y_m = np.mean(y)
      center_estimate = x_m, y_m
      center, _ = optimize.leastsq(f, center_estimate, args=(x, y))
      xc, yc = center
      Ri       = r(x, y, *center)
      R        = Ri.mean()
      residu   = np.sum((Ri - R)**2)
      return xc, yc, R, residu
  

• 绘图

  def plot_data_circle(x, y, xc, yc, R):
      """
      Plot data and a fitted circle.
      Inputs:
  
          x : data, x values (array)
          y : data, y values (array)
          xc : fit circle center (x-value) (float)
          yc : fit circle center (y-value) (float)
          R : fir circle radius (float)
  
      Output:
          None (generates matplotlib plot).
      """
      f = plt.figure(facecolor='white')
      plt.axis('equal')
  
      theta_fit = np.linspace(-pi, pi, 180)
  
      x_fit = xc + R*np.cos(theta_fit)
      y_fit = yc + R*np.sin(theta_fit)
      plt.plot(x_fit, y_fit, 'b-' , label="fitted circle", lw=2)
      plt.plot([xc], [yc], 'bD', mec='y', mew=1)
      plt.xlabel('x')
      plt.ylabel('y')
      # plot data
      plt.scatter(x, y, c='red', label='data')
  
      plt.legend(loc='best', labelspacing=0.1 )
      plt.grid()
      plt.title('Fit Circle')
  

• 圆拟合实验

  coords = np.array([[x[i][0], y[i][0]] for i in range(len(x))])
  xc, yc, r, residual = least_squares_circle(coords)
  print("least_squares: \n"
        "xc: {xc}\n"
        "yc: {yc}\n"
        "r: {r}\n"
        "residual: {residual}\n".format(xc=xc, yc=yc, r=r, residual=residual))
  plot_data_circle(coords[:, 0], coords[:, 1], xc, yc, r)
  plt.show()
  

  least_squares: 
  xc: -0.001448768226716725
  yc: 0.001868231073519271
  r: 1.0026485302748585
  s: 0.01798191226855169
  

  

对比实验

• 整圆数据拟合测试
  input

  [[-0.20222279  0.97036383]
   [-0.70938247  0.6772127 ]
   [ 0.99928194 -0.13865634]
   [ 0.43164868  0.92726211]
   [-0.98573612 -0.26547395]
   [-0.0846844   1.03034261]
   [ 0.81541324 -0.62520304]
   [ 0.63680427  0.81349025]
   [ 0.5243403   0.82797391]
   [-1.03298498 -0.13068717]
   [ 0.05791061 -1.00685492]
   [ 0.53883561 -0.85298628]
   [-0.5236286  -0.84472701]
   [-0.79429015  0.60075032]
   [-0.51160761  0.83311102]
   [-0.80068973 -0.57344023]
   [-0.71966203 -0.67572335]
   [ 0.78423462  0.57306241]
   [-0.90375463 -0.3620828 ]
   [ 1.01525985  0.25229142]]
  

  output
  有整圆的数据情况下，椭圆拟合算法和圆拟合算法得到的结果基本相同，精度都还不错。

  least_squares: 
  xc: 0.004072613875040401
  yc: -0.0018972548200443362
  r: 1.0049178131436773
  residual: 0.010274593350489435
  
  fitEllipse:
  xc: 0.010946060180664062
  yc: -0.00024581146240234376
  a: 1.998855224609375
  b: 2.0251834716796875
  

• 半圆数据拟合测试
  input

  [[ 1.00515793 -0.00581373]
   [ 0.97331632  0.14824388]
   [ 0.98822774  0.23520662]
   [ 0.95820935  0.36923693]
   [ 0.90629294  0.47303492]
   [ 0.81419798  0.56982752]
   [ 0.73716614  0.71245408]
   [ 0.61164932  0.76321955]
   [ 0.5263777   0.86008591]
   [ 0.37880738  0.90808892]
   [ 0.32055051  0.93700053]
   [ 0.20453417  0.99333937]
   [ 0.05113706  1.00244991]
   [-0.04754568  1.01810976]
   [-0.17381077  0.98401702]
   [-0.33780284  0.93704433]
   [-0.4102824   0.92013996]
   [-0.55344694  0.88235196]
   [-0.6336433   0.78243172]
   [-0.70035358  0.68997356]]
  

  output
  对于只有部分圆的数据点，圆拟合的算法比较稳定，估计圆的半径和中心都还算精确。
  椭圆拟合则会有比较大的误差，中心位置估计还算准确，但是半径的误差显著增加。

  least_squares: 
  xc: -0.005765623679872315
  yc: -0.006110242659838036
  r: 1.011230956382947
  residual: 0.005535104556472458
  
  fitEllipse:
  xc: 0.0839453125
  yc: 0.19866604614257813
  a: 1.5608851318359376
  b: 1.8778536376953125
  

  

Reference

[1] https://scipy-cookbook.readthedocs.io/index.html
[2] https://se.mathworks.com/help/curvefit/least-squares-fitting.html