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

1. Introduction

The visible light being reflected by the majority of the objects which surround us, we can apprehend our environment only by the properties of the surface of the objects which compose it. To exceed this limit and to explore the intimacy of the matter, a dedicated imaging techniques and instruments, which used penetrating radiations like X-rays, neutrons, gamma, or certain electromagnetic or acoustic rays to explore internal structure, are developed. The tomography is one of these developed techniques that allow 2D and 3D interior object examination. By combination of a set of measures and thanks to computational and images reconstruction methods, it provides cartography of attenuation parameters characteristic of the radiation/object interaction, according to one or more transversal plans (slices). It thus makes possible to see on TV monitor the interior of bodies and objects, whereas before one had access either by pure imagination, by interpreting indirect measurements, or by cutting out the objects materially. In the case of the medical imaging, for example, this direct observation requires a surgical operation. This formidable invention thus enables us to discover the interior of human bodies and different objects which surround us and their organization in space and time, without destroying them.

The tomography thus constitutes an instrument privileged to analyze and characterize matter, that it is inert or alive, static or dynamic, of microscopic or macroscopic scale. While giving access to the structure and the shape of the components, it makes possible the apprehension of the complexity of the objects studied. The computed tomography is a technique of acquisition of digital images (projections). It generates a computer coding of a digital representation of an area of interest through a patient, a structure or an object. The tomography thus provides a virtual representation of reality, in the intimacy of its composition. This numerical coding then will facilitate the exploitation, the exchanges and the information 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 position, 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 processing, analysis and visualization.

A CT system gathers several technological components. Its development requires the participation 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 industrial teams to develop, produce and market these systems.

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 ambiguous or not clear are explained and mathematically demonstrated with some illustrative examples. 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 methods for image reconstruction in tomography will be also outlined.

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 radiation transmission radiography but also for novices in the field of digital imaging. The chapter 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 presented 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 simply a curiosity about the technique.

1.1. Analytical methods for image reconstruction in tomography

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 reconstruction 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 reconstruction 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.

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 f_n(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 establishment 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 (comb-shaped) 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 scanned 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(x,y)=I_0\exp (-{\displaystyle \underset{path}{\int }\mu (x,y)ds)}$$E1

Where I₀ 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=x.\cos (\theta )+y.\sin (\theta )$$E2

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}.

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:

$${\text{P}}_{\text{θ}}\text{(t)=ln(}\raisebox{1ex}{$\text{I}$}\!\left/ \!\raisebox{-1ex}{${\text{I}}_{\text{0}}$}\right.\text{)=}{\displaystyle \underset{\text{path}}{\int }\text{μ(x,y)ds}}$$E3

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

$${\text{P}}_{\text{θ}}\text{(t)=}{\displaystyle \underset{\text{path}}{\int }\text{μ(x,y)ds}}$$

