将曲线拟合到图形上的特定范围
Fit a curve to a specific range on a graph
我从轮廓仪收集数据并尝试将收集的数据拟合到理想曲线。目的是获取曲率半径。
我编写的代码将在所有原始数据上拟合一条曲线。
适合整个范围
问题是想指定一个范围来做拟合。即：所有轮廓仪数据都不适合拟合。在下面的例子中，我想适应从 x= 180 到 x = 380
在指定范围内拟合
当然，点击指定 x_start 和 x_end 会更有帮助，但我未能从 onclick
中提取数据
奖励：如果我可以在更改拟合范围时显示 R²，那就太棒了！
感谢您的帮助
这是我想出的代码：
#importing modules
import tkinter.filedialog as tk
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
from matplotlib.widgets import Slider
#Graphic properties
fig, ax = plt.subplots()
plt.subplots_adjust(left=0.15, bottom=0.45)
#Search and open desired file:
def openfile():
opennm = tk.askopenfile()
f = open(opennm.name,"r")
data=np.genfromtxt(f,
skip_header=6,
names=True,
dtype=None,
delimiter=',')
x=[]
y=[]
for i in range(0,len(data)-1):
x.append(data[i][0])
y.append(data[i][1])
return x,y
#after opening the file, x and y data are plotted
x,y=openfile()
k,= plt.plot(x,y)
#Define the fitting function
def func(x,ROC,x_shift,y_shift):
return ((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
popt, pcov = curve_fit(func, x, y, p0=[45,250,0.8],bounds=((1,100,0),(100,300,2))) #, bounds=((-np.inf,10**-8,-np.inf,-np.inf),(np.inf,3.,np.inf,np.inf))
#The fitted curve will be overlaid on raw data
l, = plt.plot(x,func(x,*popt),'r--')
list_of_para=['ROC','X shift','Y shift']
list_of_units=['mm','a.u.','um']
for i in range (0,3):
print (list_of_para[i],'= ',"%.2f" %popt[i],' ', list_of_units[i])
#definig initial values
ROC0 = popt[0]
x_shift0 = popt[1]
y_shift0 = popt[2]
#defining sliders
axcolor = 'lightgoldenrodyellow'
axROC= plt.axes([0.15, 0.05, 0.75, 0.02], facecolor=axcolor)
axx_shift= plt.axes([0.15, 0.15, 0.75, 0.02], facecolor=axcolor)
axy_shift= plt.axes([0.15, 0.25, 0.75, 0.02], facecolor=axcolor)
#define slider limites
sROC = Slider(axROC, 'ROC', popt[0]-np.absolute(popt[0]),popt[0]+np.absolute(popt[0]), valinit=ROC0)
sx_shift = Slider(axx_shift, 'X Shift', popt[1]-np.absolute(popt[1]), popt[1]+np.absolute(popt[1]), valinit=x_shift0)
sy_shift = Slider(axy_shift, 'Y shift', popt[2]-np.absolute(popt[2]), popt[2]+np.absolute(popt[2]), valinit=y_shift0)
#define slider update values
def update(val):
ROC = sROC.val
x_shift = sx_shift.val
y_shift = sy_shift.val
l.set_ydata(((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000)
fig.canvas.draw_idle()
sROC.on_changed(update)
sx_shift.on_changed(update)
sy_shift.on_changed(update)
plt.show()
注意：您可以从此处获取轮廓仪数据示例：https://1drv.ms/u/s!AgMbHdCbxV3LatCuPSJpioDD0w0
如果您不能轻易截断数据，则需要对所有数据进行建模。我没有尝试使用您的数据（它似乎没有发布——也许会发布？）但数据似乎可以建模为
Parabola * Rectangle_with_soft_boundaries
使用抛物线模型给出您正在寻找的曲率半径，使用矩形模型给出 "truncated background"。
您可能会发现 lmfit (https://github.com/lmfit/lmfit-py) 对此很有用。它为 curve-fitting 提供了相当 high-level 的接口，并且有许多 built-in 模型可以组合。
例如，您可以这样做：
from lmfit.models import RectangleModel
#Define the fitting function
def func(x, roc, x_shift, y_shift):
return ((roc*1000)-((roc*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
# turn model function into model,
# multiply by rectangle with error-function for sides:
model = Model(func) * RectangleModel(prefix='r_', form='erf')
params = model.make_params(r_amplitude=-600, r_center1=150,
r_center2=375, r_sigma1=10, r_sigma2=10,
roc=45, x_shift=200, y_shift=1)
# you can apply bounds to any of the parameters here:
params['roc'].min = 1
params['roc'].max = 100
params['r_amplitude'].max = 0 # make sure rectangle points down
params['r_center1'].min = 50
params['r_center1'].max = 250
# then perform fit
result = model.fit(y, params, x=x)
# print report of statistics, parameter values
print(result.fit_report())
# plot results
plt.plot(x, y)
plt.plot(x, result.best_fit)
plt.show()
当然，您仍然需要了解一些关于矩形背景在何处打开和关闭的信息——我怀疑这是来自您可能确实了解（或可以弄清楚！）的某个光圈。
感谢大家的评论和回答
我可以想出另一种方法来解决这个问题（解释可能很长，所以我更喜欢 post 作为答案而不是简单地编辑问题）
请耐心等待我 只编码了 3 个月，我很高兴得到您的支持：
所以，我没有使用滑块来调整我的拟合曲线，而是使用点击来确定我想要执行拟合的范围（我最后会有 2 个问题)
import tkinter.filedialog as tk
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
fig, ax = plt.subplots()
def openfile(): #define function to open file and extract x and y
opennm = tk.askopenfile()
f = open(opennm.name,"r")
data=np.genfromtxt(f,
skip_header=6,
names=True,
dtype=None,
delimiter=',')
x=[]
y=[]
for i in range(0,len(data)-1):
x.append(data[i][0])
y.append(data[i][1])
return x,y
x,y=openfile()
k,= plt.plot(x,y)
coords = [] #I define the coordinates from where to where i want to run the fitting
def onclick(event): #this function to point and click on graph and define fitting area
global ix, iy
ix, iy = event.xdata , event.ydata
print ('x = %d, y = %d'%(
ix, iy))
global coords
coords.append((ix, iy))
if len(coords) == 2:
fig.canvas.mpl_disconnect(cid)
global l
for i in coords:
l,=ax.plot (*i,'ro' )
fig.canvas.draw_idle()
global m
def find_index(array,value): #this functions is to find index of my 2 point-and-click data
global limit
limit=(np.abs(array-value)).argmin()
return limit
limit1=find_index(x,coords[0][0]) #I extract index of Xstart
limit2=find_index(x,coords[1][0]) #I extraxt index of Xend
fitx=x[limit1:limit2] #From initial x data I extract the range where the fitting will be done
fity=y[limit1:limit2] #From initial y data I extract the range where the fitting will be done
def func(x,ROC,x_shift,y_shift): #This si the function I want to fit
return ((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
popt, pcov = curve_fit(func, fitx, fity, p0=[45,250,0.8],bounds=((0.4,100,0),(500,300,2))) #I run the fitting
m, = plt.plot(x,func(x,*popt),'r--')
print('ROC= ',round(popt[0],0),'mm| X Shift= ', round(popt[1],0),'um| Y Shift= ',round(popt[2],3),'um')
return coords, l,m
cid = fig.canvas.mpl_connect('button_press_event', onclick)
plt.show()
问题 1:
while 运行: 我 select 拟合范围的下限，在 select 上限之前，我收到此消息：
limit2=find_index(x,coords[1][0]) IndexError: list index out of range
这很自然，因为我还没有 select 上限...如何避免这个烦人的错误消息？
问题二：
运行拟合后，图形auto-scales...如何固定y轴的上限并保持自动缩放的下限？
我浏览了我在互联网上找到的所有建议解决方案（plt.gca().set_ylim(top=200)...等）但每次下限都固定为 0自动
除了这两个烦人的错误，程序运行正常
再一次，我 post link 从我的轮廓仪中获取测试数据：
我从轮廓仪收集数据并尝试将收集的数据拟合到理想曲线。目的是获取曲率半径。
我编写的代码将在所有原始数据上拟合一条曲线。
适合整个范围
问题是想指定一个范围来做拟合。即：所有轮廓仪数据都不适合拟合。在下面的例子中，我想适应从 x= 180 到 x = 380
在指定范围内拟合
当然，点击指定 x_start 和 x_end 会更有帮助，但我未能从 onclick
中提取数据奖励：如果我可以在更改拟合范围时显示 R²，那就太棒了！
感谢您的帮助
这是我想出的代码：
#importing modules
import tkinter.filedialog as tk
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
from matplotlib.widgets import Slider
#Graphic properties
fig, ax = plt.subplots()
plt.subplots_adjust(left=0.15, bottom=0.45)
#Search and open desired file:
def openfile():
opennm = tk.askopenfile()
f = open(opennm.name,"r")
data=np.genfromtxt(f,
skip_header=6,
names=True,
dtype=None,
delimiter=',')
x=[]
y=[]
for i in range(0,len(data)-1):
x.append(data[i][0])
y.append(data[i][1])
return x,y
#after opening the file, x and y data are plotted
x,y=openfile()
k,= plt.plot(x,y)
#Define the fitting function
def func(x,ROC,x_shift,y_shift):
return ((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
popt, pcov = curve_fit(func, x, y, p0=[45,250,0.8],bounds=((1,100,0),(100,300,2))) #, bounds=((-np.inf,10**-8,-np.inf,-np.inf),(np.inf,3.,np.inf,np.inf))
#The fitted curve will be overlaid on raw data
l, = plt.plot(x,func(x,*popt),'r--')
list_of_para=['ROC','X shift','Y shift']
list_of_units=['mm','a.u.','um']
for i in range (0,3):
print (list_of_para[i],'= ',"%.2f" %popt[i],' ', list_of_units[i])
#definig initial values
ROC0 = popt[0]
x_shift0 = popt[1]
y_shift0 = popt[2]
#defining sliders
axcolor = 'lightgoldenrodyellow'
axROC= plt.axes([0.15, 0.05, 0.75, 0.02], facecolor=axcolor)
axx_shift= plt.axes([0.15, 0.15, 0.75, 0.02], facecolor=axcolor)
axy_shift= plt.axes([0.15, 0.25, 0.75, 0.02], facecolor=axcolor)
#define slider limites
sROC = Slider(axROC, 'ROC', popt[0]-np.absolute(popt[0]),popt[0]+np.absolute(popt[0]), valinit=ROC0)
sx_shift = Slider(axx_shift, 'X Shift', popt[1]-np.absolute(popt[1]), popt[1]+np.absolute(popt[1]), valinit=x_shift0)
sy_shift = Slider(axy_shift, 'Y shift', popt[2]-np.absolute(popt[2]), popt[2]+np.absolute(popt[2]), valinit=y_shift0)
#define slider update values
def update(val):
ROC = sROC.val
x_shift = sx_shift.val
y_shift = sy_shift.val
l.set_ydata(((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000)
fig.canvas.draw_idle()
sROC.on_changed(update)
sx_shift.on_changed(update)
sy_shift.on_changed(update)
plt.show()
注意：您可以从此处获取轮廓仪数据示例：https://1drv.ms/u/s!AgMbHdCbxV3LatCuPSJpioDD0w0
如果您不能轻易截断数据，则需要对所有数据进行建模。我没有尝试使用您的数据（它似乎没有发布——也许会发布？）但数据似乎可以建模为
Parabola * Rectangle_with_soft_boundaries
使用抛物线模型给出您正在寻找的曲率半径，使用矩形模型给出 "truncated background"。
您可能会发现 lmfit (https://github.com/lmfit/lmfit-py) 对此很有用。它为 curve-fitting 提供了相当 high-level 的接口，并且有许多 built-in 模型可以组合。
例如，您可以这样做：
from lmfit.models import RectangleModel
#Define the fitting function
def func(x, roc, x_shift, y_shift):
return ((roc*1000)-((roc*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
# turn model function into model,
# multiply by rectangle with error-function for sides:
model = Model(func) * RectangleModel(prefix='r_', form='erf')
params = model.make_params(r_amplitude=-600, r_center1=150,
r_center2=375, r_sigma1=10, r_sigma2=10,
roc=45, x_shift=200, y_shift=1)
# you can apply bounds to any of the parameters here:
params['roc'].min = 1
params['roc'].max = 100
params['r_amplitude'].max = 0 # make sure rectangle points down
params['r_center1'].min = 50
params['r_center1'].max = 250
# then perform fit
result = model.fit(y, params, x=x)
# print report of statistics, parameter values
print(result.fit_report())
# plot results
plt.plot(x, y)
plt.plot(x, result.best_fit)
plt.show()
当然，您仍然需要了解一些关于矩形背景在何处打开和关闭的信息——我怀疑这是来自您可能确实了解（或可以弄清楚！）的某个光圈。
感谢大家的评论和回答
我可以想出另一种方法来解决这个问题（解释可能很长，所以我更喜欢 post 作为答案而不是简单地编辑问题）
请耐心等待我 只编码了 3 个月，我很高兴得到您的支持：
所以，我没有使用滑块来调整我的拟合曲线，而是使用点击来确定我想要执行拟合的范围（我最后会有 2 个问题)
import tkinter.filedialog as tk
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
fig, ax = plt.subplots()
def openfile(): #define function to open file and extract x and y
opennm = tk.askopenfile()
f = open(opennm.name,"r")
data=np.genfromtxt(f,
skip_header=6,
names=True,
dtype=None,
delimiter=',')
x=[]
y=[]
for i in range(0,len(data)-1):
x.append(data[i][0])
y.append(data[i][1])
return x,y
x,y=openfile()
k,= plt.plot(x,y)
coords = [] #I define the coordinates from where to where i want to run the fitting
def onclick(event): #this function to point and click on graph and define fitting area
global ix, iy
ix, iy = event.xdata , event.ydata
print ('x = %d, y = %d'%(
ix, iy))
global coords
coords.append((ix, iy))
if len(coords) == 2:
fig.canvas.mpl_disconnect(cid)
global l
for i in coords:
l,=ax.plot (*i,'ro' )
fig.canvas.draw_idle()
global m
def find_index(array,value): #this functions is to find index of my 2 point-and-click data
global limit
limit=(np.abs(array-value)).argmin()
return limit
limit1=find_index(x,coords[0][0]) #I extract index of Xstart
limit2=find_index(x,coords[1][0]) #I extraxt index of Xend
fitx=x[limit1:limit2] #From initial x data I extract the range where the fitting will be done
fity=y[limit1:limit2] #From initial y data I extract the range where the fitting will be done
def func(x,ROC,x_shift,y_shift): #This si the function I want to fit
return ((ROC*1000)-((ROC*1000)**2-(x-x_shift)**2)**0.5-y_shift)*1000
popt, pcov = curve_fit(func, fitx, fity, p0=[45,250,0.8],bounds=((0.4,100,0),(500,300,2))) #I run the fitting
m, = plt.plot(x,func(x,*popt),'r--')
print('ROC= ',round(popt[0],0),'mm| X Shift= ', round(popt[1],0),'um| Y Shift= ',round(popt[2],3),'um')
return coords, l,m
cid = fig.canvas.mpl_connect('button_press_event', onclick)
plt.show()
问题 1: while 运行: 我 select 拟合范围的下限，在 select 上限之前，我收到此消息：
limit2=find_index(x,coords[1][0]) IndexError: list index out of range
这很自然，因为我还没有 select 上限...如何避免这个烦人的错误消息？
问题二： 运行拟合后，图形auto-scales...如何固定y轴的上限并保持自动缩放的下限？
我浏览了我在互联网上找到的所有建议解决方案（plt.gca().set_ylim(top=200)...等）但每次下限都固定为 0自动
除了这两个烦人的错误，程序运行正常
再一次，我 post link 从我的轮廓仪中获取测试数据：