Beam hardening in micro CT

One of the first things I noticed when I started working with micro CT scans of printed specimens was that a solid region looked darker in the middle than at its edges. Same material, different grey values. That is beam hardening, and it is one of the first things I correct now, because everything after it (segmentation, porosity, mapping grey values to stiffness) treats the grey values as real physics.

Why it happens

An X-ray tube does not send out photons of one energy. It sends a whole spectrum. The low energy photons get absorbed much more easily than the high energy ones, so the deeper the beam goes into the part, the more of the soft photons are gone. What is left has a higher average energy. The beam gets "harder" on its way through.

If the beam had only one energy, attenuation would follow the Beer Lambert law:

\[ I = I_0 \exp\left(-\int_L \mu(s)\,ds\right) \]

The reconstruction works with the projection, which is the line integral of \(\mu\):

\[ p = -\ln\frac{I}{I_0} = \int_L \mu(s)\,ds \]

For a uniform material of thickness \(t\) this is just \(p = \mu t\), a straight line. Filtered back projection assumes exactly that.

With a real spectrum we have to integrate over energy instead:

\[ I = \int_0^{E_{\max}} I_0(E)\,\exp\left(-\int_L \mu(s,E)\,ds\right) dE \]

Now \(\mu\) depends on energy and the spectrum keeps shifting, so every extra millimetre attenuates a bit less than the one before. The measured \(p\) bends below the straight line. The reconstruction does not know that, so it underestimates attenuation along the longest paths, which go through the centre. That gives the cup shape, and also the dark streaks you see between two dense features.

For my work this matters directly. If the centre reads lower than the rim, one global Otsu threshold cuts the part differently in the middle than at the edge, and the pore statistics come out wrong. So I correct it on the projections, before reconstruction.

Straightening the curve

The usual software correction is linearization. You scan a step wedge of the same material, where you know each thickness, and measure the hardened projection for every step. Then you fit a polynomial that moves each measured value back onto the straight line it should be on, and apply that polynomial to every ray in the sinogram. Zero thickness has to give zero projection, so there is no constant term:

\[ p_{\text{corr}} = c_1 p + c_2 p^2 + c_3 p^3 \]

This is the logic in Python:

import numpy as np

# Step wedge of the SAME material.
# t  : known thickness of each step [mm]
# I  : measured intensity behind each step
# I0 : intensity in air (flat field)
t = np.array([2, 4, 6, 8, 10, 12], dtype=float)
p = -np.log(I / I0)                  # measured, hardened projections

# With one energy we would get p = mu * t.
# Take the slope from the thinnest steps, where hardening is smallest.
mu = np.polyfit(t[:2], p[:2], 1)[0]
p_ideal = mu * t                     # the straight line we want back

# Fit measured -> ideal with a cubic through the origin.
V = np.column_stack([p, p**2, p**3]) # no constant column
c, *_ = np.linalg.lstsq(V, p_ideal, rcond=None)

def linearize(p_raw):
    """Straighten a hardened sinogram before reconstruction."""
    return c[0] * p_raw + c[1] * p_raw**2 + c[2] * p_raw**3

sino = -np.log(proj / flat)
sino_corr = linearize(sino)          # then reconstruct as usual

To convince myself it works, I simulated a beam with three energies. Without correction, the attenuation per millimetre dropped from 0.100 to 0.084 as the thickness went from 2 to 12 mm. After linearization it stayed at 0.092 for every step.

Where it stops working

This only works well when one material dominates. With very different materials in the same scan, metal inside polymer for example, you need iterative or dual energy methods.

A lot can also be done before the software. A thin aluminium or copper filter at the tube removes the softest photons before they reach the part, so the spectrum changes less inside it. In practice both are used together. The step wedge data is also a quick way to see how much hardening your setup has: the more the raw points bend away from the straight line, the more you need to correct.

References

R. A. Brooks, G. Di Chiro, Beam hardening in X-ray reconstructive tomography, Physics in Medicine and Biology 21(3), 390–398, 1976.
G. T. Herman, Correction for beam hardening in computed tomography, Physics in Medicine and Biology 24(1), 81–106, 1979.