$$={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\delta (x\cos \theta +y\sin \theta -t)\mu (x,y)dxdy}$$E4

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_θ (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:
$$P_{\theta }(t)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\delta (x\cos \theta +y\sin \theta -t)\mu (x,y)dxdy}}$$	(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{array}{c}t=x.\cos \theta +y.\sin \theta \\ s=-x.\sin \theta +y\cos \theta \end{array}$$	(D.2)
With:
$$dsdt=J_{F}dxdy$$	(D.3)
JF is the Jacobian matrix determinant which is given by:
$$J_F=\left[\begin{array}{c}\begin{array}{ccc}& \raisebox{1ex}{$\partial $}\!\left/ \!\raisebox{-1ex}{$\partial x$}\right.& \raisebox{1ex}{$\partial $}\!\left/ \!\raisebox{-1ex}{$\partial y$}\right.\end{array}\\ \begin{array}{ccc}\begin{array}{c}t\\ s\end{array}& \begin{array}{c}\cos \theta \\ -\sin \theta \end{array}& \begin{array}{c}\sin \theta \\ \cos \theta \end{array}\end{array}\end{array}\right]={\cos }^2\theta +{\sin }^2\theta =1$$	(D.4)
So, in our case dsdt=dxdy, and we can write the following equation:
$$\begin{array}{l}{P}_{\theta }(t)={\displaystyle \underset{ray(\theta ,t)}{\int }\mu (x,y)ds=}{\displaystyle \underset{ray(\theta ,t)}{\int }\mu (s,t)ds}\\ ={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\delta (x\cos \theta +y\sin \theta -t)\mu (s,t)dtds}}\end{array}$$	(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{array}{l}R(p,\tau )\left[f(x,y)\right]={\displaystyle {\int }_{-\infty }^{+\infty }f(x+\tau +p}x)dx\\ ={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }f(x,y)\delta \left[y-(\tau +px)\right]}}dxdy\equiv U(p,\tau ),\end{array}$$	(R.1)
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:
$$f(x,y)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }\frac{d}{dy}H\left[U(p,y-px)\right]}dp,$$	(R.2)
where H is the Hilbert transform. The Radon transform can also be defined by the following expression:
$$R^{\prime }(r,\alpha )\left[f(x,y)\right]={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }f(x,y)\delta (r-x\cos \alpha -y\sin \alpha )dxdy}},$$	(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:
$$F_{\omega ,\alpha }\left[R\left[f(\omega ,\alpha )\right]\right](x,y)=F_{u,v}^2\left[f(u,v)\right](x,y),$$	(R.4)
where F is the Fourier transform, the Radon inversion formula can be expressed by:
$$f(x,y)=c{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\infty }^{+\infty }{F}_{\omega ,\alpha }\left[R\left[f(\omega ,\alpha \right]\right]}}\left|\omega \right|e^{i\omega (x\cos \alpha +y\sin \alpha )}d\omega d\alpha$$	(R.5)
This last expression can be simplified as follows:
$$f(x,y)={\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\infty }^{+\infty }R\left[f(r,\alpha )\right]}}W(r,\alpha ,x,y)drd\alpha ,$$	(R.6)
where W is a weight function given by:
$$W(r,\alpha ,x,y)=h(x\cos \alpha +y\sin \alpha -r)=F^{-1}\left[\left|\omega \right|\right]$$	(R.7)
Nievergelt (1986) determined the inversion formula as follows:
$$f(x,y)=\frac{1}{\pi }\underset{c\to 0}{\lim }{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\infty }^{+\infty }R\left[f(r+x\cos \alpha +y\sin \alpha ,\alpha )\right]}}{G}_c(r)drd\alpha ,$$	(R.8)
with:
$$G_c(r)=\{\begin{array}{cc}\frac{1}{\pi c^2}& for\left|r\right|\le c\\ \frac{1}{\pi c^2}(1-\frac{1}{\sqrt{1-c^2/r^2}}& for\left|r\right|\succ c\end{array}$$	(R.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:
$$\begin{array}{cc}p=\cot \alpha ,& \tau =\csc \alpha \end{array}$$	(R.10)
$$\begin{array}{cc}r=\frac{\tau }{1+p^2},& \alpha ={\cot }^{-1}p\end{array}$$	(R.11)
The main properties of Radon transform are the followings:
1. Superposition :
$$R(p,\tau )\left[f_1(x,y)+f_2(x,y)\right]=U_1(p,\tau )+U_2(p,\tau );$$	(R.12)
2. Linearity:
$$R(p,\tau )\left[af(x,y)\right]=aU(p,\tau );$$	(R.13)
3. Scaling:
$$R(p,\tau )\left[f_1(\frac{x}{a},\frac{y}{b})\right]=\left|a\right|U(p\frac{a}{b},\frac{\tau }{b});$$	(R.14)
4. Rotation:
$$R(p,\tau )\left[R_{\Phi }f(x,y)\right]=\frac{1}{\left|\cos \Phi +p\sin \Phi \right|}U(\frac{p-\tan \Phi }{1+p\tan \Phi },\frac{\tau }{\cos \Phi +p\sin \Phi }),$$	(R.15)
RΦ is a rotation operator;
5. Skewing:
$$R(p,\tau )\left[f(ax+by,cx+dy)\right]=\frac{1}{\left|a+bp\right|}U(\frac{c+dp}{a+bp},\tau \frac{d-b(c+bd)}{a+bp});$$	(R.16)
6. Integral along p:
$$I=\sqrt{1+p^2}U(p,\tau );$$	(R.17)
7. 1D convolution expression:
$$R(p,\tau )\left[f(x,y)\ast g(y)\right]=U(p,\tau )\ast g(\tau );$$	(R.18)
8. Equivalent equation of Plancherel theorem:
$${\displaystyle {\int }_{-\infty }^{+\infty }U(p,\tau )d\tau ={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }f(x,y)dxdy}}};$$	(R.19)
9. Equivalent equation of Parseval theorem:
$${\displaystyle {\int }_{-\infty }^{+\infty }R(p,\tau ){\left[f(x,y)\right]}^{2}d\tau ={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }{f}^2(x,y)dxdy}}};$$	(R.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 2^nd 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{array}{c}t=x.\cos \theta +y.\sin \theta \\ s=-x.\sin \theta +y.\cos \theta \end{array}$$	(D.7)
And the inverse transforms are given by:
$$\{\begin{array}{c}x=t.\cos \theta -s.\sin \theta \\ y=t.\sin \theta +s.\cos \theta \end{array}$$	(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)={\displaystyle \underset{s}{\int }f(x,y)ds}={\displaystyle \int f(t.\cos \theta -s.\sin \theta ,t.\sin \theta +s.\cos \theta })ds$$	(D.9)
The back-projected function fb(x,y) is then given by:
$$f_b(x,y)={\displaystyle {\int }_0^{\pi }{P}_{\theta }(t)d\theta }$$	(D.10)
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):

$$\text{FT}\left[{\text{P}}_{\text{θ}}\text{(t)}\right]{\text{=F}}_{\text{θ}}\text{(ω)}\subset \text{FT}\left[\text{μ(x,y)}\right]\text{=F(u,v),}$$E5

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:

$$\text{F(u,v)=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\displaystyle {\int }_{\text{-}\infty }^{+\infty }\text{μ(x,y)}{\text{.e}}^{\text{-2πi(ux+vy)}}\text{dxdy}}}\text{.}$$E6

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:

