确定绘图宽度的 Pythonic 方法

Question

我有以下情节。

我正在寻找一种 Python 方法来查找 X 和 Y 轮廓的宽度（如您所见，宽度应该约为 0.002）。

我已经考虑过的非 Python 方法包括向前和向后循环遍历 X 和 Y 配置文件列表以确定大于某个值（可能是最大值除以 3）的第一个元素。

Answer 1

免责声明：“pythonic”这个词对我来说没有太大意义；所以这只是两个解决方案，如果他们需要一个名字，我建议“erinaceidaeic”或“axolotlable”。

这个问题肯定有一个物理组成部分，即定义“宽度”的含义。我可以想象确定 Hermite-Gaussian 模式的 w 参数会很有趣。但是，出于这个问题作为编程问题的目的，假设您想找到半高全宽 (FWHM)。

这里有两个 FWHM 函数。

import numpy as np

#########
# FWHM function
#########

def cFWHM(x,y):
    """ returns coarse full width at half maximum, and the two
        xcoordinates of the first and last values above the half maximum """
    where, = np.where(y >= y.max()/2.)
    maxi = x[where[-1]]
    mini = x[where[0]]
    return maxi-mini, mini, maxi

def fFWHM(x,y):
    """ returns interpolated full width at half maximum, and the two
        xcoordinates at the (interpolated) half maximum """
    def find_roots(x,y):
        s = np.abs(np.diff(np.sign(y))).astype(bool)
        return x[:-1][s] + np.diff(x)[s]/(np.abs(y[1:][s]/y[:-1][s])+1)

    z = find_roots(x,y-y.max()/2)
    return z.max()-z.min(), z.min(), z.max()

如问题中所述，第一种方法是沿坐标轴找到最小和最大坐标，其中y 的值大于或等于数据中最大 y 值的一半。它们之间的区别是宽度。对于足够密集的点，这给出了合理的结果。

如果需要更高的精度，可以通过在数据点之间进行插值来找到y-ymax/2 的零点。 （解决方案取自my answer to this question）。

完整示例：

import matplotlib.pyplot as plt

#########
# Generate some data
#########

def H(n, x):
    # Get the nth hermite polynomial, evaluated at x
    coeff = np.zeros(n)
    coeff[-1] = 1
    return np.polynomial.hermite.hermval(x, coeff)

def E(x,y,w,l,m, E0=1):
    # get the hermite-gaussian TEM_l,m mode in the focus (z=0)
    return E0*H(l,np.sqrt(2)*x/w)*H(m,np.sqrt(2)*y/w)*np.exp(-(x**2+y**2)/w**2)

g = np.linspace(-4.5e-3,4.5e-3,901)
X,Y = np.meshgrid(g,g)
f = E(X,Y,.9e-3, 5,7)**2

# Intensity profiles along x and y direction
Int_x = np.sum(f, axis=0)
Int_y = np.sum(f, axis=1)

#########
# Plotting
#########

fig = plt.figure(figsize=(8,4.5))
ax = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,4)


dx = np.diff(X[0])[0]; dy = np.diff(Y[:,0])[0]
extent = [X[0,0]-dx/2., X[0,-1]+dx/2., Y[0,0]-dy/2., Y[-1,0]+dy/2.]
ax.imshow(f, extent=extent)


ax2.plot(g,Int_x)
ax3.plot(g,Int_y)

width, x1, x2 = cFWHM(g,Int_x)      # compare to fFWHM(g,Int_x)
height, y1, y2 = cFWHM(g,Int_y)

ax2.plot([x1, x2],[Int_x.max()/2.]*2, color="crimson", marker="o")
ax3.plot([y1, y2],[Int_y.max()/2.]*2, color="crimson", marker="o")

annkw = dict(xytext=(0,10), 
             textcoords="offset pixels", color="crimson", ha="center")
ax2.annotate(width, xy=(x1+width/2, Int_x.max()/2.), **annkw)
ax3.annotate(height, xy=(y1+height/2, Int_y.max()/2.), **annkw)

plt.tight_layout()
plt.show()

这是两个函数之间比较的视觉效果。使用第一个而不是第二个函数的最大误差是后续数据值之间的差异。在这种情况下，这可能是相机分辨率。 （但请注意，真正的误差当然需要在确定最大值时考虑误差，因此可能会大得多。）

Answer 2

您可以对配置文件应用硬阈值。 您可以通过查看均值+一些标准差计算的背景强度来得出阈值（例如，您可以对图像的 10% 外边界进行采样）