高斯拟合计算成像分辨率

这里简单展示如何对高斯分布的曲线进行拟合以计算成像分辨率。

由于中心极限定理,我们采集到的数据经常呈现正态分布(也叫高斯分布)。成像数据中,有时我们需要评估成像的分辨率。一个比较直接简单的做法,是找到单个光斑进行划线测量,或者是对线性结构(如细胞微管)的进行垂直划线测量。

中心极限定理说的啥?它说样本的平均值约等于总体的平均值。然后实际测量过程中我们通常观测到的不是单分子单颗粒这样的事件,所以得到的第一手数据其实是若干个事件的平均值,这些平均值的分布,必然呈现正态分布。

Pasted-image-20240905145810.png-e493ca849f.png

这里以在 ImageJ 中划线测量得到的 profile 数据为例,保存的 csv 文件,可以使用python中的 pandas 模块打开,为了方便,直接打包了一个函数,可以同时显示和读取 profile。

ImageJ划线测量详见:

图像划线测量

本文介绍对图像进行划线测量的操作 通过前面我们已经知道有各种图像选区的工具,针对不同的感兴趣对象可以使用不同的工具。上篇中使用了 freehand 来对细胞进行选择,这里则进一步展示使用 line 来对凝胶电泳成像的结果进行划线测量,获取其 profile。 通过拍摄得到的原始图像是黑底亮带,然而在...

import pandas as pd
def show_profile(fp):
'''
fp: imagej中划线测量保存的 data.csv,包含 distance 和 gray_value两列
'''
c = pd.read_csv(fp)
c.columns = ['distance', 'value']
c.plot(x='distance', y='value', figsize=(4,3))
plt.show()
return c
Pasted-image-20240905150547.png-a204a31b58.png

读取一个 profile,直接作图,方便高斯拟合的时候,提供比较接近的初始估计值。

import numpy as np
# 定义高斯函数
def gaussian(x, a, x0, sigma, b):
'''
x: 输入变量
a: 参数,缩放量
x0:高斯峰所在位置,相当于高斯分布平均值
sigma:高斯分布一个标准差
b: 参数,平移量
'''
return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))+b
initial_guess = [1000, 0.5, 0.2, 10]
x_fit, y_fit = single_gauss_fit(df, initial_guess)

这里我们使用的高斯函数,具有四个变量,可以大致估计下,当前要拟合的数据部分,大致的参数。比较接近的参数,有利于拟合算法快速收敛。为了方便,我们这里封装了一个单个高斯峰拟合的函数,执行后会直接显示拟合的效果,并且打印各个参数,以及参数的误差,还有R方值,最后把拟合的数据返回回来。

Pasted-image-20240905151027.png-a59e44f8df.png

封装好的拟合函数代码如下:

from scipy.optimize import curve_fit
def single_gauss_fit(a:pd.DataFrame, initial_guess:list):
'''
a: 包含 distance 和 value 两列,其中distance以μm为单位
initial_guess: 对高斯分布各个参数的估计值
返回拟合后的曲线的xy值
'''
# 进行拟合
x = a.distance.values
y = a.value.values
popt, pcov = curve_fit(gaussian, x, y, p0=initial_guess)
x2 = np.linspace(x.min(),x.max(), 100) # 根据x范围均匀插值
y2 = gaussian(x2, *popt)
plt.figure(figsize=(4,3))
plt.scatter(x, y, c='gray')
plt.plot( x2, y2, c='red')
plt.show()
# 计算各个参数的拟合误差
perr = np.sqrt(np.diag(pcov))
# 计算拟合的R方值
mean = np.mean(y)
ss_tot = np.sum((y - mean) ** 2)
y2 = gaussian(x, *popt)
ss_res = np.sum((y - y2) ** 2)
r_squared = 1 - (ss_res / ss_tot)
print(f'R_squared:{r_squared}')
# 计算半峰高全宽
# FWHM=2*sqrt(2ln2)*σ=2.355σ
FWHM = 2 * np.sqrt(2 * np.log(2)) * abs(popt[2])
# 打印拟合参数和误差
print("参数:", popt)
print("参数误差:", perr)
print(f"Peak location: {popt[1]*1000:.1f} nm\nSigma: {abs(popt[2])*1000:.1f} nm\nFWHM(分辨率): {FWHM*1000:.1f} nm")
return x2, y2