$${\text{P}}_{\text{θ=0}}\text{(t)=}{\displaystyle \underset{\text{-}\infty }{\overset{\text{+}\infty }{\int }}\text{μ(x,y)dy}}$$E7

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

$$\begin{array}{l}\text{F(u,0)=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\text{μ(x,y)e}}^{\text{-2πi(ux)}}\text{dxdy}}}\\ \text{=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }\{{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }\text{μ(x,y)}}\text{dy}\}{\text{e}}^{\text{-2πi(ux)}}\text{dx}}\\ \text{=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\text{P}}_{\theta =0}{\text{(x)e}}^{\text{-2πi(ux)}}\text{dx}}\end{array}$$E8

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

$${\text{F}}_{\text{θ}}\text{(ω)=}{\displaystyle {\int }_{-\infty }^{+\infty }{\text{P}}_{\text{θ}}\text{(t}\equiv {\text{x)e}}^{\text{-2πiωt}}\text{dt}}$$E9

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 transmission 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).

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.

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:

$$\text{S=}\frac{\text{D}}{\text{Δx}}$$E10

For a scanning of the object over 360 °

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

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

$$\text{P=}\frac{\pi \text{D}}{\text{Δy}}$$E11

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

$$\text{P}\ge \frac{\text{π}}{\text{2}}\text{S}$$E12

Indeed, S and P are the most important parameters that determine the quality of the reconstruction. 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]:
$$\text{Δθ=}\frac{\text{π}}{\text{P}}$$	(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]:
$${\text{ω}}_{\text{max}}\text{=}\frac{\text{1}}{\text{2}\Delta \text{x}}$$	(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:
$${\text{Δf=ω}}_{\text{max}}\text{Δθ=}\frac{\text{1}}{\text{2Δx}}\frac{\text{π}}{\text{P}}$$	(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:
$$\text{ε=}\frac{{\text{2ω}}_{\text{max}}}{\text{S}}\text{=}\frac{\text{1}}{\text{DS}}$$	(D.21)
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{\text{1}}{\text{2D}}\frac{\text{π}}{\text{P}}\approx \frac{\text{1}}{\text{DS}}\Rightarrow \text{P}\approx \text{S}\frac{\text{π}}{\text{2}}$$	(D.22)
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

1.1.2. 2D and 3D image reconstruction

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:

$$\text{μ(x,y)=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }\text{F(u,v)}{\text{.e}}^{\text{+2πi(ux+vy)}}\text{dudv}}}$$E13

The same inverse transform in polar coordinates is given by:

$$\text{μ(x,y)=}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }\text{F(ω,θ)}{\text{.e}}^{\text{+2πiω(xcosθ+ysinθ)}}\omega \text{d}\omega \text{dθ}}}$$E14

The substitution of x cos (θ) + y sin (θ) by t and the application of Fourier Central Slice theorem, allow us to get the following expression:

