将曲线拟合到图形上的特定范围

Question

我从轮廓仪收集数据并尝试将收集的数据拟合到理想曲线。目的是得到曲率半径。

我编写的代码将拟合所有原始数据的曲线。

适合整个范围

问题是要指定一个范围来进行拟合。即：所有轮廓仪数据都不适合拟合。在下面的例子中，我想从 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

Answer 1

如果您无法轻松截断数据，则需要对所有数据进行建模。我没有尝试过使用您的数据（似乎没有发布——也许会这样做？）但似乎数据可以建模为

Parabola * Rectangle_with_soft_boundaries

抛物线模型给出您正在寻找的曲率半径，矩形给出“截断背景”。

您可能会发现 lmfit (https://github.com/lmfit/lmfit-py) 对此很有用。它为曲线拟合提供了一个相当高级的接口，并具有许多可以组合的内置模型。

例如，您可以执行以下操作：

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

当然，你仍然需要知道一些关于矩形背景在哪里打开和关闭的信息——我怀疑这是来自一些你可能知道（或者可以弄清楚！）的光圈。

Answer 2

谢谢大家的cmets和回答

我可以想出另一种解决方法（解释可能很长，所以我更喜欢将帖子作为答案而不是简单地编辑问题）

请耐心等待，我正在编写代码仅 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： 运行时：我选择拟合范围的下限，在选择上限之前我收到此消息：

limit2=find_index(x,coords[1][0]) IndexError: list index out of range

这很自然，因为我还没有选择上限...如何避免这个烦人的错误信息？

问题 2： 运行拟合后，图形自动缩放...如何固定y轴的上限并保持自动缩放的下限？

我浏览了我在互联网上找到的所有建议的解决方案 (plt.gca().set_ylim(top=200)...等)，但每次下限都会自动固定为 0

除了这两个恼人的错误，程序运行正常

再一次，我将链接发布到我的轮廓仪中的测试数据：

https://1drv.ms/u/s!AgMbHdCbxV3LatCuPSJpioDD0w0