Filtered rear projection

Sectional image of a person, reconstructed with the filtered rear projection.

The filtered back projection (also FBP for filtered back projection ) is a method for image reconstruction based on the Radon transformation , which is primarily used in computer tomography . The FBP reconstructs an original 2D image from a set of one-dimensional projections in different directions. For this purpose, the projections are first filtered and then wiped across the screen in the respective direction ("back-projected"). The method has the great advantage that it is fast because it requires little computing power. However, errors in the imaging system cannot be corrected. In SPECT as well as in PET it has meanwhile been replaced by iterative reconstruction methods.

overview

Recording of projections by means of fluoroscopy with X-rays.

In numerous tomographic applications, the projection of an object in different directions can be measured, but not the interior of the object directly. Examples of this are the fluoroscopy of objects using X-rays , neutrons or ultrasound , or projective measurements in quantum mechanics . The set of projections in all directions contains the complete information about the interior of the object. However, it is difficult to interpret directly because each point consists of a superposition of properties that were added up during projection. The original image and the projections are linked by a one-to-one transformation that is mathematically described by the Radon transformation . The filtered back projection is an implementation of the transformation that allows the original image to be calculated from a set of projections.

algorithm

Schematic representation of the filtered back projection of a circular disk.

The filtered back projection is carried out on discrete data, i.e. pixel images of the projections for a finite number of evenly distributed (typically a few hundred) projection angles. The one-dimensional projections can also be combined to form a two-dimensional image, the so-called sinogram.

The filtered rear projection is based on the simple idea that the projections, similar to how they were projected out of the image, can also be projected back again. The naive implementation of the unfiltered rear projection does not work, however, because every point of the projection is smeared over the entire image area instead of just at its original 2D position. An image would be obtained in which each original point had a point spread function of the form 1 / | r | widened and thus corresponds to a very blurred version of the original image.

The functioning algorithm is obtained from the inverse Radon transformation . Be the original image and be the projection in the direction of the angle . Then ${\ displaystyle f (x, y)}$ ${\ displaystyle p _ {\ phi} (z)}$ ${\ displaystyle \ phi}$

