Python 圆拟合

发布于:2023-01-24 ⋅ 阅读:(26) ⋅ 点赞:(0) ⋅ 评论:(0)

上一篇文章介绍了关于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


网站公告

欢迎关注微信公众号

今日签到

点亮在社区的每一天
签到