*1.1.1. Projection and scanning of the object*

The methodological basis used for describing 2D and 3D image reconstruction was presented in detail in the work of Kak and Slantay [1] and Rosenfeld and Kak [2]. Other important work such as that of Herman [3], describes the reconstruction methods such as algebraic reconstruc‐ tion technique (ART). In this work, we focus on the work and analytical methods of Kak and Stanlay to explain the different steps of tomography from the measurement to the reconstruc‐ tion of a single layer. The superposition of such layers will constitute the 3D image or volume representing the studied object. By means of image processing tools, parts of the reconstructed 3D volume can be extracted and analyzed separately from the rest of the data set.

mation storage associated. It becomes thus possible by appropriate processing to detect the presence of defects, to identify the internal structures and to study their form and their posi‐ tion, to quantify the variations of density, to model the internal components and to guide the instruments of intervention in medicine. Finally, the user will be also able to take benefit from a large variety of existent software and algorithms for the tomography digital images

82 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

A CT system gathers several technological components. Its development requires the partic‐ ipation of the end users - as the doctors, the physicists or the biologists - to specify the needs, of the engineers and researchers to develop the novel methods and, finally, of the in‐

In this work, we are interested in the physics and mathematics related to the main phases of computed tomography, namely: the scanning or projection phase and the phase of 2D or 3D image reconstruction by various analytical methods such filtered back-projection method (FBP). In this context, all the mathematical equations and relations which seemed ambigu‐ ous or not clear are explained and mathematically demonstrated with some illustrative ex‐ amples. This chapter is organized in the form of a main body text with sub-sections presenting a brief explanation of most important concept or the detailed demonstration of any equation or expression, and this, each time that seemed necessary. The algebraic meth‐

This chapter is intended as an introduction to Computed Tomography. It was written not only for those persons who have some familiarity with other imaging techniques such radia‐ tion transmission radiography but also for novices in the field of digital imaging. The chap‐ ter begins with some simple, yet fundamental, concepts regarding computed tomography and the physics and mathematics at the origin of CT. As one progresses through the chapter, more detail regarding the CT technique and methodology is meet. The reader should not be alarmed if his or her particular problem or preoccupation is not mentioned here. So many different CT developments have been achieved in the last thirty-three years that it would be difficult to describe them all in a single chapter. Some practical information are also present‐ ed that can be vital to obtaining good analytical results; it is sometimes difficult to find.

I hope that this introduction to the CT technique will provide useful information to those persons who are about to get involved with CT as well as present TC users and those with

The methodological basis used for describing 2D and 3D image reconstruction was presented in detail in the work of Kak and Slantay [1] and Rosenfeld and Kak [2]. Other important work such as that of Herman [3], describes the reconstruction methods such as algebraic reconstruc‐ tion technique (ART). In this work, we focus on the work and analytical methods of Kak and Stanlay to explain the different steps of tomography from the measurement to the reconstruc‐

processing, analysis and visualization.

simply a curiosity about the technique.

*1.1.1. Projection and scanning of the object*

dustrial teams to develop, produce and market these systems.

ods for image reconstruction in tomography will be also outlined.

**1.1. Analytical methods for image reconstruction in tomography**

An object O (x, y, z) is considered as a superposition of n layers of the same thickness along z axis, all located in planes parallel to the plane (x, y) and perpendicular to z (Fig. 1). Each layer represents a section in the object to be reconstructed. It is considered as a 2D function fn(x, y) that describes, for example, the distribution of linear attenuation coefficients as a function of the position or any other 2D function that can be measured and whose measurement signal is described by a full line. Any function maybe considered in tomography in condition to be limited and finite in a given region and equal to zero outside this region. The condition of a finite size is easily achievable for solid samples, liquid or gas contained in a box. The estab‐ lishment of condition is very difficult when it will be used for the tomography of electric or magnetic fields. The generalization of such function for the tomographic measurements is closely related to the determination of the nature of the interaction of the scanning beam (combshaped) with the object under examination. The purpose of tomography is the reconstruction of this 2D function, representing a layer or slice of the object, from the measured projections in a unique way.

**Figure 1.** The geometry of a studied object scanning in the {x,y,z} coordinates system. A layer in the plane (x, y) is scan‐ ned along the angle θ and the transmitted intensity is stored in a system (t, s) of rotational coordinates.

The layer is scanned (scanned) in an angle θ varying from 0 ° to 180 ° for the transmission tomography. The intensity of the transmitted beam is recorded like a translation function of of the position parameter t (Fig. 2). The transmitted intensity is given by Lambert's law given by the following expression:

$$I(\mathbf{x}, y) = I\_0 \exp(-\int\_{path} \mu(\mathbf{x}, y) ds) \tag{1}$$

**Figure 2.** Scanning of a single layer in the plane (x, y). Note that z is the axis of rotation which coincides with the z-axis of the rotational coordinate's system {s, t, z}.

Where I0 is the incident beam intensity and μ(x, y) = f(x,y) is the 2D function to rebuild. A new square and rotational coordinates system (t, s) is defined to express the detection system rotatable in comparison to the fixed object coordinates system (or vice versa, if the object is rotating and the detection system is fixed). In the transformation from the system (x, y) to the system (t, s), t is given by (see demonstration 2):

$$t = \text{x.} \cos(\theta) + \text{y.} \sin(\theta) \tag{2}$$

d

extension to the case of diverging beam (fun beam) is quiet easy.

*P*θ(*t*)= ∫ −∞

*JF is the Jacobian matrix determinant which is given by:*

