用Python在二维平面拟合圆

输入离散点的数据

将在excel工作簿中的离散点数据录入xi，yi列表

import xlwings as xw

count = 8 # 拟合点的数量
#将点的XY坐标数据录入xi，yi列表
wb = xw.Book('path')
sht = wb.sheets["二维平面拟合圆"]
xi = []
yi = []
for i in range(count): 
    xi.append(sht.range((i+2,2)).value) # xlwings
    yi.append(sht.range((i+2,3)).value)

拟合方法

首先假设拟合的圆的圆心为点（A,B），半径为R，那么其方程为

拟合后的圆离各离散点的距离之和应为最小，考虑使用导数=0进行求解
各点距圆的距离之和：

上式值最小等价于下式值最小：



令Q=上式，Q最小则满足Q对a，b，c的偏微分等于0，即：

即：

联立矩阵，得：

求解矩阵可得a,b,c

用Python实现

import numpy as np
import operator
from functools import reduce

sum_x = sum(xi)
sum_x2 = sum(list(i*i for i in xi))
sum_y = sum(yi)
sum_y2 = sum(list(i*i for i in yi))
sum_xy = sum(list(i*j for (i,j) in zip(xi,yi)))
sum_x2y = sum(list(i*i*j for (i,j) in zip(xi,yi)))
sum_xy2 = sum(list(i*j*j for (i,j) in zip(xi,yi)))
sum_x3 = sum(list(i*i*i for i in xi))
sum_y3 = sum(list(i*i*i for i in yi))
jz = np.array([[sum_x2,sum_xy,sum_x],[sum_xy,sum_y2,sum_y],[sum_x,sum_y,count]])
njz = np.linalg.inv(jz)#求矩阵的逆矩阵
jie = njz.dot([[-(sum_x3 + sum_xy2)],[-(sum_y3 + sum_x2y)],[-(sum_x2 + sum_y2)]])#两边同乘逆矩阵得到解
jie = jie.tolist()#先将numpy矩阵转换为嵌套列表
jie =reduce(operator.add,jie)#使用reduce转换为一维列表
a,b,c = jie
A = -a/2
B = -b/2
R = 0.5 * np.sqrt(a*a+b*b-4*c) 
print("拟合圆的方程为(x"+"{:+.2f}".format(-A)+")^2+(y"+"{:+.2f}".format(-B)+")^2 = "+f"{round(R,2)}^2")
#计算圆度,与拟合圆最大的距离减去最小的距离
yuandu = []
for (i,j) in zip(xi,yi):
    yuandu.append(np.sqrt((i-A)**2+(j-B)**2)-R)

画图显示

import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 正常显示中文字体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号

#画图查看实际拟合效果
fig = plt.figure()
axes = fig.add_subplot(111)
axes.axis('equal')
plt.scatter(xi,yi)
theta = np.arange(0, 2*np.pi, 0.01)
cir_x = A + R*np.cos(theta)
cir_y = B + R*np.sin(theta)
plt.plot(cir_x,cir_y,label = '圆度：'+str(round(max(yuandu)-min(yuandu),2)))
plt.legend()
plt.title("拟合圆的方程为(x"+"{:+.2f}".format(-A)+")^2+(y"+"{:+.2f}".format(-B)+")^2 = "+f"{round(R,2)}^2") #设置标题
plt.show()

结果展示