{\ displaystyle f (x, y) = \ int _ {0} ^ {\ pi} \ left (\ int _ {- \ infty} ^ {\ infty} p _ {\ phi} (z ') \ cdot g ( z-z ') \, \ mathrm {d} z' \ right) \, \ mathrm {d} \ phi}

.

Here the coordinates and were expressed in polar coordinates . The integral over expresses that the contributions of all angles must be added up. The integral over is a convolution of the projection with a suitable high-pass filter , which gives the method its name. Indeed, inverse radon transformation is a poorly posed problem that calls for an irregular filter kernel . In order to obtain a discrete filter kernel, a slight soft focus must be accepted as regularization , which, depending on the parameter, leads to a different width . ${\ displaystyle x = z \ sin (\ phi)}$ ${\ displaystyle y = z \ cos (\ phi)}$ ${\ displaystyle \ phi}$ ${\ displaystyle z '}$ ${\ displaystyle p _ {\ phi} (z)}$ ${\ displaystyle g (z)}$ ${\ displaystyle g (z)}$ ${\ displaystyle g (z)}$

In practice, filtering is implemented as a fast convolution in Fourier space, where it takes on a particularly elegant form. There the (discrete) Fourier transform is multiplied by the filter in Fourier space , which has the particularly simple form there ${\ displaystyle {\ tilde {p}} _ {\ phi} (k) = {\ mathcal {F}} p _ {\ phi} (z)}$ ${\ displaystyle {\ tilde {g}} (k) = {\ mathcal {F}} g (z)}$

{\ displaystyle {\ tilde {g}} (k) = | k |}

accepts. All frequency components are weighted proportionally to the amount of their frequency. To avoid excessive noise in the result, an additional low-pass filter in the form of a window function is also used here , which suppresses the very high frequencies. The literature has produced dozens of different examples for the selection of suitable windows, for example Hamming, Hanning or Blackman. You also have to be careful with the DC component , which determines the mean value of . For discrete frequency bins the width must correspond to the bin mean value of , i.e. instead of 0. Using the Fourier transform and its inverse one obtains ${\ displaystyle w (k)}$ ${\ displaystyle k = 0}$ ${\ displaystyle g}$ ${\ displaystyle \ Delta k}$ ${\ displaystyle {\ tilde {g}} (0)}$ ${\ displaystyle | k |}$ ${\ displaystyle {\ tilde {g}} (0) = \ Delta k / 4}$ ${\ displaystyle {\ mathcal {F}}}$ ${\ displaystyle {\ mathcal {F}} ^ {- 1}}$

{\ displaystyle f (x, y) = \ int _ {0} ^ {\ pi} {\ mathcal {F}} ^ {- 1} {\ bigl (} {\ mathcal {F}} (p _ {\ phi }) (k) \ cdot | k | \ cdot w (| k |) {\ bigr)} \, \ mathrm {d} \ phi}

.

The implementation of the Fourier transform and its inverse is done efficiently as an FFT .

So the steps of filtered back projection are:

Transform each projection into Fourier space. ${\ displaystyle p _ {\ phi} (z)}$
Multiply the transform by and a window function . Set . ${\ displaystyle {\ tilde {p}} _ {\ phi} (k)}$ ${\ displaystyle {\ tilde {g}} (k) = | k |}$ ${\ displaystyle w (k)}$ ${\ displaystyle {\ tilde {g}} (0) = \ Delta k / 4}$
The filtered projection is transformed back into spatial space . ${\ displaystyle p '_ {\ phi} (z)}$
Back projection of into the image plane of , by stretching the 1D projection in a second dimension, then rotating it by the angle and finally adding it to the image. ${\ displaystyle p '_ {\ phi} (z)}$ ${\ displaystyle f (x, y)}$ ${\ displaystyle - \ phi}$

The sum of the back projections of all angles gives the picture . ${\ displaystyle f (x, y)}$

example

An example should clarify the algorithm :

Imagine an aquarium illuminated from the front, in which the contours of the fish are depicted on a canvas behind it. If you move the light source and the screen around the aquarium at the same angle, you get many projections of the fish on the screen. In order to determine the position of the fish in the aquarium afterwards from this data, one takes each of the recorded projections and projects them back onto the volume of the aquarium (hence the name of the procedure). It is clear that the depth information is not taken into account here; H. the finally deep image of the fish on the screen is smeared over the image in the projection direction. However, this error can be effectively suppressed by using a suitable image filter during reprojection.

The point spread function of the unfiltered back projection is , where the magnitude is in location space; This means that if the object to be imaged consists of only one point with the coordinates ( delta distribution ), the unfiltered back projection receives a signal at the location that is proportional to the image . The "filtering" corresponds mathematically to a convolution . With the help of the convolution theorem , unfolding can be carried out by transforming the back-projected image into Fourier space, multiplying it with the amount in Fourier space and then transforming it back into spatial space. According to the sampling theorem , the filter can be cut off at a certain spatial frequency. It should also be noted that data from a computer tomography are available discretely and not continuously, as is actually mathematically necessary. Therefore, the CT is "exactly" and there is not not the ideal filter, as for example, averaged between the discrete points (Shepp-Logan filter) or averaged (Ram-Lak filter) can be. Depending on which filter is used, the image is either richer in contrast, but noisier, or poorer in contrast, but reduced in noise. ${\ displaystyle 1 / r}$ ${\ displaystyle r}$ ${\ displaystyle \ left (x ', y', z '\ right)}$ ${\ displaystyle \ left (x, y, z \ right)}$ ${\ displaystyle {\ frac {1} {\ sqrt {(x-x ') ^ {2} + (y-y') ^ {2} + (z-z ') ^ {2}}}}}$ ${\ displaystyle | k | = {\ sqrt {k_ {x} ^ {2} + k_ {y} ^ {2} + k_ {z} ^ {2}}}}$