Full width at half maximum is the width of a peak measured halfway up it. In image analysis it is the standard way to describe how sharp a point source is — a star in astronomy, a bead in microscopy, a spot in a beam profile. This guide covers three ways to measure it in Python, from a quick interpolation to a proper two-dimensional Gaussian fit.
What it means#
Take a peak, find its maximum, halve that value, and measure how wide the peak is at that height. That width is the FWHM. It is preferred over other width measures because it is easy to read off real data and does not depend on where you decide the peak “ends” — a decision that noise makes arbitrary.
For a Gaussian there is an exact relationship between FWHM and the standard deviation:
FWHM = 2 * sqrt(2 * ln(2)) * sigma
= 2.3548 * sigma
Setting up#
python -m pip install numpy scipy matplotlib pillow
import numpy as np
from PIL import Image
image = np.array(Image.open("star.png").convert("L"), dtype=float)
print(image.shape, image.min(), image.max())
Converting to float matters. Integer arrays wrap around on subtraction, which silently corrupts background removal.
Method 1: interpolate the crossings#
Fast, no fitting, and good enough when the signal is clean:
import numpy as np
def fwhm_from_profile(x, y):
"""FWHM of a single peak, by linear interpolation at half maximum."""
y = np.asarray(y, dtype=float)
x = np.asarray(x, dtype=float)
baseline = np.median(y)
y = y - baseline
peak_index = int(np.argmax(y))
half = y[peak_index] / 2.0
# Walk left until the profile drops below half
left = peak_index
while left > 0 and y[left] > half:
left -= 1
right = peak_index
while right < len(y) - 1 and y[right] > half:
right += 1
if left == peak_index or right == peak_index:
raise ValueError("Peak does not fall below half maximum within the profile")
x_left = np.interp(half, [y[left], y[left + 1]], [x[left], x[left + 1]])
x_right = np.interp(half, [y[right], y[right - 1]], [x[right], x[right - 1]])
return abs(x_right - x_left)
The interpolation is what makes this usable. Without it you would be limited to whole-pixel resolution, and typical FWHM values are only a few pixels wide.
row = image[image.shape[0] // 2, :]
x = np.arange(row.size)
print("FWHM:", round(fwhm_from_profile(x, row), 3), "pixels")
Method 2: fit a Gaussian#
Noise makes the interpolation method jump around. Fitting uses every point, so it is far steadier:
import numpy as np
from scipy.optimize import curve_fit
def gaussian(x, amplitude, centre, sigma, offset):
return amplitude * np.exp(-((x - centre) ** 2) / (2 * sigma ** 2)) + offset
def fwhm_by_fit(x, y):
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
guess = [
y.max() - np.median(y), # amplitude
x[int(np.argmax(y))], # centre
len(x) / 10.0, # sigma
np.median(y), # offset
]
params, covariance = curve_fit(gaussian, x, y, p0=guess, maxfev=10000)
amplitude, centre, sigma, offset = params
fwhm = 2 * np.sqrt(2 * np.log(2)) * abs(sigma)
sigma_error = np.sqrt(np.diag(covariance))[2]
return fwhm, 2.3548 * sigma_error, centre
fwhm, error, centre = fwhm_by_fit(x, row)
print(f"FWHM = {fwhm:.3f} +/- {error:.3f} pixels, centred at {centre:.2f}")
The covariance matrix gives you an uncertainty for free, which the interpolation method cannot provide. Reporting a width without an uncertainty is rarely useful in real work.
Method 3: two-dimensional fitting#
Taking a single row assumes the peak is centred on that row and symmetric. Fitting the whole spot is more honest:
import numpy as np
from scipy.optimize import curve_fit
def gaussian_2d(coords, amplitude, x0, y0, sigma_x, sigma_y, offset):
x, y = coords
value = amplitude * np.exp(
-(((x - x0) ** 2) / (2 * sigma_x ** 2) + ((y - y0) ** 2) / (2 * sigma_y ** 2))
) + offset
return value.ravel()
def fwhm_2d(patch):
ny, nx = patch.shape
y, x = np.mgrid[0:ny, 0:nx]
guess = [
patch.max() - np.median(patch),
nx / 2.0, ny / 2.0,
nx / 6.0, ny / 6.0,
np.median(patch),
]
params, _ = curve_fit(
gaussian_2d, (x, y), patch.ravel(), p0=guess, maxfev=20000
)
_, x0, y0, sigma_x, sigma_y, _ = params
factor = 2 * np.sqrt(2 * np.log(2))
return factor * abs(sigma_x), factor * abs(sigma_y), (x0, y0)
// Cut a small window around the brightest pixel first
cy, cx = np.unravel_index(np.argmax(image), image.shape)
half = 15
patch = image[cy - half:cy + half + 1, cx - half:cx + half + 1]
fx, fy, centre = fwhm_2d(patch)
print(f"FWHM x = {fx:.2f}, y = {fy:.2f} pixels")
Cropping around the peak matters. Fitting a whole image containing several sources will not converge on any of them.
Converting to physical units#
pixel_scale = 0.31 # arcseconds per pixel, or microns per pixel
print(f"FWHM = {fwhm * pixel_scale:.3f} arcsec")
The pixel value is only meaningful with a scale attached. In astronomy that is the plate scale; in microscopy it is the calibrated pixel size, which you should measure rather than assume.
Checking the result#
import matplotlib.pyplot as plt
fit_y = gaussian(x, *params)
plt.plot(x, row, "o", markersize=3, label="data")
plt.plot(x, fit_y, "-", label="fit")
plt.axhline(offset + amplitude / 2, linestyle="--", label="half maximum")
plt.legend()
plt.xlabel("pixel")
plt.ylabel("intensity")
plt.show()
Always plot the fit over the data. A number from curve_fit that looks plausible can still come from a fit that missed the peak entirely, and one glance at the plot tells you.
rng = np.random.default_rng(42)
x = np.arange(100)
truth = gaussian(x, 1000, 50, 4.0, 100)
noisy = truth + rng.normal(0, 20, x.size)
print("expected:", round(2.3548 * 4.0, 3))
print("measured:", round(fwhm_by_fit(x, noisy)[0], 3))
Questions people ask#
Which method should I use?
Interpolation for a quick look at clean data. A one-dimensional fit for a proper measurement of a symmetric peak. Two-dimensional fitting when the shape is elongated or you need the centre position accurately.
What if my peak is not Gaussian?
Fit the shape that matches. Astronomical point sources are often better described by a Moffat profile; spectroscopic lines are sometimes Lorentzian or Voigt. The conversion from a shape parameter to FWHM differs for each.
Why does the same star give different values on different frames?
Genuine variation — seeing, focus, tracking — plus measurement noise. Measure several sources per frame and report the median, not a single value.
Does saturation affect it?
Badly. A flat-topped peak has no meaningful maximum, so half maximum is wrong and the fit is biased. Check whether the peak reaches the detector’s maximum value and exclude it if so.
Where to go next#
- Data visualization — plotting profiles and fits properly.
- Python and data science — NumPy and SciPy foundations.
- Exploring np.stack() — combining image frames.