$$\begin{array}{l}\text{μ(x,y)=}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\infty }^{+\infty }\left[\text{F(ω,θ)}{\text{.e}}^{\text{+2πiωt}}\left|\omega \right|\text{d}\omega \right]\text{dθ}}}\\ \text{=}{\displaystyle {\int }_0^{\pi }\left[{\displaystyle {\int }_{-\infty }^{+\infty }\text{F(ω,θ)}{\text{.e}}^{\text{+2πiωt}}\left|\omega \right|\text{d}\omega }\right]}\text{dθ}\end{array}$$E15

The integral in brackets can be regarded as the inverse Fourier transform P_θ(t) of P_θ(ω). However, we note that it is multiplied by the function |ω| which plays the role of a special ramp filter function in the frequency domain. So, we define the filtered projection by:

$${\text{Q}}_{\text{θ}}\text{(t)=}{\displaystyle {\int }_{\text{-}\infty }^{\text{+}\infty }{\text{P}}_{\text{θ}}\text{(ω)}\left|\text{ω}\right|}{\text{e}}^{\text{+2πiωt}}\text{dω}$$E16

Thus, the object function will be simply given by:

$$\text{μ(x,y)=}{\displaystyle {\int }_0^{\pi }\left[Q_{\theta }(t)\right]}\text{dθ=}{\displaystyle {\int }_0^{\pi }\left[{\text{Q}}_{\text{θ}}\text{(xcosθ+ysinθ)}\right]}\text{dθ}$$E17

Knowing previously that a product (multiplication) in Fourier space (frequency domain) corresponds to a convolution of inverse Fourier transforms in real space (spatial domain), the following relation can be written:

$${\text{FT}}^{\text{-1}}\left[{\text{P}}_{\text{θ}}\text{(ω)}\times \left|\text{ω}\right|\right]{\text{=P}}_{\text{θ}}\text{(t)}\otimes {\text{FT}}^{\text{-1}}\left[\left|\text{ω}\right|\right]$$E18

The convolution operator is denoted by the symbol ⊗. The function |ω| is not a square-integrable function, thus it has no inverse Fourier transform. However, the inverse transform FT^-1[|ω|] can be approximated by different filter response functions (convolution with gains, filter functions). Considering these last remarks and regarding the result of Eq.17, The reconstruction process of μ (x, y) can be performed from the projections P_θ(t) obtained as follows:

“In FBP reconstruction process, each projection P_θ(t) is first convoluted with a specific and suitable filtering function (eg Shepp-Logan, ω/2π sinc(ω)) and the measured values are recorded in (x, y) plans as illustrated in Fig. 4. To control the measurement process of tomography, the projections are arranged in a sinograms

the number of sinograms is equal to the sampling number of S

recorded as shown in Fig. 5. On the sinogram, each vertical line is a projection at an angle θ and represents the variation in the gray level as a function of the pixel’s or detector’s position. On Fig. 6 are shown 4 different reconstructions. The projections P_θ(t) are converted to grayscale, convoluted with a filter with a specific gain and back-projected onto the entire plan (x, y). The addition of all projections result in the reconstruction of the 2D layer desired. More the number of projections is high more the reconstruction plans is well covered in terms of data and less the star-shaped artefacts are present on the reconstructed image. With such reconstructions, 3D images can be obtained by stacking all the 2D layers and a reconstructed 3D volume data (details of the object) can be extracted from the stack obtained”.

Figure 4.

Scanning a 2D object (details: circle and square) and the corresponding projections [7].

Figure 5.

Example of projection Sinogram of the object of Fig. 4: the y-axis shows the projections (gray level variation as a function of the pixel’s or detector’s position) and on the x-axis, these one pixel width projections are arranged from the first projection (θ=0) to the last one (θ=180°). Each sinogram will be used for the reconstruction of one pixel layer of the object.

Figure 6.

Reconstructions obtained for 4 different numbers of projections to illustrate the appearance of star-shaped artefacts when the number of projection is not sufficient.

1.1.3. Discretization of the analytical methods