object function μ (x, y) ((μ (x, y) = μ (s, t)).

*JF* =

*So, in our case dsdt=dxdy, and we can write the following equation:*

*t s*

*P*θ(*t*)= ∫ *ray*(θ,*t*)

+∞

= ∫ −∞

+∞ ∫ −∞

∂ ∂ *x* ∂ ∂ *y*

cosθ −sinθ

+∞ ∫ −∞

+∞

*Eq.D.1:*

*With:*

*expression.*

+¥ +¥


 q

 qm( cos sin ) ( , ) *x y t x y dxdy*

This last expression is just one of the different forms of the Radon transform basically used for determining a function from its integral according to certain directions (see math reminder 1). The two-dimensional Radon transform projects an object f (x, y) to get its projections P<sup>θ</sup> (t). The projections values depend on of the integral of the object values along the line of integral according to a direction θ. In this work, we are interested in the case of a parallel beam; the

*Demonstration 1 Demonstrating that the projection P*θ *(t) as defined by Eq.1 can be written in one form of the Radon transform given by*

*For a convenience and simplicity reasons and for the fact that tomography is a process based on a rotational scanning, a new square and rotational coordinates system (t, s) is defined as presented above to take into consideration the rotating aspect of the source and the detection system around the fixed object. The coordinates of this system are given by:*

{ *<sup>t</sup>* <sup>=</sup> *<sup>x</sup>*.cos<sup>θ</sup> <sup>+</sup> *<sup>y</sup>*.sin<sup>θ</sup>

sinθ cosθ

μ(*x*, *y*)*ds* = ∫

=cos<sup>2</sup>

*ray*(θ,*t*)

*Note that the coordinates system transformation from (x, y) system to (t, s) system has no influence on the values of the*

*End of Demonstration 1.*

δ(*x*cosθ + *y*sinθ −*t*)μ(*s*, *t*)*dtds*

*Thus, it was demonstrated that the tomography projection can be expressed by one of the Radon's transform*

θ + sin<sup>2</sup>

μ(*s*, *t*)*ds*

= +- ò ò (4)

http://dx.doi.org/10.5772/52351

85

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

δ(*x*cosθ + *y*sinθ −*t*)μ(*x*, *y*)*dxdy* (D.1)

*<sup>s</sup>* <sup>=</sup> <sup>−</sup> *<sup>x</sup>*.sin<sup>θ</sup> <sup>+</sup> *<sup>y</sup>*cos<sup>θ</sup> (D.2)

*dsdt* = *JFdxdy* (D.3)

θ =1 (D.4)

(D.5)

The expression of the ray (path) through the sample expressed in terms of t and θ and by the substitution of s is: δ(t-x.cos (θ) + x.sin (θ)). The Dirac function δ ensures that only the points in obedience to the equation (2) that are related to the beam (comb-shaped) contribute to the projection Pθ (t). Such a projection can be defined by:

$$\mathbf{P}\_{\partial}(\mathbf{t}) = \ln(\underset{\mathbf{I}\_{0}}{\overset{\text{Id}}{\prod}}\_{\text{path}}) = \int \boldsymbol{\mu}(\mathbf{x}, \mathbf{y}) \, \mathrm{ds} \tag{3}$$

Note here that the expression of Pθ (t) can be given by (see demonstration 1):

$$\mathcal{P}\_{\Theta}(\mathbf{t}) = \int\_{\text{path}} \mu(\mathbf{x}, \mathbf{y}) \, \text{ds}.$$

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples http://dx.doi.org/10.5772/52351 85

$$=\int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} \delta(x\cos\theta + y\sin\theta - t)\mu(x,y)dxdy\tag{4}$$

This last expression is just one of the different forms of the Radon transform basically used for determining a function from its integral according to certain directions (see math reminder 1). The two-dimensional Radon transform projects an object f (x, y) to get its projections P<sup>θ</sup> (t). The projections values depend on of the integral of the object values along the line of integral according to a direction θ. In this work, we are interested in the case of a parallel beam; the extension to the case of diverging beam (fun beam) is quiet easy.

#### *Demonstration 1*

*Demonstrating that the projection P*θ *(t) as defined by Eq.1 can be written in one form of the Radon transform given by Eq.D.1:*

$$\mathcal{P}\_{\theta}(t) = \int\_{-\infty}^{+\infty} \int \delta(\chi \cos \theta + \chi \sin \theta - t) \mu(\chi, \chi) d\chi d\chi \tag{D.1}$$

*For a convenience and simplicity reasons and for the fact that tomography is a process based on a rotational scanning, a new square and rotational coordinates system (t, s) is defined as presented above to take into consideration the rotating aspect of the source and the detection system around the fixed object. The coordinates of this system are given by:*

$$\begin{cases} \mathbf{t} = \times \cos \theta + \chi \sin \theta\\ \mathbf{t} = -\times \cos \theta + \chi \cos \theta \end{cases} \tag{\text{D.2}}$$

*With:*

Where I0 is the incident beam intensity and μ(x, y) = f(x,y) is the 2D function to rebuild. A new square and rotational coordinates system (t, s) is defined to express the detection system rotatable in comparison to the fixed object coordinates system (or vice versa, if the object is rotating and the detection system is fixed). In the transformation from the system (x, y) to the

**Figure 2.** Scanning of a single layer in the plane (x, y). Note that z is the axis of rotation which coincides with the z-axis

y

84 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

t

θ

*tx y* = + .cos( ) .sin( ) q

The expression of the ray (path) through the sample expressed in terms of t and θ and by the substitution of s is: δ(t-x.cos (θ) + x.sin (θ)). The Dirac function δ ensures that only the points in obedience to the equation (2) that are related to the beam (comb-shaped) contribute to the

 q

(2)

x

μ(x,y)ds (3)

system (t, s), t is given by (see demonstration 2):

of the rotational coordinate's system {s, t, z}.

Pθ(t)= ln (I/I0)

projection Pθ (t). Such a projection can be defined by:

Pθ(t)= *∫* path

μ(x,y)ds

<sup>P</sup>θ(t)=ln(<sup>I</sup>

I0 )= *∫* path

Note here that the expression of Pθ (t) can be given by (see demonstration 1):

$$\mathbf{r} \times \mathbf{c} \mathbf{d} \times \mathbf{t} = \mathbb{J}\_F \mathbf{d} \times \mathbf{c} \mathbf{t} \tag{\text{D.3}}$$

*JF is the Jacobian matrix determinant which is given by:*

$$\mathbf{J}\_{\mathcal{A}} = \begin{bmatrix} \cos \Theta & \cos \Theta \end{bmatrix} \tag{2.4}$$
 
$$\mathbf{J}\_{\mathcal{A}} = \begin{bmatrix} \cos \Theta & \cos \Theta \end{bmatrix} \tag{2.5}$$

*So, in our case dsdt=dxdy, and we can write the following equation:*

$$\begin{aligned} \mathcal{P}\_{\theta}(t) &= \int\_{s\rho(\theta,t)} \mu(\mathbf{x}, \mathbf{y}) d\mathbf{x} = \int\_{s\rho(\theta,t)} \mu(\mathbf{x}, \mathbf{t}) d\mathbf{x} \\ &\quad \mathop{\rm str\underset{\ast}{\to}}\limits\_{s\to\infty} \\ &= \int\_{-\infty} \int\_{-\infty} \mathcal{S}(\mathbf{x}\cos\Theta + \jmath s\sin\Theta - t)\mu(\mathbf{x}, \mathbf{t}) d\mathbf{t} d\mathbf{x} \end{aligned} \tag{D.5}$$

*Note that the coordinates system transformation from (x, y) system to (t, s) system has no influence on the values of the* object function μ (x, y) ((μ (x, y) = μ (s, t)).

*Thus, it was demonstrated that the tomography projection can be expressed by one of the Radon's transform expression.*

#### *End of Demonstration 1.*

Before continuing our mathematical development and demonstration specific to the object scanning in transmission tomography, it is useful to give some math reminder on Radon transform and its main properties.

#### *Math Reminder 1.*

*The Radon transform is defined by:*

$$\begin{aligned} \mathcal{R}(\boldsymbol{\mu}, \mathbf{r}) \mathbf{I} \{ f \left( \boldsymbol{\chi}, \boldsymbol{\chi} \right) \} &= \Big|\_{ - \infty }^{ + \infty } f \left( \boldsymbol{\chi} + \mathbf{r} + \boldsymbol{\mu} \boldsymbol{\chi} \right) d \boldsymbol{\chi} \\\\ \mathbf{I} = \Big|\_{ - \infty }^{ + \infty } f \left( \boldsymbol{\chi}, \boldsymbol{\chi} \right) \mathsf{d} \{ \boldsymbol{\chi} - \left( \mathbf{r} + \boldsymbol{\mu} \boldsymbol{\chi} \right) \} \mathsf{d} \boldsymbol{\chi} d \boldsymbol{\chi} \succeq \mathsf{U} \left( \boldsymbol{\mu}, \ \mathsf{r} \right), \end{aligned} \tag{8.1}$$

*p* =cotα, τ =cscα (R.10)

<sup>1</sup> <sup>+</sup> *<sup>p</sup>* <sup>2</sup> , <sup>α</sup> =cot−1*<sup>p</sup>* (R.11)

*R*(*p*, τ) *f* 1(*x*, *y*) + *f* <sup>2</sup>(*x*, *y*) =*U*1(*p*, τ) + *U*2(*p*, τ); (R.12)

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

*R*(*p*, τ) *af* (*x*, *y*) =*aU*(*p*, τ); (R.13)

*d* −*b*(*c* + *bd*)

*R*(*p*, τ) *f* (*x*, *y*)∗*g*(*y*) =*U*(*p*, τ)∗*g*(τ); (R.18)

*U*(*p*, τ); (R.17)

*f* (*x*, *y*)*dxdy*; (R.19)

(*x*, *y*)*dxdy*; (R.20)

*<sup>b</sup>* ); (R.14)

http://dx.doi.org/10.5772/52351

87

cos<sup>Φ</sup> <sup>+</sup> *<sup>p</sup>*sin<sup>Φ</sup> ), (R.15)

*<sup>a</sup>* <sup>+</sup> *bp* ); (R.16)

*<sup>r</sup>* <sup>=</sup> <sup>τ</sup>

