I have the following plot.
I am looking for a pythonic way to find the width of the X and Y Profiles (as you can see, the width should be about 0.002).
Non-pythonic methods I have already considered include looping through the X and Y Profile Lists both forwards and backwards to determine the first element that is greater than a certain value (maybe the maximum value divided by 3).
Disclaimer: The word "pythonic" does not make too much sense to me; so those are simply two solutions, if they need a name, I'd suggest "erinaceidaeic" or "axolotlable".
There is sure a physics component to this problem which would be to define what "width" means. I could imagine it to be interesting to determine the w parameter of the Hermite-Gaussian mode. However, for the purpose of this question as a programming question, let's just say you want to find the full width at half maximum (FWHM).
Here are two FWHM functions.
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()
The first would, as mentionned in the question, find the minimal and maximal coordinate along the axes where the value of y is larger or equal than half the maximal y value in the data. The difference between those is the width. This gives reasonable results for sufficiently dense points.
If a larger accuracy is needed, one may find the zeros of y-ymax/2 instead, by interpolating between the data points. (Solution taken from my answer to this question).
Complete example:
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()
And here is a visual of the comparisson between the two functions. The maximum error of using the first instead of the second function is the difference between subsequent data values. In this case this would probably be the camera resoluton. (Note however that the true error will of course need to take the error in determining the maximum value into account and will hence probably be much larger.)
you could apply a hard threshold on the profiles. you could derive the value of the threshold by looking at mean+some standard deviations calculated background intensity (you can sample the 10% outer border of the image for instance)