In the last sections, many times we have considered the ideal continuous case to model the projection and the reconstruction processes of tomography and to explain the principle of the analytical reconstruction methods such as FBP. Indeed, we have worked in an infinite continuous domain (R²), and tried to reconstruct a continuous function f from continuous projections P_θ defined for all angle θ in the interval [0, π [. These conditions are obviously not possible in practice. Acquisition systems allow just obtaining a number of projections P_θ for a finite number of angles, denoted θ_k. The limited number of detectors makes these projections sampled and known only at discrete points u_k. It would be unrealistic with such data to try to reconstruct a continuous function f in R². Therefore, the object function f(x,y) will be reconstructed on a discrete grid for a finite number of points (f(x_i,y_j)). This is also the limit imposed by the used numerical reconstruction algorithms. The reconstruction problem then will be approached as follows: for given a set of projection measurements:

$$\{{\text{P}}_{\text{θ}}{}_{\text{k}}{\text{(u}}_{\text{l}}\text{),0}\le \text{l}\le \text{S,0}\le \text{k}\le \text{P}\}\text{,}$$E19

the reconstruction problem is reduced to finding f at any point of a finite discrete grid:

$$\{{\text{f(x}}_{\text{i}}{\text{,y}}_{\text{j}}\text{)}\text{,0}\le \text{i}\le \text{S,0}\le \text{j}\le \text{S}\}$$E20

with:

$$\begin{array}{ccc}\begin{array}{c}{\theta }_k=k\Delta \theta ,\\ x_i=i\Delta x,\end{array}& \begin{array}{c}\Delta \theta =\pi /P,\\ y_j=j\Delta y,\end{array}& u_l=ld,\end{array}$$E21

where Δθ is the sampling step of the rotation angles, P is the number of angles (projections), d is the sampling step on each projection line (ray), Δx and Δy are the sampling step on x and y in the reconstruction plan.

The discrete reconstruction methods according to this definition can be divided into two categories:

1. The first class method consists in the definition of discrete operators and functions equivalent to those defined in the continuous case (Radon transform, back-projection, Fourier transform, etc.) and all the inversion formulas and theorems defined for the analytical methods already presented.

2. The second category is based on a completely different approach: in these methods the projection equation (p=Rf) is directly discredited and a linear equations system is built. The resolution of this system is possible only by iterative methods called algebraic methods. The algebraic method will be presented in the next section.

2. Algebraic methods of image reconstruction in tomography

The algebraic methods are especially iterative methods. These methods are less used when compared to the easy and well known filtered back-projection method. In this method the reconstruction problem is see differently and no longer refers to the Radon transform. The image of the object consists of a number k of pixels whose values f_k are unknown. Similarly, the projections are discrete and formed by a number of l dexels (depth pixels) whose values "p_l" are known since they correspond to measurements in each line of projection. Reconstruction of the object image by the iterative method is based on the following hypothesis: each detected values in a dexel is a linear combination of pixel’s values to be reconstructed [6]. The reconstruction problem is formulated by a discrete expression of matrix (p = R.f) describing the projection process. The set of values of the projection lines (dexels) is arranged in projection vector p. All pixels of the image to be reconstructed are also grouped in a image vector f. The coefficients which characterize the contribution of each pixel in each line of projection is determined and stored in a matrix R. The projection of an object is modeled in the case of an original image of n x n pixels, a number n of projection directions and a number n of dexels in each projection line, by the following equation [6]:

$$\text{p=R}\text{.f}$$E22

$$\left[\begin{array}{c}{p}_1\\ p_2\\ ⋮\\ ⋮\\ p_n\end{array}\right]=\left[\begin{array}{ccccc}{r}_{11}& \dots & \dots & \dots & r_{14}\\ ⋮& \ddots & & & \\ ⋮& & \ddots & & \\ ⋮& & & \ddots & \\ r_{41}& \dots & \dots & \dots & r_{44}\end{array}\right].\left[\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ ⋮\\ f_n\end{array}\right]$$E23

This last equation expresses the fact that what is detected (p) is the result of values (f) of the image to be reconstructed, subject to a projection operation represented by the projection operator R [6]). Through this modeling of the projection process, one looks in practice to find f according to p by solving the inverse problem f = R^-1.p. Because of the size of this equations system, the resolution cannot be performed except by successive iterations [7]. Fig. 7 shows an example of projection and back-projection. In the algebraic reconstruction method based on iterative process, the back-projection is modelled by a back-projection operator R^t who is none other than the transposed matrix of R. Thus, the reconstruction problem is limited in solving the inverse problem f = R^t.p.

The resolution of the inverse problem by iterative methods consists in finding a solution f minimizing the distance d between p and R.f where p and R are known. Here, it is question to start from an arbitrary estimate of the image solution and to proceed schematically to the correction the first estimate basing on a principle of trial and error. Each next estimate is projected and the result obtained is compared to the measured projection. The returned error is used to improve the next estimate. This method leads to build gradually low-to-high frequencies of the image solution. The results of the first iterations are smooth because of the predominance of low frequencies (internal structure) of the object. Subsequently, more iteration are applied more high frequencies (overall shape and background noise) are represented. The image produced by iterations approximate gradually the image solution (the algorithm converges) [6]. However, we show that when using such iterative methods and after a number of iterations, the process begins to diverge (under the influence of noise) and the image moves away from the true solution. To overcome this inconvenient, we impose a constraint to the reconstruction process to interrupt the iteration after a certain number of iterations. This is equivalent to use a low pass filter as in the case of filtered back-projection method.

Figure 7.

Illustration of projection and back-projection processes.

2.1. Projection modeling

The projection process is modeled by considering the coefficients of the projection matrix R that generate the acquisition data. Some special geometric and physical considerations are necessary for this modeling. In the iterative reconstruction method the modeling concerns the following points:

1. Modeling of the pixel (detector) intensity distribution: it is necessary to specify the real conditions of image pixels projection [6]. It is based on the evaluation of the contribution of each pixel in the corresponding projection line (ray). Knowing that the most accurate model (perfect) and the more complicated to be applied consists in considering square pixel (uniform); simpler models are possible. Among these simple models, there is the model called Dirac model in which all the pixel intensity is concentrated in the centre of the pixel. Thus, the whole intensity of the pixel contributes to the projection line (ray) if and only if it passes through the Dexel. There is also another model called the model of concave disc which is considered as a compromise between the two previous models. In this model, the intensity of the pixel is geometrically limited to a disk included in the pixel and distributed so that its projection is rectangular regardless of the direction of projection [6].

2. Geometric modeling of the projection operator: for the determination of the coefficients of the projection matrix R, we must consider the number of projections and their angular distributions and the projection beam geometry which can be parallel or fan. If, for example, the intensity model distribution is that of Dirac, a given pixel with an index k crossed by a ray of an index p_l generates a projection coefficient r_lk equal to 1; if this is not the case, zero value is assigned.

3. Physical modeling of the projection operator: this model is based on the distance between the position of the object pixel to the detector. A pixel located away from the detector will see its contribution to the projection ray reduced compared to a nearest one. With this modeling aspect the beam attenuation will be included in the resulting reconstructed image.

2.2. Main iterative algorithms

There are many algorithms that have been developed that used the iterative method for tomographic reconstruction (see practical example 1). Among these algorithms, the main ones are: the Algebraic Reconstruction Technique (ART), the Simultaneous Iterative Reconstruction Technique (SIRT), and the Iterative Least Squared Technique (ILST). Currently, the most commonly used algorithms are: the Expectation Maximization algorithm (EM) and the Conjugate Gradient algorithm (CG).

1. EM algorithm

This algorithm was developed and proposed by Lange and Carson. The formula of this algorithm is the following [6]:

$$f^{n+1}=f^n{R}^t\frac{p}{R.f^n}$$E24

Where n is the number of the actual iteration. This algorithm is characterized by the fact that it keeps the number of iterations for each projection. Moreover its multiplicative form gives a positivity constraint although it implies a slow convergence.

1. Conjugate gradient algorithm (CG)

It is an algorithm which is used because it converges very quickly. It is based on a classical descent method. Its formula of iterative updating can be, roughly, given by the following expression [7]:

$$f^{n+1}=f^n+{\alpha }^n{d}^n$$E25

We can verify that the correction is not multiplicative as for EM algorithm, but it is additive. This formula is characterized by a descent direction d (Eq. 26) and a descent speed α (Eq.27) that are recalculated in a conjugated manner for each iteration and this to optimize the speed of convergence [6].

$${\alpha }^n=\frac{{\Vert R^{t}p-R^{t}R{f}^n\Vert }^2}{{(R^{t}p-R^{t}R{f}^n)}^{t}R(\text{Rt p - Rt R fn)}}$$E26

$$\begin{array}{cc}iteration.1:& d^1=R^{t}p-R^{t}R{f}^0\\ iteration.2:& d^n=(R^{t}p-R^{t}R{f}^n)+b^n{d}^{n-1}\end{array}$$E27

Through this process, the error between the measured projections and those calculated is minimized progressively. This error evaluation formula is given by:

$$e={\Vert p-Rf\Vert }^2$$E28

1. ART algorithm

The formula of the iterative updating of this algorithm is given by the following expression:

$$f_i^{(n+1)}=f_i^{(n)}+R_{ji}\frac{p_j-R_j{f}^{(n)}}{{\Vert R_j\Vert }^2}$$E29

The correction in this algorithm can be additive (e=p_k-p_kⁿ) or multiplicative (e= p_k/p_kⁿ).

1. SIRT algorithm

The formula of the iterative updating of this algorithm is given by the following expression:

$$f_i^{(n+1)}=f_i^{(n)}+\frac{{\displaystyle {\sum }_j{p}_j}}{{\displaystyle {\sum }_j{\displaystyle {\sum }_i{R}_{ji}}}}-\frac{{\displaystyle {\sum }_j{R}_j{f}^{(n)}}}{{\displaystyle {\sum }_j{\Vert R_J\Vert }^2}}$$E30

2.3. Iteration stopping rules

The iterative process for image reconstruction must be stopped at the end of the convergence phase before it starts diverging. This can be done using some regularization methods. Indeed, well selected regularization parameter can controls the amount of stabilization imposed on the solution. In iterative methods one can use the stopping index as regularization parameter. When an iterative method is employed for image reconstruction, the user can also study on-line adequate visualizations of the iterates as soon as they are computed, and simply halt the iteration when the approximations reach the desired quality. This may actually be the most appropriate stopping rule in many practical applications, but it requires a good intuitive imagination of what to expect. If this is not the case, a computer can give us some aid to determine the optimal approximation. The stopping rules are divided into two categories: rules which are based on knowledge of the norm of the errors, and rules which do not require such information. If the error norm is known within reasonable accuracy, the perhaps most well known stopping rule is the discrepancy principle Morozov [8]. Examples of the second category of methods are the L-curve criterion [9], and the generalized cross-validation criterion [10].

3. Conclusions

After a substantial effort, major breakthroughs have been achieved in the last fourteen years in the mathematical modeling of CT. The aim of this chapter is to survey this progress and to describe the relevant models, mathematical problems and reconstruction procedures used in CT. We give a summary on the mathematics of computed tomography. We start with a short introduction to integral computed tomography. We then go over to projection and inversion algorithms. We give a detailed analysis of the filtered back-projection algorithm in the light of the sampling theorem. We also describe the convergence properties of iterative algorithms. We shortly mention Fourier based algorithms used in CT.

Good reconstructions from the interlaced lattice can also be obtained by using the direct algebraic reconstruction algorithm, or by increasing the amount of data through the interpolation according to the sampling theorem. As we have seen, the interpolation step can introduce significant errors in certain cases. It has also been shown that the interpolation can be avoided by choosing the points x where the reconstruction is computed on a polar grid rather than on a square Cartesian grid, and interchanging the order of the two summations. This algorithm should work well for the interlaced lattice and is particularly beneficial in case of the fan-beam sampling geometry, since the method also avoids the homogeneous approximation, whose influence on the reconstruction is difficult to estimate.

The cone and fan beam scanning are the standard scanning modes in present day’s clinical practice. The methods described above assumed that the geometry of the acquisition was a parallel geometry, like that of first generation systems. In the case a conical geometry (or fan), three methods are possible:

1. The first is to ignore the discrepancy. The error induced by this approximation is considered negligible if the beam angle is low (typically below 15 degrees). This method is applicable to systems of second generation, but more to the following where the beam must cover the whole section of the patient to remove the translation.

2. The second method is to rearrange the data in parallel projections, mainly by interpolation.

3. The last method is to reformulate the problem completely. It becomes apparent that the projection theorem cannot be generalized, which does not have direct inversion method. The projection theorem can be adapted to different geometries, and then the algorithm is the same as for a parallel geometry, with the same calculations. Finally, the filtered back projection formulas can be corrected and result in a slightly different algorithms. The algebraic methods have some fewer additional problems, but the flexibility in the choice of basic functions and thus the coefficients of the matrix R, allows these methods to be adapted to different geometries.

Finally, I want just to add that the purpose of this chapter is to give an introduction to the studied topic and treat some related aspects in more detail. The reader interested in a broader overview of the field, its relation to various branches of pure and applied mathematics, and its development over the years may wish to consult the appropriate bibliography.