*R*(*p*, τ) *f* 1(

*<sup>R</sup>*(*p*, <sup>τ</sup>) *<sup>R</sup>*<sup>Φ</sup> *<sup>f</sup>* (*x*, *<sup>y</sup>*) <sup>=</sup> <sup>1</sup>

*<sup>R</sup>*(*p*, <sup>τ</sup>) *<sup>f</sup>* (*ax* <sup>+</sup> *by*, *cx* <sup>+</sup> *dy*) <sup>=</sup> <sup>1</sup>

∫−∞ +∞

∫−∞ +∞

*x a* , *y*

*<sup>b</sup>* ) <sup>=</sup> <sup>|</sup>*<sup>a</sup>* <sup>|</sup>*U*(*<sup>p</sup> <sup>a</sup>*

<sup>|</sup> cos<sup>Φ</sup> <sup>+</sup> *<sup>p</sup>*sinΦ<sup>|</sup> *<sup>U</sup>*( *<sup>p</sup>* <sup>−</sup>tan<sup>Φ</sup>

<sup>|</sup>*<sup>a</sup>* <sup>+</sup> *bp* <sup>|</sup> *<sup>U</sup>*(

*I* = 1 + *p* <sup>2</sup>

*<sup>U</sup>*(*p*, <sup>τ</sup>)*d*<sup>τ</sup> <sup>=</sup>∫−<sup>∞</sup>

*R*(*p*, τ) *f* (*x*, *y*) <sup>2</sup>

+∞ ∫−∞ +∞

*<sup>d</sup>*<sup>τ</sup> <sup>=</sup>∫−<sup>∞</sup> +∞ ∫−∞ +∞ *f* 2

*End of Math Reminder 1.*

At this level of mathematical development a very important question must be asked: Why the switching between different coordinates systems, presented above (Fig.1), is so important to establish the mathematical basis of image projection and reconstruction in CT in the case of analytical method such as FBP? The response to this question is satisfactory explained in the

*b* , τ

<sup>1</sup> <sup>+</sup> *<sup>p</sup>*tan<sup>Φ</sup> , <sup>τ</sup>

*c* + *dp <sup>a</sup>* <sup>+</sup> *bp* , <sup>τ</sup>

*The main properties of Radon transform are the followings:*

*1. Superposition :*

*2. Linearity:*

*3. Scaling:*

*4. Rotation:*

*5. Skewing:*

*R*Φ *is a rotation operator;*

*6. Integral along p:*

*7. 1D convolution expression:*

*8. Equivalent equation of Plancherel theorem:*

*9. Equivalent equation of Parseval theorem:*

2nd demonstration (see Demonstration 2).

where p is the slope of the projection line of and τ is its intersection with the y axis. The inverse Radon transform is given *by:*

$$\text{If } f(\chi, \chi) = \frac{1}{2\pi} \left[ \int\_{-\infty}^{+\infty} \frac{\mathcal{O}}{\mathfrak{d}\mathcal{Y}} \# \coprod \mathsf{U} \, \mathsf{U} \, (\mathsf{p}, \, \mathsf{y} - \mathsf{p}\mathsf{x}) \, \mathsf{Id}\mathfrak{p} \,. \tag{\mathsf{R.2}}$$

where H is the Hilbert transform. The Radon transform can also be defined by the following expression:

$$\left[\mathcal{R}\left(r,\mathbf{a}\right)\right]f\left(\mathbf{x},\,\,\mathbf{y}\right)\mathbf{j}=\Big|\_{r=0}^{+\infty}\Big|\_{r=0}^{+\infty}f\left(\mathbf{x},\,\,\mathbf{y}\right)\mathbf{\hat{S}}\left(r-\mathbf{x}\cos\mathbf{a}-\mathbf{y}\sin\mathbf{a}\right)\mathbf{\hat{k}}\,\mathbf{c}\,\mathbf{\hat{y}}\,\mathbf{j}\,\,\tag{\text{R.3}}$$

Where r is the perpendicular distance of the integral line with respect to the origin and α is the angle formed by this line and the x axis.

Using the following identification:

$$\mathbb{P}\_{\omega,\mathsf{a}}\mathsf{f}\mathsf{R}\mathsf{f}\mathsf{f}\mathsf{f}\left(\omega,\mathsf{a}\right)\mathsf{T}\mathsf{f}\left(\mathsf{x},\mathsf{y}\right) = \mathsf{F}\_{\omega,\mathsf{b}}^{2}\mathsf{f}\mathsf{f}\left(\omega,\mathsf{v}\right)\mathsf{f}\left(\mathsf{x},\mathsf{y}\right),\tag{\mathsf{R}.4}$$

where F is the Fourier transform, the Radon inversion formula can be expressed by:

*<sup>f</sup>* (*x*, *<sup>y</sup>*)=*c*∫<sup>0</sup> π ∫−∞ +∞ *<sup>F</sup>*ω,<sup>α</sup> *<sup>R</sup> <sup>f</sup>* (ω, <sup>α</sup> <sup>|</sup><sup>ω</sup> <sup>|</sup>*<sup>e</sup> <sup>i</sup>*ω(*x*cosα+*y*sinα) *d*ω*d*α (R.5)

*This last expression can be simplified as follows:*

$$f \begin{pmatrix} \chi & \chi \end{pmatrix} = \left[ \int\_0^\pi \begin{bmatrix} \* & \mathbf{a} \\ - \infty & \mathbf{I} \end{bmatrix} f \begin{pmatrix} r & \mathbf{a} \end{pmatrix} \mathbf{M} \begin{pmatrix} r & \mathbf{a} \end{pmatrix} \chi \right] d\mathbf{r} d\mathbf{a},\tag{8.6}$$

where W is a weight function given by:

$$\mathcal{W}\{r,\mathfrak{a},\times,\chi\} = \hbar \left(\times \cos \mathfrak{a} + \mathcal{y} \sin \mathfrak{a} - r\right) = \mathbb{F}^{-1} \mathfrak{f} \left| \mathfrak{w} \right| \mathfrak{f} \tag{\mathbb{R} \ \mathcal{I}}$$

Nievergelt (1986) determined the inversion formula as follows:

$$f \text{ ( $\chi$ ,  $\chi$ )}=\lim\_{\Pi} \lim\_{\varepsilon \to 0} \Big|\_{\mathfrak{U}} \Big|\_{\mathfrak{U}} \Big|\_{\mathfrak{U}} \stackrel{\star \infty}{\mathfrak{R}} f \, (\mathfrak{r} + \chi \text{cos}\mathfrak{a} + \chi \text{sin}\mathfrak{a}, \mathfrak{a}) \Big| \! \mathfrak{G}\_{\varepsilon} \langle r \rangle \text{d}r \text{d}a,\tag{8.8}$$

*with:*

$$\mathbf{G}\_c(r) = \begin{cases} \frac{1}{\pi \mathbf{c}^2} & \text{for } |r| \le c \\\\ \frac{1}{\pi \mathbf{c}^2} \{1 - \frac{1}{\sqrt{1 - \mathbf{c}^2/r^2}}\} & \text{for } |r| > c \end{cases} \tag{8.9}$$

The Ludwig's inversion formula presents the relations between the two forms of the Radon transform of a given function: R (r, α) and R*'* (r, α), which are given by:

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples http://dx.doi.org/10.5772/52351 87

$$
\mathfrak{p} = \mathfrak{c} \circledast \mathfrak{a}, \quad \mathfrak{r} = \mathfrak{c} \mathsf{c} \mathsf{c} \mathsf{a} \tag{\mathfrak{R}.10}
$$

$$\sigma = \frac{\pi}{1 + \rho^2}, \quad \sigma = \cot^{-1} \rho \tag{8.11}$$

*The main properties of Radon transform are the followings:*

#### *1. Superposition :*

$$\mathcal{R}(\mathfrak{p}, \mathfrak{r}) \Big[ f\_1(\mathbb{X}, \mathcal{Y}) + f\_2(\mathbb{X}, \mathcal{Y}) \Big] = \cup\_1 \langle \mathfrak{p}, \mathfrak{r} \rangle + \cup\_2 \langle \mathfrak{p}, \mathfrak{r} \rangle;\tag{8.12}$$

*2. Linearity:*

Before continuing our mathematical development and demonstration specific to the object scanning in transmission tomography, it is useful to give some math reminder on Radon

86 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

*Math Reminder 1.*

*f* (*x* + τ + *px*)*dx*

where p is the slope of the projection line of and τ is its intersection with the y axis. The inverse Radon transform is given

Where r is the perpendicular distance of the integral line with respect to the origin and α is the angle formed by this

*<sup>F</sup>*ω,<sup>α</sup> *<sup>R</sup> <sup>f</sup>* (ω, <sup>α</sup> <sup>|</sup><sup>ω</sup> <sup>|</sup>*<sup>e</sup> <sup>i</sup>*ω(*x*cosα+*y*sinα)

<sup>π</sup>*<sup>c</sup>* <sup>2</sup> *for* <sup>|</sup>*<sup>r</sup>* <sup>|</sup> <sup>≤</sup>*<sup>c</sup>*

<sup>1</sup>−*<sup>c</sup>* <sup>2</sup> / *<sup>r</sup>* <sup>2</sup> *for* <sup>|</sup>*<sup>r</sup>* <sup>|</sup> <sup>≻</sup>*<sup>c</sup>*

The Ludwig's inversion formula presents the relations between the two forms of the Radon transform of a given

*H U*(*p*, *y* − *px*) *dp*, (R.2)

*<sup>f</sup>* (*x*, *<sup>y</sup>*)δ(*<sup>r</sup>* <sup>−</sup> *<sup>x</sup>*cos<sup>α</sup> <sup>−</sup> *<sup>y</sup>*sinα)*dxdy*, (R.3)

<sup>2</sup> *f* (*u*, v) (*x*, *y*), (R.4)

*R f* (*r*, α) *W* (*r*, α, *x*, *y*)*drd*α, (R.6)

*R f* (*r* + *x*cosα + *y*sinα, α) *Gc*(*r*)*drd*α, (R.8)

*W* (*r*, α, *x*, *y*)=*h* (*x*cosα + *y*sinα −*r*)= *F* <sup>−</sup><sup>1</sup> |ω | (R.7)

*d*ω*d*α (R.5)

*f* (*x*, *y*)δ *y* −(τ + *px*) *dxdy* ≡*U*(*p*, τ),

*<sup>R</sup>*(*p*, <sup>τ</sup>) *<sup>f</sup>* (*x*, *<sup>y</sup>*) <sup>=</sup>∫−<sup>∞</sup>

*<sup>f</sup>* (*x*, *<sup>y</sup>*)= <sup>1</sup>

(*r*, <sup>α</sup>) *<sup>f</sup>* (*x*, *<sup>y</sup>*) <sup>=</sup>∫−<sup>∞</sup>

*<sup>f</sup>* (*x*, *<sup>y</sup>*)=*c*∫<sup>0</sup>

Nievergelt (1986) determined the inversion formula as follows:

*<sup>f</sup>* (*x*, *<sup>y</sup>*)= <sup>1</sup>

*This last expression can be simplified as follows:*

where W is a weight function given by:

2π ∫−<sup>∞</sup> <sup>+</sup><sup>∞</sup> *d dy*

+∞ ∫−∞ +∞

*F*ω,<sup>α</sup> *R f* (ω, α) (*x*, *y*)= *Fu*,<sup>v</sup>

where F is the Fourier transform, the Radon inversion formula can be expressed by:

π ∫−∞ +∞

*<sup>f</sup>* (*x*, *<sup>y</sup>*)=∫<sup>0</sup> π ∫−∞ +∞

> <sup>π</sup> lim *<sup>c</sup>*→<sup>0</sup> ∫<sup>0</sup> π ∫−∞ +∞

(r, α), which are given by:

*Gc*(*r*)={ <sup>1</sup>

1 <sup>π</sup>*<sup>c</sup>* <sup>2</sup> (1<sup>−</sup> <sup>1</sup>

where H is the Hilbert transform. The Radon transform can also be defined by the following expression:

<sup>=</sup>∫−<sup>∞</sup> +∞ ∫−∞ +∞

*R* ′

+∞

transform and its main properties.

*The Radon transform is defined by:*

*by:*

line and the x axis.

*with:*

function: R (r, α) and R*'*

Using the following identification:

*R*(*p*, τ) *af* (*x*, *y*) =*aU*(*p*, τ); (R.13)

*3. Scaling:*

(R.1)

(R.9)

$$\mathcal{R}(\boldsymbol{\rho}, \,\mathrm{tr}\Big[\boldsymbol{f}, \,\mathrm{i}\Big(\frac{\boldsymbol{\chi}}{\mathrm{a}}, \,\frac{\boldsymbol{\chi}}{\mathrm{b}}\Big)\Big] = \boldsymbol{\parallel} \,\mathrm{a} \,\big|\,\boldsymbol{\cup} \,\mathrm{(\boldsymbol{\rho}} \,\frac{\mathrm{a}}{\mathrm{b}}, \,\frac{\mathrm{r}}{\mathrm{b}}\Big);\tag{8.14}$$

*4. Rotation:*

$$\mathcal{R}(\boldsymbol{\rho},\,\boldsymbol{\tau}) \Big\| \mathcal{R}\_{\Phi} f(\boldsymbol{\nu},\,\boldsymbol{\nu}) \Big\| = \frac{1}{|\cos\mathfrak{\boldsymbol{\sigma}} + \rho\sin\mathfrak{\boldsymbol{\sigma}}|} \cup (\frac{\mathfrak{\boldsymbol{\sigma}} - \tan\mathfrak{\boldsymbol{\sigma}}}{1 + \rho\tan\mathfrak{\boldsymbol{\sigma}}}, \frac{\mathfrak{\boldsymbol{\tau}}}{\cos\mathfrak{\boldsymbol{\sigma}} + \rho\sin\mathfrak{\boldsymbol{\sigma}}}),\tag{8.15}$$

*R*Φ *is a rotation operator;*

#### *5. Skewing:*

$$\mathcal{R}\{\mathfrak{p},\,\mathsf{r}\}\mathsf{f}\left(\mathsf{a}\mathsf{x}+\mathsf{b}\mathsf{y},\,\mathsf{c}\mathsf{x}+\mathsf{c}\mathsf{b}\mathsf{y}\right)\mathsf{j} = \frac{1}{\left|\mathsf{a}+\mathsf{b}\mathsf{p}\right|}\mathsf{U}\left(\frac{\mathsf{c}+\mathsf{d}\mathsf{p}}{\mathsf{a}+\mathsf{b}\mathsf{p}},\,\mathsf{r}\frac{\mathsf{d}-\mathsf{b}\mathsf{c}+\mathsf{b}\mathsf{d}\mathsf{d}}{\mathsf{a}+\mathsf{b}\mathsf{p}}\right);\tag{\mathsf{R}.16}$$

*6. Integral along p:*

$$\mathcal{U} = \sqrt{1 + \mathcal{P}^2} \cup \langle \mathbf{p}, \mathbf{r} \rangle;\tag{8.17}$$

*7. 1D convolution expression:*

$$\mathcal{R}(\mathfrak{p}, \mathfrak{r}) \mathsf{f} f(\chi, \mathfrak{y}) \* \mathfrak{g}(\mathfrak{y}) \mathsf{f} = \mathsf{U}(\mathfrak{p}, \mathfrak{r}) \* \mathfrak{g}(\mathfrak{r});\tag{8.18}$$

*8. Equivalent equation of Plancherel theorem:*

$$\left[\int\_{-\infty}^{+\infty} \bigcup \{ \mathfrak{p}, \ \mathfrak{r} \} d\mathfrak{r} \right] = \left[\int\_{-\infty}^{+\infty} \bigwedge \{ \mathfrak{x}, \ \mathfrak{y} \} d\mathfrak{x} d\mathfrak{y} \right] \tag{\mathbb{R}.19}$$

#### *9. Equivalent equation of Parseval theorem:*

$$\int\_{-\infty}^{+\infty} \Re \langle \mathfrak{p}, \mathfrak{r} \rangle \mathbf{\hat{f}} \, f \, \langle \mathfrak{x}, \chi \rangle \mathbf{\hat{J}}^2 d\mathfrak{r} = \int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} f \, \prescript{2}{}{\mathbf{x}} \, f^2(\mathbf{x}, \chi) d\mathbf{x} d\mathfrak{y} \, \mathbf{\hat{y}} \,. \tag{8.20}$$

#### *End of Math Reminder 1.*

At this level of mathematical development a very important question must be asked: Why the switching between different coordinates systems, presented above (Fig.1), is so important to establish the mathematical basis of image projection and reconstruction in CT in the case of analytical method such as FBP? The response to this question is satisfactory explained in the 2nd demonstration (see Demonstration 2).

#### *Demonstration 2*

*Because the CT process is likely rotational in its projection phase, a rotational coordinates system (Fig.D.1) is more suitable to describe the object scanning and projection. This the first reason for using different coordinates systems.*

*Fig. D.1. Considered coordinates system: fixed (x, y) and rotational (t, s), used to identify a projection point P of the object.*

*According to figure (D.1), for a given point P of the object, the transforms used to switch between a fixed coordinates system (x,y) to a rotational coordinates system (t,s) are given by:*

$$\begin{cases} \mathfrak{t} = \times \arccos \mathfrak{t} + \mathsf{y} \cdot \mathsf{s} \mathsf{l} \mathsf{n} \Theta\\ \mathfrak{s} = -\times \mathsf{s} \mathsf{l} \mathsf{n} \mathsf{t} \theta + \mathsf{y} \cdot \mathsf{c} \mathsf{c} \mathsf{s} \mathsf{o} \mathsf{d} \end{cases} \tag{\mathsf{D}.\mathcal{T}}$$

In our case of CT, the set of all projections Pθ(t) of the object function μ (x, y) is called the Radon transform of μ(x, y). Using all these projections, a 2D image can be analytically reconstructed by exploiting a theorem called "Fourier Central Slice" (see demonstration 3). This theorem announces that the data of the 1D Fourier transform of a projection Pθ(t) is a subset of the data

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

where u = ω.cosθ and v = ω.sinθ. To demonstrate this relationship, we will follow the method presented in references [1] and [2]. Assuming that the function μ (x, y) is finite and limited, so

To demonstrate the central slice theorem of Fourier, we consider the case of θ = 0. In such a case the two coordinates systems (x, y) and (u, v) coincide and the projection Pθ(t) is simply given by:

The 2D Fourier transform of the object function becomes (by taking into account that v = ω

μ(x,y)dy}e-2πi(ux)dx

Thus, we have demonstrate that Fθ=0(ω) is just a subset of F(u, 0) for a particular case of θ=0. Because the orientation of the object in the coordinates system (u, v) is arbitrary with respect to the coordinates system (x, y), this particular case can be extended to all angles θ by keeping in mind that the Fourier transform is conserved by rotation which is necessary process in transmis‐

<sup>P</sup>*θ*=0(x)e-2πi(ux)dx

Note here that the 1D Fourier transform of a projection Pθ(t) is given by:

Fθ(ω)=*∫* −*∞* +*∞* μ(x,y)e-2πi(ux)dxdy

Pθ=0(t)= *∫*


+*∞*

FT Pθ(t) =Fθ(ω)⊂FT μ(x,y) =F(u,v), (5)

μ(x,y).e-2πi(ux+vy)dxdy. (6)

http://dx.doi.org/10.5772/52351

89

μ(x,y)dy (7)

<sup>P</sup>θ(t≡x)e-2πiωtdt (9)

(8)

of 2D Fourier transform F(u, v) of the object function μ (x, y):

F(u,v)=*∫* -*∞* +*∞ ∫* -*∞* +*∞*

> F(u,0)=*∫* -*∞* +*∞ ∫* -*∞* +*∞*

=*∫* -*∞* +*∞* {*∫* -*∞* +*∞*

=*∫* -*∞* +*∞*

that it has a Fourier transform (u, v) given by:

sin(θ = 0) = 0):

*And the inverse transforms are given by:*

{ *<sup>x</sup>* <sup>=</sup>*t*.cos<sup>θ</sup> <sup>−</sup>*s*.sin<sup>θ</sup> *<sup>y</sup>* <sup>=</sup>*t*.sin<sup>θ</sup> <sup>+</sup> *<sup>s</sup>*.cos<sup>θ</sup> (D.8)

According the physical principle of the projection generation in transmission CT, a projection at an angle θ can *mathematically be expressed as integration of the object function (attenuation coefficient) across the line s. Thus it can be given by:*

$$P\_{\theta}(t) = \left[ f \left( \times, \,\,\nu\right) \text{cls} = \right] f \left( t, \cos\theta - \text{s.s.in}\theta, \,\,t\sin\theta + \text{s.c.cos}\theta \right) \text{cls} \tag{D.9}$$

*The back-projected function fb(x,y) is then given by:*

$$\mathcal{F}\_{\clubsuit}(\times,\ \chi) = \left[ \prescript{\pi}{}{\mathcal{P}}\_{\theta}(\mathfrak{t}) \mathscr{C} \mathfrak{O} \tag{\mathsf{D}} \tag{\mathsf{D}}$$

*Where t is to be determined for each projection using equation (D.7). Thus the second reason for using rotational coordinates system is the easy representation and manipulation of Fourier transform in this system that facilitate the understanding and the implementation of the analytical reconstruction process of tomography.*

#### *End of Demonstration 2*

In our case of CT, the set of all projections Pθ(t) of the object function μ (x, y) is called the Radon transform of μ(x, y). Using all these projections, a 2D image can be analytically reconstructed by exploiting a theorem called "Fourier Central Slice" (see demonstration 3). This theorem announces that the data of the 1D Fourier transform of a projection Pθ(t) is a subset of the data of 2D Fourier transform F(u, v) of the object function μ (x, y):

*Demonstration 2*

*Because the CT process is likely rotational in its projection phase, a rotational coordinates system (Fig.D.1) is more suitable to describe the object scanning and projection. This the first reason for using different coordinates systems.*

88 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

*Fig. D.1. Considered coordinates system: fixed (x, y) and rotational (t, s), used to identify a projection point P of the*

*According to figure (D.1), for a given point P of the object, the transforms used to switch between a fixed coordinates*

*<sup>s</sup>* <sup>=</sup> <sup>−</sup> *<sup>x</sup>*.sin<sup>θ</sup> <sup>+</sup> *<sup>y</sup>*.cos<sup>θ</sup> (D.7)

*<sup>y</sup>* <sup>=</sup>*t*.sin<sup>θ</sup> <sup>+</sup> *<sup>s</sup>*.cos<sup>θ</sup> (D.8)

*P*θ(*t*)*d*θ (D.10)

*f* (*x*, *y*)*ds* =∫*f* (*t*.cosθ −*s*.sinθ, *t*.sinθ + *s*.cosθ)*ds* (D.9)

{ *<sup>t</sup>* <sup>=</sup> *<sup>x</sup>*.cos<sup>θ</sup> <sup>+</sup> *<sup>y</sup>*.sin<sup>θ</sup>

{ *<sup>x</sup>* <sup>=</sup>*t*.cos<sup>θ</sup> <sup>−</sup>*s*.sin<sup>θ</sup>

*<sup>f</sup> <sup>b</sup>*(*x*, *<sup>y</sup>*) =∫<sup>0</sup>

*understanding and the implementation of the analytical reconstruction process of tomography.*

According the physical principle of the projection generation in transmission CT, a projection at an angle θ can *mathematically be expressed as integration of the object function (attenuation coefficient) across the line s. Thus it can*

π

*Where t is to be determined for each projection using equation (D.7). Thus the second reason for using rotational coordinates system is the easy representation and manipulation of Fourier transform in this system that facilitate the*

*End of Demonstration 2*

*system (x,y) to a rotational coordinates system (t,s) are given by:*

*P*θ(*t*) =∫ *s*

*The back-projected function fb(x,y) is then given by:*

*And the inverse transforms are given by:*

*object.*

*be given by:*

$$\mathsf{FT}[\mathsf{P}\_{\theta}(\mathsf{t})] \mathsf{F}\_{\theta}(\omega) \subset \mathsf{FT}[\mathsf{\mu}(\mathsf{x}, \mathsf{y})] \mathsf{=} \mathsf{F}(\mathsf{u}, \mathsf{v}),\tag{5}$$

where u = ω.cosθ and v = ω.sinθ. To demonstrate this relationship, we will follow the method presented in references [1] and [2]. Assuming that the function μ (x, y) is finite and limited, so that it has a Fourier transform (u, v) given by:

$$\mathbf{F}(\mathbf{u}, \mathbf{v}) = \left[ \int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} \mu(\mathbf{x}, \mathbf{y}) \, \mathrm{e}^{-2\pi \mathrm{i}(\mathbf{u}\mathbf{x} + \mathbf{v}\mathbf{y})} \mathrm{d}\mathbf{x} d\mathbf{y} \, \mathrm{d}\mathbf{x} \right] \tag{6}$$

To demonstrate the central slice theorem of Fourier, we consider the case of θ = 0. In such a case the two coordinates systems (x, y) and (u, v) coincide and the projection Pθ(t) is simply given by:

$$\mathbf{P}\_{\Theta \to 0}(\mathbf{t}) = \int\_{\cdots}^{\cdots} \mu(\mathbf{x}, \mathbf{y}) \, \mathrm{d}\mathbf{y} \tag{7}$$

The 2D Fourier transform of the object function becomes (by taking into account that v = ω sin(θ = 0) = 0):

$$\begin{aligned} \mathbf{F}(\mathbf{u}, \mathbf{0}) &= \int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} \mu(\mathbf{x}, \mathbf{y}) \mathbf{e}^{-2\pi \mathbf{i}(\mathbf{x}\mathbf{x})} \mathbf{d} \mathbf{x} d\mathbf{y} \\ &= \int\_{-\infty}^{+\infty} \left\{ \int\_{-\infty}^{+\infty} \mu(\mathbf{x}, \mathbf{y}) d\mathbf{y} \right\} \mathbf{e}^{-2\pi \mathbf{i}(\mathbf{x}\mathbf{x})} d\mathbf{x} \\ &= \int\_{-\infty}^{+\infty} \mathbf{P}\_{\theta \rightarrow 0}(\mathbf{x}) \mathbf{e}^{-2\pi \mathbf{i}(\mathbf{x}\mathbf{x})} d\mathbf{x} \end{aligned} \tag{8}$$

Note here that the 1D Fourier transform of a projection Pθ(t) is given by:

$$\mathbf{F}\_{\Theta}(\omega) = \Big|\_{\ast \to \ast}^{\ast \ast \ast} \mathbf{P}\_{\Theta}(\mathbf{t} \cong \mathbf{x}) \mathbf{e}^{-2\pi \mathbf{h}\omega \mathbf{t}} \mathbf{d} \mathbf{t} \tag{9}$$

Thus, we have demonstrate that Fθ=0(ω) is just a subset of F(u, 0) for a particular case of θ=0. Because the orientation of the object in the coordinates system (u, v) is arbitrary with respect to the coordinates system (x, y), this particular case can be extended to all angles θ by keeping in mind that the Fourier transform is conserved by rotation which is necessary process in transmis‐ sion tomography in the real space. Indeed, the 1D Fourier transform of Pθ(t) produces a values which are a part (the same) of the values produced by the 2D Fourier transform of μ (x, y).

*<sup>P</sup>*θ(ω) =*TF* (*P*θ(*t*)) =∫−<sup>∞</sup>

*<sup>F</sup>* (*u*, *<sup>v</sup>*) =∫−<sup>∞</sup> +∞ ∫−∞ +∞

<sup>=</sup>∫−<sup>∞</sup> +∞ ∫0 2π ∫0 ∞

=∫0 2π ∫0 ∞

known transformation equations:

we obtain the following expression:

+∞

μ(*r*, φ)*e* <sup>−</sup>2π*i*ω*r*cos(φ−θ) |*r* |*drd*φ

Similarly, we proceed to the calculation of the 2D Fourier transform of the object function μ (x, y).

{<sup>x</sup> <sup>=</sup>*r*cos<sup>φ</sup> *y* =*r*sinφ

*F* (*u*, *v*) = *F* (ωcosψ, ωsinψ)

=∫0 2π ∫0 ∞

=∫0 2π ∫0 ∞

may allow the reconstruction of a 2D layer (slice) of the object.

Fig. D.2. Illustration of Central Slice Theorem of Fourier.

*and*{

μ(*r*, φ)*e* <sup>−</sup>2π*i*ω*r*cos(φ−ψ) |*r* |*drd*φ

The comparison between the equations (D.13) and (D.16), allow us to write the following equivalence:

*<sup>P</sup>*θ(*t*)*<sup>e</sup>* <sup>−</sup>2π*i*ω*<sup>t</sup>*

*dt*

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

(D.14)

91

(D.17)

*d*x*dy* (D.15)

http://dx.doi.org/10.5772/52351

*<sup>v</sup>* <sup>=</sup>ωsin<sup>ψ</sup> , (D.16)

δ(*r*cos(φ −θ)−*t*)μ(*r*, φ)*e* <sup>−</sup>2π*i*ω*<sup>t</sup>* |*r* |*drd*φ*dt*

μ(x, *y*)*e* <sup>−</sup>2π*i*(*u*x+*yv*)

If we transform this last expression to polar coordinates in both spatial and frequency domains using the following well-

*u* =ωcosψ

μ(*r*, φ)*e* <sup>−</sup>2π*i*ω*r*(cosφcosψ+sinφsinψ) |*r* |*drd*φ

Thus, it was demonstrated that the 1D Fourier transform of a projection is a subset (central slice) of the 2D Fourier transform of the object function. Therefore, a simple back-transformation (inverse transform) of the Fourier transform

Image reconstruction using the central slice theorem is theoretically possible for an infinite number of projections. For real data case, we have only a finite number of projections. In this case, the Fourier transform function F(u,v) is known only on number points along radial lines. In practice, the number of samples (points) taken is the same for each projection direction. As well, in the Fourier domain, the sampling is constant regardless of the direction of projection.

*T F*1*<sup>D</sup> P*θ(*t*) = *P*θ(ω) =*T F*2*<sup>D</sup>* μ(*r*, φ) = *F* (ωcosψ, ωsinψ)θ=<sup>ψ</sup> (D.18)

A dense set of values in the Fourier space can be approximated and simplified by considering square coordinates that will easily back-transformed into real space. This can, sometimes, lead to a confused reconstruction which is the consequence of an unachieved adjustment of the high frequencies in the Fourier space (Fig. 3).

**Figure 3.** Fourier transform of a projection Pθ(t) (left) and the interpolated data in the square coordinates system (right). The interpolation of high frequencies (high values of (u, v)) is inaccurate.

#### *Demonstration 3.*

*Here, we demonstrate the Fourier Central Slice Theorem announcing that the Fourier transform of a projection P*θ*(t) is a* subset of a of two-dimensional Fourier transform of the μ (x, y) object function (Fig. D.2).

*The Fourier transform of a projection P*θ(t) may be given by the following expression:

$$\mathcal{P}\_{\theta}(\omega) = \int\_{-\infty}^{\omega} \mathcal{D}\_{\theta}(t) \Theta^{-2\pi i \omega t} dt \tag{9.11}$$

The projection is given as mentioned before by:

$$P\_{\theta}(t) = \bigcap\_{-\infty \to \infty}^{+\infty \to \infty} \left\{ \mathcal{S}(\times \cos \Theta + \mathcal{Y} \sin \Theta - t) \mu(\times, \mathcal{Y}) d\mathbf{x} d\mathbf{y} \right\} \tag{D.12}$$

The last expression of Pθ(t) can be whiten in the polar coordinates system (r, φ) be the following expression:

$$\begin{aligned} \mathcal{P}\_{\theta}(t) &= \begin{bmatrix} 2\pi \\ 0 \end{bmatrix} \begin{bmatrix} \cos\varphi\cos\theta + r\sin\varphi\sin\theta - t \end{bmatrix} \mu(r,\,\varphi) \begin{bmatrix} r \mid \,\mathrm{d}r d\varphi \end{bmatrix} \\\\ &= \begin{bmatrix} 2\pi \\ 0 \end{bmatrix} \begin{bmatrix} \cos(r\cos(\varphi-\theta)-t)\mu(r,\,\varphi) \mid r \mid \,\mathrm{d}r d\varphi \end{bmatrix} \end{aligned} \tag{D.13}$$

If we take the 1D Fourier transform of last expression of Pθ(t), we get the following expression:

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples http://dx.doi.org/10.5772/52351 91

$$\begin{cases} \mathcal{P}\_{\theta}(\omega) = \text{TF}\left(\mathcal{P}\_{\theta}(t)\right) = \Big[\Big.]\_{-\infty}^{+\infty} \mathcal{P}\_{\theta}(t) \text{e}^{-2\pi i \omega t} \text{c/t} \\\\ \mathcal{P}\_{\theta} = \Big[\Big.]\_{-\infty}^{+\infty} \Big[\Big/ \int\_{0}^{\infty} \delta(r \cos(\varphi - \theta) - t) \mu(r, \cdot \varphi) \text{e}^{-2\pi i \omega t} \mid r \mid dr \, d\varphi \, d\varphi \,\text{t} \\\\ \mathcal{P}\_{\theta} = \Big[\Big/ \int\_{0}^{\infty} \mu(r, \cdot \varphi) \text{e}^{-2\pi i \omega r \cos(\varphi - \theta)} \mid r \mid dr \, d\varphi \,\text{} \end{cases} \tag{D.14}$$

Similarly, we proceed to the calculation of the 2D Fourier transform of the object function μ (x, y).

$$\mathcal{F}\{\boldsymbol{\mu},\ \mathbf{y}\} = \left[ \int\_{-\infty}^{+\infty} \right]\_{-\infty}^{+\infty} \mu(\boldsymbol{\chi},\ \mathbf{y}) \mathbf{e}^{-2\pi i(\boldsymbol{\chi}\boldsymbol{\omega} + \boldsymbol{\chi}\boldsymbol{\chi})} d\boldsymbol{\omega} d\boldsymbol{\chi} \tag{\mathbb{D}.15}$$

If we transform this last expression to polar coordinates in both spatial and frequency domains using the following wellknown transformation equations:

$$\begin{cases} \chi = r \csc \mathfrak{sp} \\ \chi = r \sin \mathfrak{sp} \end{cases} \not\equiv \omega \csc \mathfrak{sp} \tag{\Box \mathfrak{s} \Omega \mathfrak{s}} \tag{\Box \mathfrak{s} \Omega \mathfrak{s}}$$

we obtain the following expression:

sion tomography in the real space. Indeed, the 1D Fourier transform of Pθ(t) produces a values which are a part (the same) of the values produced by the 2D Fourier transform of μ (x, y).

90 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

A dense set of values in the Fourier space can be approximated and simplified by considering square coordinates that will easily back-transformed into real space. This can, sometimes, lead to a confused reconstruction which is the consequence of an unachieved adjustment of the high

Approximation

**Figure 3.** Fourier transform of a projection Pθ(t) (left) and the interpolated data in the square coordinates system

*Demonstration 3. Here, we demonstrate the Fourier Central Slice Theorem announcing that the Fourier transform of a projection P*θ*(t) is a*

*p*θ(*t*)*e* <sup>−</sup>2π*i*ω*<sup>t</sup>*

The last expression of Pθ(t) can be whiten in the polar coordinates system (r, φ) be the following expression:

δ(*r*cos(φ −θ)−*t*)μ(*r*, φ)|*r* |*drd*φ

If we take the 1D Fourier transform of last expression of Pθ(t), we get the following expression:

δ(*r*cosφcosθ + *r*sinφsinθ −*t*)μ(*r*, φ)|*r* |*drd*φ

*dt* (D.11)

(D.13)

δ(xcosθ + ysinθ −*t*)μ(x, y)*d*x*d*y (D.12)

(right). The interpolation of high frequencies (high values of (u, v)) is inaccurate.

subset of a of two-dimensional Fourier transform of the μ (x, y) object function (Fig. D.2). *The Fourier transform of a projection P*θ(t) may be given by the following expression:

> *P*θ(ω) = <sup>∫</sup> −∞

*P*θ(*t*) = ∫ −∞

*<sup>P</sup>*θ(*t*) =∫<sup>0</sup> 2π ∫0 ∞

=∫0 2π ∫0 ∞ +∞ ∫ −∞

+∞

The projection is given as mentioned before by:

+∞

frequencies in the Fourier space (Fig. 3).

$$\begin{aligned} \mathcal{F}(\boldsymbol{\omega}, \boldsymbol{\nu}) &= \mathcal{F}(\boldsymbol{\omega}\cos\varphi\boldsymbol{\upmu}, \boldsymbol{\omega}\sin\varphi\boldsymbol{\upmu}) \\ &= \begin{bmatrix} 2\pi \\ \mathbf{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{\mu} \\ \boldsymbol{\mu} \end{bmatrix} \boldsymbol{\mu}(\boldsymbol{r}, \boldsymbol{\upmu}) \mathbf{e}^{-2\pi i \boldsymbol{\upmu} \mathbf{r} (\cos\varphi\cos\varphi + \sin\varphi\sin\varphi)} \mid \boldsymbol{r} \mid \boldsymbol{\upmu}\boldsymbol{r} d\varphi \\ &= \begin{bmatrix} 2\pi \\ \mathbf{0} \end{bmatrix} \boldsymbol{\upmu}(\boldsymbol{r}, \boldsymbol{\upmu}) \mathbf{e}^{-2\pi i \boldsymbol{\upmu} \mathbf{r} \cos(\boldsymbol{\upmu}-\boldsymbol{\upmu})} \mid \boldsymbol{r} \mid \boldsymbol{\upmu}\boldsymbol{r} d\varphi \end{bmatrix} \end{aligned} \tag{D.17}$$

The comparison between the equations (D.13) and (D.16), allow us to write the following equivalence:

$$\mathsf{T}\,\mathsf{F}\_{1\mathfrak{L}}\mathsf{I}\,\mathsf{P}\_{\theta}(t)\Big| = \mathsf{P}\_{\theta}(\omega) = \mathsf{T}\,\mathsf{F}\_{2\mathfrak{L}}\mathsf{I}\,\mu(r,\,\mathsf{q}\,\mathsf{I})\Big| = \mathsf{F}\,(\omega\circledcirc\mathsf{cos}\,\mathsf{y},\,\omega\sin\mathsf{i}\,\mathsf{n}\,\mathsf{y})\_{\theta\sim\mathsf{q}}\tag{\mathsf{D}.18}$$

Thus, it was demonstrated that the 1D Fourier transform of a projection is a subset (central slice) of the 2D Fourier transform of the object function. Therefore, a simple back-transformation (inverse transform) of the Fourier transform may allow the reconstruction of a 2D layer (slice) of the object.

#### Fig. D.2. Illustration of Central Slice Theorem of Fourier.

Image reconstruction using the central slice theorem is theoretically possible for an infinite number of projections. For real data case, we have only a finite number of projections. In this case, the Fourier transform function F(u,v) is known only on number points along radial lines. In practice, the number of samples (points) taken is the same for each projection direction. As well, in the Fourier domain, the sampling is constant regardless of the direction of projection.

*Thus, the digitalization and sampling is a very important process in practical tomography. To be able to reconstruct the object function, these samples (points) must be interpolated from a polar to a Cartesian coordinates (Fig. D.3). Generally, this interpolation is done by taking the nearest neighbour value or by a linear interpolation between known points. The density of points (frequencies) in the polar reference becomes smaller when we move away from low frequencies (i.e., the origin). So the interpolation error is larger at high frequencies than at low frequencies, and this causes the degradation of the image details.*

<sup>D</sup> P= Δy p

<sup>π</sup> P S 2

Indeed, S and P are the most important parameters that determine the quality of the recon‐ struction. If P violates the last condition (Eq.11) when for example a number less than necessary P projections is recorded, a new reduced diameter (D = D \*) which satisfies the Shannon condition must be calculated. This reduced dimater determines the reduced volume (area) of the object that can be well reconstructed for this reduced number of projections although the

*Demonstration 4.*

*In computed tomography, it is recommended that the number of projections (P) must be in the same order as the number of rows of pixels in a single projection (S) [3]. The condition on P and s given by Eq.12 which is established when*

*considering the Nyquist-Shannon sampling theorem can be proven by the following manner. Let considering P* projections over 180° and S rows of pixels, the angular increment between two successive projections Δθ in the Fourier

Δθ= <sup>π</sup>

<sup>ω</sup>max<sup>=</sup> <sup>1</sup>

Δf=ωmaxΔθ= <sup>1</sup>

For a distance Δx between two adjacent rows, the highest spatial frequency measured (ω*max) in a projection is given*

The scanned object representation in the frequency domain (2D Fourier Transform) correspond to disk's radius that contains the measured values (Fig. D. 4). The distance Δf between two consecutive values of the furthest circle from the

> 2Δx π

³ (12)

<sup>P</sup> (D.18)

<sup>2</sup>Δ<sup>x</sup> (D.19)

<sup>P</sup> (D.20)

demonstration 4):

*space is given by [4]:*

*origin is given by:*

whole object of real diameter D is scanned.

*according to the Nyquist-Shannon theorem by [5]:*

According to Nyquist-Shannon's theorem it is required for a good reconstruction that Δy≥ 2Δx. If we put Δy ≥ 2Δx, we obtain a relationship between the number of scanned points S in one projection line (sampling of the object) and the number of necessary projections P (see

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

(11)

93

http://dx.doi.org/10.5772/52351

*Fig. D.3. Transition from a polar grid to a Cartesian square grid.*

#### *End of Demonstration 3*

How many projections are needed for good reconstruction? This question has been answered by the Nyquist-Shannon theorem. This theorem announced that a unique reconstruction of an object sampled in space is obtained if the object was sampled with a frequency greater than twice the highest frequency of the object details. In the parallel scan mode, for a sampling for S points per projection line, a number of P projections (angles θ) is necessary to accomplish the Shannon theorem in tomography. If D is the diameter of the object to be scanned, and Ax is the difference between two points of scanning, then the number of points (sampling) in each projection line to be scanned is given by:

$$\mathbf{S} = \frac{\mathbf{D}}{\Delta \mathbf{x}} \tag{10}$$

For a scanning of the object over 360 °1 , each point is scanned again after a path equal to πD and this for each point situated on the surface of a circular object of a diameter D (the highest frequency). For this case, the number of projections must be equal to:

<sup>1</sup> Generally if the object is homogenous and symmetric, a scanning over 180° is sufficient.

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples http://dx.doi.org/10.5772/52351 93

$$\mathbf{P} = \frac{\pi \mathbf{D}}{\Delta \mathbf{y}} \tag{11}$$

According to Nyquist-Shannon's theorem it is required for a good reconstruction that Δy≥ 2Δx. If we put Δy ≥ 2Δx, we obtain a relationship between the number of scanned points S in one projection line (sampling of the object) and the number of necessary projections P (see demonstration 4):

*Thus, the digitalization and sampling is a very important process in practical tomography. To be able to reconstruct the object function, these samples (points) must be interpolated from a polar to a Cartesian coordinates (Fig. D.3). Generally, this interpolation is done by taking the nearest neighbour value or by a linear interpolation between known points. The density of points (frequencies) in the polar reference becomes smaller when we move away from low frequencies (i.e., the origin). So the interpolation error is larger at high frequencies than at low frequencies, and this causes the*

92 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

*End of Demonstration 3*

How many projections are needed for good reconstruction? This question has been answered by the Nyquist-Shannon theorem. This theorem announced that a unique reconstruction of an object sampled in space is obtained if the object was sampled with a frequency greater than twice the highest frequency of the object details. In the parallel scan mode, for a sampling for S points per projection line, a number of P projections (angles θ) is necessary to accomplish the Shannon theorem in tomography. If D is the diameter of the object to be scanned, and Ax is the difference between two points of scanning, then the number of points (sampling) in each

<sup>D</sup> S=

frequency). For this case, the number of projections must be equal to:

1 Generally if the object is homogenous and symmetric, a scanning over 180° is sufficient.

and this for each point situated on the surface of a circular object of a diameter D (the highest

Δx (10)

, each point is scanned again after a path equal to πD

*degradation of the image details.*

*Fig. D.3. Transition from a polar grid to a Cartesian square grid.*

projection line to be scanned is given by:

For a scanning of the object over 360 °1

$$\mathbf{P} \ge \frac{\pi}{2}\mathbf{S} \tag{12}$$

Indeed, S and P are the most important parameters that determine the quality of the recon‐ struction. If P violates the last condition (Eq.11) when for example a number less than necessary P projections is recorded, a new reduced diameter (D = D \*) which satisfies the Shannon condition must be calculated. This reduced dimater determines the reduced volume (area) of the object that can be well reconstructed for this reduced number of projections although the whole object of real diameter D is scanned.

#### *Demonstration 4.*

*In computed tomography, it is recommended that the number of projections (P) must be in the same order as the number of rows of pixels in a single projection (S) [3]. The condition on P and s given by Eq.12 which is established when considering the Nyquist-Shannon sampling theorem can be proven by the following manner. Let considering P* projections over 180° and S rows of pixels, the angular increment between two successive projections Δθ in the Fourier *space is given by [4]:*

$$
\Delta\Theta = \frac{\Pi}{\overline{\mathsf{P}}} \tag{D.18}
$$

For a distance Δx between two adjacent rows, the highest spatial frequency measured (ω*max) in a projection is given according to the Nyquist-Shannon theorem by [5]:*

$$
\omega\_{\text{max}} = \frac{1}{\overline{\Delta \Delta \chi}} \tag{D.19}
$$

The scanned object representation in the frequency domain (2D Fourier Transform) correspond to disk's radius that contains the measured values (Fig. D. 4). The distance Δf between two consecutive values of the furthest circle from the *origin is given by:*

$$
\Delta \mathbf{f} = \omega\_{\text{max}} \Delta \mathbf{f} = \frac{1}{2\Delta \mathbf{x}} \frac{\mathbf{r}}{\mathbf{F}} \tag{\text{D.20}}
$$

*Fig. D.4. Density of the measured values in the frequency domain.*

*For the S values of each projection in spatial domain correspond S measured values in the frequency space (measured* for each line). Therefore, the distance ε between two consecutive values measured on a radial line in the frequency domain (Fig. D.4) is given by:

$$\varepsilon = \frac{2\,\omega\_{\text{max}}}{\mathcal{S}} = \frac{1}{\overline{\mathcal{D}}\,\mathcal{S}}\tag{\text{D.2.1}}$$

*Fig. D.5. Results of 2D reconstruction of an object by FBP method of tomography and induced artifacts for different number of projections: (A) original image, (B) 1 projection, (C) 3 projections, (D) 4 projections, (E) 16 projections, (F) 32*

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

http://dx.doi.org/10.5772/52351

95

*The number of projections is not the only parameter that affects the reconstructed image. The sampling of projections S (number of rows) and the dimensions of the grid of reconstruction will also affect the reconstructed image dramatically*

*End of Demonstration 4*

As already mentioned, the inverse transformation (simple inversion) of the Fourier transform Pθ(ω) can be directly used to produce the 2D reconstructed layer of μ (x, y). However, a more efficient way and more elegant has been developed called " Filtered Back-Projection (FBP)" making the reconstruction process less expensive and less complicated when compared to the direct calculation of the inverse Fourier transform [1,2]. The basic idea of this method is derived from the tomography principle and the scanning mode: the object (set of layers) is scanned projection by projection from 0 ° to 180 ° which means that Pθ(t) = Pθ+180(-t) and suggesting the use of a polar coordinates rather than Cartesian square ones. As it was demonstrated the digitization and sampling details by the right selection of P and S are very important in Computed tomography. If the projections are recorded in polar coordinates, the projections Fourier transforms will be discrete values of a polar function. Thus, it seems appropriate to write the 2D Fourier transform of the object function μ (x, y) in polar coordinates. The inverse

Fourier transform μ (x, y) of F(u, v) written in Cartesian coordinates is given by:

The same inverse transform in polar coordinates is given by:

p

allow us to get the following expression:

+ + +2πi(ux+vy)

2 +2πiω(xcosθ+ysinθ) <sup>0</sup> μ(x,y)= F(ω,θ).e d dθ

The substitution of x cos (θ) + y sin (θ) by t and the application of Fourier Central Slice theorem,


w w+¥ ò ò-¥ (14)

*projections and (G) 64 projections [7] .*

*1.1.2. 2D and 3D image reconstruction*

*if they are not well optimized*

A sufficient condition to obtain a good reconstruction is to ensure that the worst azimuth resolution (s) in the frequency domain (Fig. D.4) is in the same order of radial resolution (ε). This condition can be expressed as follows:

$$\frac{1}{\text{L2D}} \frac{\text{m}}{\text{P}} \approx \frac{1}{\text{D} \text{S}} \implies \text{P} \approx \text{S} \frac{\text{m}}{\text{Z}} \tag{\text{D}}$$

Thus, the ratio between the number of projections (P) and the number of rows (S) must be in the order of π / 2.

An insufficient number of projections may produce a 2D or 3D reconstructed image which present many undesirable artefacts. In practice, most of the tomography detectors cannot measure below the nominal Nyquist resolution determined by their dimensions or pixel size. The exploration beam is not perfectly parallel, too; this is another source of errors on the measured data obtained especially when using a reconstruction method which assumes that the beam is *parallel.*

In practice, insufficient number of projection generates artefacts such as those shown in figure (D.5) [3]. The projections of figure (D.5) have a dimension of 64x64 pixels. So the number of rows S is 64. As it has already been demonstrated, a number P of projections around a value of 100 (P ≈ Sπ / 2) is sufficient to produce an acceptable reconstructed image. For the case of this example, it is clear that the reconstruction image obtained from 64 projections is the closest one to the original image. The 2D reconstruction is performed using the Filtered Back-Projection method (FBP) that will be described in the next section.

*Fig. D.5. Results of 2D reconstruction of an object by FBP method of tomography and induced artifacts for different number of projections: (A) original image, (B) 1 projection, (C) 3 projections, (D) 4 projections, (E) 16 projections, (F) 32 projections and (G) 64 projections [7] .*

*The number of projections is not the only parameter that affects the reconstructed image. The sampling of projections S (number of rows) and the dimensions of the grid of reconstruction will also affect the reconstructed image dramatically if they are not well optimized*

#### *End of Demonstration 4*
