Filtered backprojection

Background

The most common reconstruction method is the so called filtered backprojection (FBP). The idea is simple, the basic reconstruction using only a backprojection is blurry, hence perform a filtering step. That’s it, but there is on key aspect. Filter the projection data, i.e. the sinogram, instead of the backprojection. Then backproject the filtered sinogram.

Implementation

Let us implement this method quickly. First, we need a fast Fourier transformation (FFT) implementation. elsa has one, but it’s a little more cumbersome than what is necessary here. Hence, we will use SciPy’s FFT. Import it with:

from scipy.fft import fft, ifft

Now we need to perform some padding to ensure the FFT works properly:

# Convert to numpy array
np_sinogram = np.array(sinogram)
sinogram_shape = np_sinogram.shape[0]

# Add padding
projection_size_padded = max(64, int(2 ** np.ceil(np.log2(2 * sinogram_shape))))
pad_width = ((0, projection_size_padded - sinogram_shape), (0, 0))
padded_sinogram = np.pad(np_sinogram, pad_width, mode="constant", constant_values=0)

Note that this step is to some extent optional. It would be sufficient to round to the next even number for the padding. However, the reconstruction will be worse. Try it for yourself, replace the line regarding padding with something like: int(np.ceil(sinogram_shape / 2)) * 2, and check the reconstruction.

The next step is to create the filter. This will determine which frequencies one wishes to suppress and by how much. The simplest filter, is the ideal Ramp filter. It will suppress high frequency noise, compared to lower frequency noise. Interestingly, from a mathematical point of view, this is the exact solution. So let us create that

def ramp_filter(size):
    n = np.concatenate(
        (
            # increasing range from 1 to size/2, and again down to 1, step size 2
            np.arange(1, size / 2 + 1, 2, dtype=int),
            np.arange(size / 2 - 1, 0, -2, dtype=int),
        )
    )
    f = np.zeros(size)
    f[0] = 0.25
    f[1::2] = -1 / (np.pi * n) ** 2

    # See "Principles of Computerized Tomographic Imaging" by Avinash C. Kak and Malcolm Slaney,
    # Chap 3. Equation 61, for a detailed description of these steps
    return 2 * np.real(fft(f))[:, np.newaxis]

# Ramp filter
fourier_filter = ramp_filter(projection_size_padded)

Interested readers are encouraged to read the given reference to see more details. These details would be out of scope for this guide. Now, we need to perform the convolution between the sinogram and the filter. In Fourier space, this is just a component wise multiplication, then apply the inverse Fourier transform:

projection = fft(padded_sinogram, axis=0) * fourier_filter
filtered_sinogram = np.real(ifft(projection, axis=0)[:sinogram_shape, :])

Now, finally we perform the backprojection with the filtered sinogram.

# Go back into elsa space
filtered_sinogram = elsa.DataContainer(filtered_sinogram, sino_descriptor)

# Do the backprojection:
filtered_backprojection = projector.applyAdjoint(filtered_sinogram)

plt.imshow(np.array(filtered_backprojection))
plt.show()

This is a good enough looking reconstruction. Different filters, can improve the results further, especially certain scenarios of noise might be handled better.

Flaws of the Filtered Backprojection

The FBP is the de facto standard in commercial X-ray CT scanners. This is due to its efficiency. The first algorithms developed for X-ray tomographic reconstruction were iterative reconstruction algorithms. However, due to the immense amount of computation necessary for iterative reconstructions to be successful, they were not practical. Hence, the efficiency of the FBP was key in the beginning.

However, the FBP breaks down in so-called low dose scenarios. These are, for example, scenarios where the number of projection angles is low. You can try this out, change the trajectory given above, to instead of 420 angles, only use, say, 60. The artifacts in the reconstruction take over quickly, and as you will see in the next chapters, iterative reconstruction algorithms can deal with these kinds of scenarios better.

Another low dose scenario is reducing the tube current, resulting in very noise projection images. Also for this scenario, iterative reconstruction algorithms, particular ones with regularization, typically perform better than FBP.

Low dose reconstruction is tremendously important, as a high dosage of X-ray radiation is dangerous and has health risks associated with it. Hence, in the medical field it is important to reduce the dosage to a necessary minimum.

Another aspect is that the FBP as described here is 2D only. There exists a 3D extension (FDK) to it, which works similarly to the FBP and is still quite efficient. However, iterative algorithms are agnostic of the dimensions. And with the rise of GPU programming, iterative algorithms have become practically more feasible. Much research in the field is focused in iterative reconstructions, however, as mentioned above, until today, commercial CT scanners still mostly rely on the FBP (FDK).