A Pythonic method that determines the width of a drawing

I have the following episode.

I’m looking for a pythonic way to find the width of the X and Y sections (as you can see, the width should be around 0.002).

The non-pythonic approach I’ve considered involves looping forward and backward through the X and Y profile lists to determine the first element greater than a certain value (perhaps the maximum divided by 3).

Disclaimer: The word “pythonic” doesn’t mean much to me; So these are just two solutions, and if they need a name, I recommend using “erinaceidaeic” or “axolotlable”.

There must be a physical component to this problem, which is the meaning of defining “width”. I can imagine it would be interesting to determine the w parameter of the Hermite-Gaussian schema. However, for the purposes of this problem as a programming problem, suppose you want to find half-peak full width (FWHM).

There are two FWHM functions here.

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

As described in the question, the first method finds the minimum and maximum coordinates along the axis, where the value of y is greater than or equal to half of the maximum y value in the data. The difference between them is in the width. For sufficiently dense points, this gives reasonable results.

If you need more precision, you can find the zero point of y-ymax/2 by interpolating between data points. (Solution taken from my answer to this question).

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

Here is a comparison diagram between two functions. The largest error of using the first function instead of the second is the difference between subsequent data values. In this case, this could be the camera resolution. (Note, however, that the true error certainly needs to be considered when determining the maximum value, so it can be much larger.) )