Introduction to PET Physics



[Contents] [Section 1] [Section 2] [Section 3] [Section 4] [Section 5] [Section 6]


4. Image reconstruction

4.1 Introduction
4.2 Notation and mathematical theorems used
4.3 Analytic image formation in 2D PET
4.4 Filtered Back-Projection in 3D and 3D-RP

4.1 Introduction

After all corrections (e.g. for scatter, randoms and the effects of attenuation) have been applied to data acquired in a PET camera, the number of counts assigned to an LOR joining a pair of crystals is proportional to a line integral of the activity along that LOR. Parallel sets of such line integrals are known as projections. Reconstruction of images from projections is a problem to which much attention has been paid over the last 30 years, and many analytical and iterative reconstruction schema exist. For 2D reconstruction, the most commonly used algorithm is the analytical method called Filtered Back-projection (FBP). FBP is straight-forward to implement but does have the property of amplifying noise in the signal (see section 4.3). Recently, considerable interest has been shown in iterative reconstruction schema, such as the Ordered Subsets - Expectation Maximisation (OSEM) algorithm (Hudson and Larkin 1994), which possess different noise properties to FBP. For 3D reconstruction, the Reprojection and Filtered Back-projection (3D-RP) method of Kinahan and Rogers (1990) has been the most popular, in part because of the significant computational burden of newer 3D iterative reconstruction methods. 3D-RP itself is computationally expensive, and this has led to the development of approximate 3D reconstruction algorithms. Of these, Fourier Rebinning (Defrise et al 1997), which reduces the 3D problem to a series of 2D problems without significantly distorting the image and results in a significant reduction in the computational burden, is stimulating particular interest.

In this section FBP in 2 and 3 dimensions will be briefly summarised and the basic elements of 3D-RP will be described.
 

[top of page]



4.2 Notation and mathematical theorems used

In this work the following notation will be adopted:
 

 

 

Spatial variables

 

 

Fourier variables

 

 

vector
Cartesian components
polar components
vector
Cartesian components
polar components
1D
-
x
-
-
k
-
2D
x
x,y
r,q 
k
kx,ky
rk,qk
3D
x
x,y,z
r,q,j 
k
kx,ky,kz
rk,qk,jk
Table 3. Notation for spatial and Fourier quantities.
 
 

Fourier transform

inverse Fourier transform 

Back-projection 

Convolution operator

 

The spatial co-ordinates for a full-ring PET camera are shown in figure 11.

 

Figure 11. 3D co-ordinate system for a full-ring PET camera

Several mathematical theorems will be used without proof. These include the following:

Fourier's theorem:

                                             (5)

Fourier addition theorem:

                              (6)
 

 
 

Convolution theorem:

                                  (7a)

and

                                 (7b)

Shannon's sampling theorem:

If a function is sampled with a sampling distance d then it may be fully recovered from its samples (apart from harmonic components which are zero-valued at the sample points) if its Fourier transform has no non-zero components at frequencies beyond a particular value kc where kc = (1/2d).

The critical frequency kc = (1/2d) is known as the Nyquist frequency. A more detailed treatment of this material may be found in, for example,Bracewell (1986).
 

[top of page]



4.3 Analytic image formation in 2D PET

In 2D PET, data from LORs are arranged into sets of 1-dimensional parallel projections. For a point source distribution, this gives rise to a series of intensity profiles as shown in figure 12.
 
 

Figure 12. Projections generated from a single central point source (3 projections shown).

 

An estimate of the original source distribution may be obtained by a process known as back-projection. In this process, the magnitude of each value in a projection is added to every point in image space corresponding to the relevant line of integration in object space. Back-projections for a single point-source are shown in figure 13. When a small number of projections are used, the resultant image contains star-shaped artefacts. In the limit of an infinite number of projection angles, this process is the equivalent of convolving the original source distribution with the function f(r,q ) = |r|-1.

An extended source distribution can be considered as being made up of a series of points of varying intensity, each of which, after back-projection, is convolved with |r|-1. Since convolution is distributive over addition, the relationship between the density function and the back projection b(r,q ) may be written:

                                                       (8)

where r (r,q ) is the extended source distribution. Application of the convolution theorem gives:

                                       (9)

Now the Fourier transform of |r|-1 is just |rk|-1, so that

                                          (1.10)

and we have recovered the extended source distribution function. In practice the source distribution is sampled with a finite sample-width, so to enable application of the Sampling theorem it is necessary cut off the Fourier transform of the back-projection at the Nyquist frequency. This is achieved by multiplying it with a gate function of the appropriate width. In real space this is equivalent to convolution with the sinc function. This can result in "undershoot" artefacts in the image (where image values are underestimated or even negative), particularly in regions close to sharp edges in the source distribution.

 

Figure 13. Back-projections of a point source. With finite numbers of

back-projection angles, "star" artefacts are seen.

 

Filtered Back-projection is a mathematically equivalent process to that described above. Fourier transforms of the projections are first calculated and multiplied by the cut-off 1D version of |rk| (the ramp filter). The inverse Fourier transform of the result is then calculated and back-projected to create the image.

Use of the ramp filter amplifies high-frequency components in the back-projection. Unfortunately statistical noise in the data manifests in Fourier space as high-frequency components. So the process of FBP amplifies noise in the image. In order to reduce this effect, a range of modifications to the ramp filter can be used. Perhaps the most commonly used of these is the Hanning filter (figure 14). Such filters are equivalent to some form of smoothing in image space.

Figure 14. The Ramp and Hanning filters

 
A more detailed mathematical treatment of analytic 2D tomographic reconstruction is given by Brooks and Di Chiro (1976).
 
 

[top of page]



4.4 Filtered Back-Projection in 3D and 3D-RP
 

In three dimensions, the data from the LORs may be arranged into 2D sets of parallel projections (figure 15(b)). FBP generalises to 3D directly if the projections can be obtained over all j as well as q . Unfortunately in real cameras projections cannot easily be obtained over the full range of j , and a different filter must be found for the de-convolution step (Colsher,1980).
 

 

 

Figure 15(a). Parallel projections in 2D. Note that the LORs become closer together towards the edge of the FOV. To correct for this, the data must be re-sampled (arc corrected) prior to reconstruction.
 
 

 

Figure 15(b). Parallel projections in 3D

 

Another problem is the fact that the limited axial extent of the camera causes some of the 2D projections to be truncated, and the degree of truncation depends on position (figure 16). As a result there is no single filter which is appropriate for every projection set. The 3D-RP algorithm circumvents this problem by performing a 2D reconstruction on a subset of the 3D projections first. The resulting image volume is then reprojected in order to obtain estimates of the truncated projections and a 3D FBP is then performed on the combined real and synthesised data.

 

Figure 16. Axial cut-away diagram of a PET camera operating in 3D mode, showing the extent of the projection sets as a function of angle j . As j increases, the measurable extent of the projection set decreases, requiring the reconstruction filter to change with position. To avoid this, an initial 2D reconstruction is performed on the j = 0 projection set, and estimates of the missing parts of the truncated projections obtained by reprojecting through the image volume.
 
 
 

 
 
 
 

[top of page]

Last revised by:

Ramsey Badawi

Revision date:

12 Jan 1999