Fundaments of the 3-D Emission Imaging

In real clinical applications, beyond detecting the functional lesions, it may also be necessary to determine it’s size and location in the examined organ (e.g. brain, liver, heart, etc.). Therefore, the task is to image in 3-D the γ emitting observed organ. Technical realization of the proposed problem is based on the most novel imaging procedure of medical diagnostics: called Emission Computer Tomography. 2-D emission type of projection data sets (images) have to be acquired and registered from various angles around the 3D γ emitting object to be imaged - this is the investigated organ of the patient -. (see Figure 1.)

Figure 1.

The projection data set (in other words: raw data) is the fundament of the 3-D spatial imaging. If adequate number of projection data has been collected i.e. adequate fine angular sampling has been performed, then it is possible to produce in trans-axial plane so called cross-sectional images that are perpendicular to the rotation axis.
The operation is called image reconstruction procedure. The image reconstruction method can be considered as a mathematical problem and may be interpreted through purely mathematical models. Let us consider the model shown in Figure 2., where a “v” angle projection data, i.e. raw data of the spatial object described by bivariate function z=f(x,y) is presented in the (x,y,z) coordinate system (Cartesian coordinate system).

Figure 2.

The function g(s,v) is a “v” angle projection of f(x,y), to be derived as the line integral along “v” angle straight lines of f(x,y) function. [3], [4]

$$ g\left( s,v \right)=\int\limits_{-\infty }^{\infty }{f\left( x,y \right)du=}\int\limits_{-\infty }^{\infty }{F\left( s\cdot \cos v-u\cdot \sin v,s\cdot \sin v+u\cdot \cos v \right)du}$

where, due to the coordinate transformation:

$$ x=s\cdot \cos v-u\cdot \sin v\,\,\,\left( s=x\cdot \cos v+y\cdot \sin v \right) \\$

$$ y=s\cdot \sin v+u\cdot \cos v\,\,\,\left( s=y\cdot \cos v-x\cdot \sin v \right) \\$

$$ -\infty <s<\infty ;\,\,\,0\le v<\pi \\$

Consequently, the projection and back projection problem of a spatial object described by above mentioned bivariate function can be summarized on a clearly mathematical model as follow.
Let f(x,y) be a bivariate real function, where the Radon transformation $[\mathfrak{Rf}](s,v)$ exists and defined on the 0≤v<π interval:

$$ [\mathfrak{Rf}](s,v) = \mathfrak{R}\{f(x,y) \} = \int\limits^\infty_{-\infty } f(x,y) du \\$

, where s and u denote the variables in the v angle rotated coordinate system

From another approach, Radon transform assigns the projections of the function f(x,y) to the function. The various imaging systems - SPECT, PET, CT - perform transformations based on this model in various geometries.

Thus, the image reconstruction problem from a purely mathematical approach can be described as follow:

$[\mathfrak{Rf}](s,v)$ function is known and exists and let us determine -i.e. reconstruct- function z=f(x,y), where the $[\mathfrak{Rf}](s,v)$ function is the v angle projection in the 0≤v< $\pi$ interval.

Practically it is needed to find the inverse Radon transformation, i.e. which will transform the
$[\mathfrak{Rf}](s,v)$ function that is defined in a $\infty$ <s< $\infty$ ; 0≤v< $\pi$ set, into the set $\infty$ <x< $\infty$ ; $\infty$ <y< $\infty$ ; where the function is existing and defined:

$$ \mathfrak{B}\{ g(s,v) \} = \int\limits_0^\pi \int\limits_{-\infty}^{\infty} [\mathfrak{Rf}](s,v) g(s,v) ds dv = \int\limits_0^\pi \left( \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(s \ cos v-u \ sin v, s \ sin v + u \ cos v) du g(s,v)ds\right) dv$

Next change from the (s,u) coordinate to (x,y) coordinate and change the sequence of integrals as follow:

$$ \mathfrak{B}\{ g(s,v) \} = \int\limits_0^\pi \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x,y) g(x \ cos v+y \ sin v, v) dvdxdy = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x,y) \left( \int\limits_0^\pi g(x \ cos v+y \ sin v, v) dv\right) dxdy$

$$ \mathfrak{B}\{ g(s,v) \} = [\mathfrak{B}g](x,y) = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x,y) \left( \int\limits_0^\pi g(x \ cos v+y \sin v, v) dv\right) dxdy$

$$ \mathfrak{B}\{ g(x,y) \} = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x,y)dxdy \int\limits_{0}^{\pi} g(x \cdot cosv+y \cdot sinv,v)dv = BPI \int\limits_{0}^{\pi} g(x \cdot cosv+y \cdot sinv,v)dv$

The obtained formalism shows the inverse transform of $[\mathfrak{R}$ transforms the g(s,v) function defined on 0≤v<π set to the real value set defined ℬg(x,y) function:
$[\mathfrak{B}g](x,y)$ függvénybe viszi át

$$ [\mathfrak{B}g](x,y) = \mathfrak{B} \{ g(s,v) \}=BPI \int\limits_{0}^{\pi} g(x \cdot cosv+y \cdot sinv,v)dv$ (2)

, ahol (2) formulában $BPI = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x,y)dxdy$$ (3)
The meaning of inverse Radon transformation in formula (2) is the superimposition of the function g(s,v) that belongs to the individual v angles along the straight lines

$$ s=x\cos v+y\sin v$

ℬ transformation is called back projection transformation, i.e. LSBP (Linear Superposition of Back Projection).

Figure 3. presents the back projection of the function g(s,v) with respect to three angles. The “Figure” attempts to present how LSBP is working in case of such simple projection data. [3]

Figure 3.

