Loading...
 
PDF Print

The Filtered Backprojection

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
Image of a foal in (x,y) space
Image
Radon transform of the picture of a foal

 
Let us backproject the Radon transform without the filtering step

Image
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:
Image

Not its sinogram represents the acquired data:
Image
With operators $  \mathfrak{R}f

Filtering must be done for each $\vartheta. as an illustration we took projection number 200, that looks like:
Image
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:

Image
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:
Image
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:
Image
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:
Image
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:
Image
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):
Image
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:

  1. The filtering in Fourier domain is often replaced by a convolution
  2. 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 .


Site Language: English

Log in as…