Learning Tomography: Weighted Back-Projection

In my last posts, we were able to compute a 2D-reconstruction from a series of 1D-projections ... but the result appears blurred. How can we fix this?

1- Why is my reconstruction blurry?

As shown in Fig.1, each back-projection of Fig.1A (used to calculate the temporary 2D-reconstruction of Fig.1B) corresponds in Fourier Space to a central line resembling to spokes of a bicycle wheel  (Fig.1C). The mathematics behind this phenomenon is called the central slice theorem [Wikipedia] stating that a 2D-reconstruction  can be calculated in Fourier Space from FFTs of each 1D-projection correctly oriented and inserted as a central slice.

Fig.1: Central slices in Fourier Space. A) Ten back-projections calculated from a sinogram. B) Temporary 2D-reconstruction calculated from series (A). C) FFT of image (B). Each back-projection appears as a central line (a 1D-slice) rotated according to its projection direction.

Note: There is a family of reconstruction techniques using the central slice theorem called Direct Fourier Reconstruction (a subject for a new post? Yes, I do it, it's available  [here]).

We now know that a 2D-reconstruction appears as a bicycle wheel in Fourier space!?! What about the blurring artifact?
In the central area of the FFT image (called a power spectrum) corresponding to the low frequencies, the various spokes are close to each other and thus, there is no "empty space". However, at the periphery (the high frequencies area), the spokes are distant. In other words, we have much more missing data at the periphery than at the center.

Fig.2: xxx TODO xxx

In conclusion, we have more data in low frequencies than in high frequencies. This uneven distribution between the low and high frequencies is responsible for the blurring effect observed in the simple back-projection algorithm.

2- Filtering the projections.

We saw previously that there is an imbalance between the low and high frequencies. We must design a filter to counterbalance this effect. This filter is called a ramp filter and is characterized by this typical profile (Fig.3).

Fig.3: Filter Profile

To build such a filter, create a 256x256 32-bit image entitled "Filter" with the command File > New > Image... and fill it using the Process > Math > Macro... with the formula  'v = d'.

Fig.4: Formula allowing the creation of the filter

The resulting image is presented in Fig. 5.

Fig. 5: Ramp filter

Now, we have to filter the sinogram in Fourier Space, click on the sinogram image (to make it active) and go to the Process > FFT > Custom Filter..., choose in the menu "Filter", the sinogram now appears more contrasted (the edges are brighter as shown in Fig. 6B).

Fig.6: Sinogram before and after filtration.

Finally, run the IJ script backProj.js [Link] described in my previous posts (use the improved version, because we don't need intermediate back-projections).

Fig.7: 2D-reconstruction calculated from sinogram of Fig.3B. A) Before  and (B) after normalization. Lena is now correctly reconstructed.

Note: If you have some problem to normalize the rec2D, select a central area with the rectangle or ellipse tool (Tools #1 and 2) and run again the Process > Enhance Contrast (only check 'Normalize' and 0.0% of 'Saturated Pixels'). This should be better.

3- Other filters ...

::: Work in Progress :::

The filtering is done in Fourier Space, however there are equivalent filters in real space (or spatial domain).

3-1- Laplacian filters

Here are 3x3 and 5x5 Laplace filters.

        
-1 -1 -1
-1  8 -1
-1 -1 -1
             
-1 -1 -1 -1 -1
-1 -1 -1 -1 -1
-1 -1 24 -1 -1
-1 -1 -1 -1 -1
-1 -1 -1 -1 -1

 Note: Laplacian 5x5 is available by default in Process > Filters > Convolve...

3-2-  'Mexican Hat' family

Mexican hat 7x7

 
 0  0 -1 -1 -1  0  0
 0 -1 -3 -3 -3 -1  0
-1 -3  0  7  0 -3 -1
-1 -3  7 24  7 -3 -1
-1 -3  0  7  0 -3 -1
 0 -1 -3 -3 -3 -1  0
 0  0 -1 -1 -1  0  0

Mexican hat 9x9

 0  0  0 -1 -1 -1  0  0  0
 0 -1 -1 -3 -3 -3 -1 -1  0
 0 -1 -3 -3 -1 -3 -3 -1  0
-1 -3 -3  6 13  6 -3 -3 -1
-1 -3 -1 13 24 13 -1 -3 -1
-1 -3 -3  6 13  6 -3 -3 -1
 0 -1 -3 -3 -1 -3 -3 -1  0
 0 -1 -1 -3 -3 -3 -1 -1  0
 0  0  0 -1 -1 -1  0  0  0
 

Mexican hat 13x13

 0  0  0  0  0 -1 -1 -1  0  0  0  0  0
 0  0  0 -1 -1 -2 -2 -2 -1 -1  0  0  0
 0  0 -2 -2 -3 -3 -4 -3 -3 -2 -2  0  0
 0 -1 -2 -3 -3 -3 -2 -3 -3 -3 -2 -1  0
 0 -1 -3 -3 -1  4  6  4 -1 -3 -3 -1  0
-1 -2 -3 -3  4 14 19 14  4 -3 -3 -2 -1
-1 -2 -4 -2  6 19 24 19  6 -2 -4 -2 -1
-1 -2 -3 -3  4 14 19 14  4 -3 -3 -2 -1
 0 -1 -3 -3 -1  4  6  4 -1 -3 -3 -1  0
 0 -1 -2 -3 -3 -3 -2 -3 -3 -3 -2 -1  0
 0  0 -2 -2 -3 -3 -4 -3 -3 -2 -2  0  0
 0  0  0 -1 -1 -2 -2 -2 -1 -1  0  0  0
 0  0  0  0  0 -1 -1 -1  0  0  0  0  0

Fig.8: 2D-reconstruction obtained from a sinogram filtered with a Mexican Hat 7x7.

Hope that helps!

Other crazybiocomputing posts

Further readings are available in ...

• Learning Tomography Series  [Link]
• Image Processing TOC [Link]