Figure 4. illustrates a point source as a three-dimensional object, and shows the results of the back projection (LSBP) in case of three then more angular sampling.

Figure 4.

It is possible to read from the figure that with increasing the number of the views i.e. angular sampling the LSBP improves the tomography effect. However image blurring can’t be reduced below a certain value. The image will be blurred even in the case of infinite views. The phenomenon can be interpreted by system Point Spread Function (PSF) (see Chapter 3.3.1). The degree of blurring is proportional to 1/r, where r is the distance from the point source (Figure 5.). 2, 3

Figure 5.

Taking blurring into consideration, the following mathematical relation is true. [2], [3], [4]

$$ [\mathfrak{B}g](x,y)\{ I \} = RI *(1/r)$ (3)

where * denotes the operation of convolution. In formula (3) the {I}is the - 2D - raw data set to be back projected, while RI represents the real, blur less reconstructed object.

Henceforth, through formal mathematical operations it is attempted to show how to best approximate the RI blurless reconstructed image from the image ℬg(x,y){I} taking the system image blur factor into consideration.

Let us perform the Fourier transform of both sides of formula (3).
$$ \mathfrak{F} \{ [\mathfrak{B}g](x,y)\{ I \} \} = \mathfrak{F} \{ RI *(1/r) \}$
$$ \mathfrak{F} \{ [\mathfrak{B}g](x,y)\{ I \} \} = \mathfrak{F} \{ RI \} \mathfrak{F} \{ (1/r) \}$\\$ /due to the convolution operation rule/
$$ \mathfrak{F} \{ RI \} = \mathfrak{F} \{ [\mathfrak{B}g](x,y)\{ I \} \} \mathfrak{F} \{ (1/r) \}$

The Fourier transform of 1/r is tipical of the image blurring factor to be nothing else than the spatial frequency.
$$ \mathfrak{F} \{ (1/r) \} = 1/\Omega$
thus $$\mathfrak{F} \{ RI \} = \Omega \ \mathfrak{F} \{ [\mathfrak{B}g](x,y)\{ I \} \}$ (4)::
Performing inverse Fourier transformation on both side of equation (4), the blur less image is obtained in the real spatial domain as follow:
$$ RI = \mathfrak{F^{-1}}\{ \Omega \} * [\mathfrak{B}g](x,y)\{ I \} = h * [\mathfrak{B}g](x,y)\{ I \}$ (5)::
where $\mathfrak{F^{-1}}\{ \Omega \} = h$ is called the convolution filter factor.

Whether the formula (4) or formula (5) is used, the image blurring effects can be significantly reduced depending on the raw image data noise level and the degree of distortion. The image reconstruction method described by formulae (4) and (5) is called Linear Superposition of Filtered Back Projection (LSFBP).

Nowadays the LSFBP method is the most frequently used image reconstruction method in SPECT imaging technique and reference as well. By taking into consideration the noises and distortions originating from the whole system as well as from the image acquisition technique, image quality can be considerably improved by well-designed filter functions. Even today extra research and development work is necessary to find the optimal dedicated filter function for a particular SPECT system with the investigated organ.

Every final filter function has to contain the simplest filter function that belongs to the ideal case i.e. considers only the PSF blurring. This is called RAMP FILTER component. Figure 6. presents the frequency response. Figure 7. shows the three-dimensional spatial imaging of a point source activity, and its reconstruction by LSBP and LSFBP algorithms with the selection of an ideal – RAMP FILTER – filter function.

$$ FINAL\_ Filt[\Omega] = RAMP \_ Filt[\Omega] \{ USR\_ APOD \_ Filt[\Omega ] \}$ (6)

In case of point source object: Dirac Delta $\delta(\vec{r})$ – then appr. $$ USR\_ APOD \_ Filt[\Omega ] = 1$ (7)

Figure 6. (see more detailed in the chapter 2.6.6 ”Problems (8)” and 2.6.9 “Appendix” 2.6.9.2.8 S8 Solution of problems)

Figure 7.

The following animation (movie like presentations) figures (Figure 8., 9, 10.) show how LSFBP is working in case of a real 3D measured object (e.g. SPECT imaging). Both 2D raw data i.e. projection image set and the result of reconstruction as well as the way of reconstruction process are demonstrated.


Test iframe

Figure 8. 2D projection data of Tc99m-HMPAO labelled Brain SPECT imaging. SPECT image was acquired by parallel collimator with 6° angular sampling (i.e. 60pcs projections) and 64x64 matrix size. The “red marked” curve is the resultant filter curve in frequency domain for LSFBP (the convolution of the applied apodization and RAMP filters). The white curve is the same resultant filter curve, but presented in the real space domain.


Test iframe ::Figure 9. LSFBP in progress:	Figure 10. Result of the reconstruction displayed by optimized windowing level.

Finally, some typical evaluated and processed result of SPECT imaging are presented in the following figures (Figure 11., 12, 13.).

Figure 11.

Figure 12.

Figure 13.

The results of Brain SPECT stress/rest study is presented in Figure 11. in the three principal cross-sectional planes (transvere, frontal, sagittal). Figure 12. shows the summarized results of a myocardial perfusion (²⁰¹Tl) stress/rest study for the three main directions in relation to the heart (short axis, vertical long axis and horizontal long axis). Figure 13. displays the three-dimensional i.e. multi-dimensional functional parametric presentation (see detailed at 4.1, 4.3.5, 4.3.6, 4.4.6 chapters) of an ECG-gated “blood-pool” heart SPECT study. It is interesting to note that the result images practically contain five-dimensional information (multi-dimension): the three spatial coordinates, the phase and amplitude.

The original document is available at http://549552.cz968.group/tiki-index.php?page=Fundaments+of+the+3-D+Emission+Imaging