用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()