The Filtered Backprojection

The Fourier inversion formula

Hoping for a readily usable expression, let us start with the 2Dinverse Fourier-transform:

$$ f(x,y)=\left (\frac{1}{2\pi} \right )^{2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \mathfrak{F}\left [ f \right ]\left ( \xi_{1} ,\xi_{2}\right )e^{i\mathbf{x} \boldsymbol{\xi }} d\xi_{1} d\xi_{2} = \left (\frac{1}{2\pi} \right )^{2}\int_{0}^{2\pi}\int_{-\infty}^{\infty} \mathfrak{F}\left [ f \right ]\left ( r\boldsymbol{\omega} )\left | r \right | e^{ir\mathbf{x} \boldsymbol{\omega }}dr d\vartheta$

Filling in the Fourier transform of the Radon transform regarding the t variable:
$$f(x,y) = \left (\frac{1}{2\pi} \right )^{2}\int_{0}^{2\pi}\int_{-\infty}^{\infty} \mathfrak{F}_{t}\left [ \mathfrak{R}f \right ]\left ( r\boldsymbol{\omega} )\left | r \right | e^{ir\mathbf{x} \boldsymbol{\omega }}dr d\vartheta$

The resulting expression gains its usefulness from the fact that the inverse Fourier transform is only carried out in 1D and the sampling points are evenly distribuited.

The filter term of the Fourier inversion formula

For the interpretation of this formula let us note that the integral regarding $$\vartheta$ only plays a role in the last transform. Until this point the variable $$\vartheta$ can be regarded as a simple parameter. It can be regarded as such until the point when the inverse Fourier transform delivers new spatial coordinates to bereplaced by the $$\boldsymbol{x\omega}$ expression. Now we write the state before the integration regarding $$\vartheta$ :

$$ \mathfrak{F}_{r}^{-1}\left [\mathfrak{F}_{t}\left [ \mathfrak{R}f\left ( t \right ) \right ]\left ( r )\left | r \right | \right ]$

This simplified expressions just the Fourier transfrom of the 1D function $$ \mathfrak{R}f\left ( t \right )$ , which is multiplied by the frequency variable |r| , and we do the inverse Fourier transform. This process describes a filtering in the Fourier-space, it is actually a high pass filtering.

The backprojection part of the Fourier inversion formula

After the filtering is done, we can start integrating according to $$ \vartheta$ . Let us use the notation g for the function obtained after the filtering: $$ (\boldsymbol{x}\boldsymbol{\omega},\boldsymbol{\omega})$ . Now the integral:
$$\int_{0}^{2\pi} g\left (\boldsymbol{\omega} \boldsymbol{x} \right ,\vartheta) d\vartheta =\int_{0}^{2\pi} g\left (x\cos\vartheta+y\sin\vartheta ,\vartheta) d\vartheta$

By definition $$(t=x\cos\vartheta+y\sin\vartheta,\vartheta)$ describes lines in the space of $$(t,\vartheta)$ , that cross point (x,y). The integral thus integrates the projections belonging to lines going through (x,y). This does not contradict our intuition since we collect all the information that belongs to the point (x,y). We can try how the backprojection of the Radon transform looks like without doing the filtering step. Let us start from distributions:

Image of a foal in (x,y) space	Radon transform of the picture of a foal

Let us backproject the Radon transform without the filtering step

Backprojection without filtering

You can see that the image is unacceptably soft, lacks details.
Let us use the following notation for the operator of the backprojection:
$$\mathfrak{R^{+}}$ .
Now we can write the steps of the filtered backprojection in operator form:
$$ {\mathfrak{R^{+}}\mathfrak{F}_{r}^{-1}\left [\mathfrak{F}_{t}\left [ \mathfrak{R}f \right ]\left ( r )\left | r \right | \right ]$

Example: steps of the filtered backprojection

Now let us go through the steps one-by-one. We have again the same photo for the 2D distribution we are trying to reconstruct from its sinogram:

Not its sinogram represents the acquired data:

With operators $$ \mathfrak{R}f$

Filtering must be done for each $$\vartheta$ . as an illustration we took projection number 200, that looks like:

With operator notation: $$ \mathfrak{R}f(t, 2\pi*\frac{200}{360})$

Now let us take the Fourier transform of this projection that yielded an amplitude spectrum as follows:

With operator notation: $$ \mathfrak{F}_{t}\left [\mathfrak{R}f(t, 2\pi*\frac{200}{360})\right ]\left ( r \right )$

Now multiplied by |r| we obtain the next amplitude spectrum:

With operator notation: $$ \mathfrak{F}_{t}\left [\mathfrak{R}f(t, 2\pi*\frac{200}{360})\right ]\left ( r \right )\left | r \right |$

Its inverse Fourier transform yields the filtered projection:

With operator notation: $$ \mathfrak{F}_{r}^{-1}\left [\mathfrak{F}_{t}\left [\mathfrak{R}f(t, 2\pi*\frac{200}{360})\right ]\left ( r \right )\left | r \right | \right ](t, 2\pi*\frac{200}{360})$
Now we can see that the DC (constant) component vanishes.

Now let us carry out the Fourier transform for every projection, then we get an angle-frequency sinogram that looks like:

With operator notation $$ \mathfrak{F}_{t}\left [\mathfrak{R}f(r,\vartheta )\right ]\left ( r \right ) ](t, \vartheta )$

Let us multiply every projection with |r|, then do inverse Fourier transform:

With operator notations: $$\mathfrak{F}_{r}^{-1}\left [\mathfrak{F}_{t}\left [\mathfrak{R}f(t,\vartheta )\right ]\left ( r \right ) ](r, \vartheta )\left | r \right | \right ](t,\vartheta )$
It is apparent that the image became sharper an more contrasty compared to the original image.

Taken the integral regarding $$\vartheta$ we obtain the distribution we sought for (the foal):

In operator notation: $$\mathfrak{R^{+}}\left [\mathfrak{F}_{r}^{-1}\left [\mathfrak{F}_{t}\left [\mathfrak{R}f(t,\vartheta )\right ]\left ( r \right ) ](r, \vartheta )\left | r \right | \right ](t,\vartheta ) \right ]\left ( x,y \right )$

Practical realization

In practice there are two main differences in the realization of the algorithm:

The filtering in Fourier domain is often replaced by a convolution
the filter |r| is not ideal, as it amplifies noises so different filter terms should be applied

These aspects are dealt with in the next section .

The original document is available at http://537999.nhjqzg.asia/tiki-index.php?page=The+Filtered+Backprojection