## Abstract

This paper presents a brief overview of some basic iterative algorithms, and more sophisticated methods are presented in the research papers in this issue. A range of algebraic iterative algorithms are covered here including ART, SART and OS-SART. A major limitation of the traditional iterative methods is their computational time. The Krylov subspace based methods such as the conjugate gradients (CG) algorithm and its variants can be used to solve linear systems of equations arising from large-scale CT with possible implementation using modern high-performance computing tools. The overall aim of this theme issue is to stimulate international efforts to develop the next generation of X-ray computed tomography (CT) image reconstruction software.

## 1. Basic principle of computed tomography

X-ray computed tomography (CT) plays a key role in many scientific discoveries [1]. X-ray CT image recons- truction is an exciting mathematical and computational challenge. There is growing interest in iterative reconstruction algorithms mainly due to recent surge in computational power. In this paper, a brief overview of iterative algorithm in X-ray CT is presented [2–7].

When the monochromatic rays penetrate through a homogeneous object as shown in figure 1*a*, the ray attenuation can be explained by the Lambert–Beer Law as follows: *I*=*I*_{0}⋅e^{−μx}, where *I*_{0} and *I* are an incident photon intensity and an intensity of photons after having passed an object with thickness *x*. The purpose of CT is to determine the linear attenuation coefficient *μ* of the object and convert those values to a grey-scale image. The attenuation coefficient is measured in units per metre. Unfortunately, the body organs which photons pass through are not homogeneous tissues as shown in figure 1*b*. The Lambert–Beer Law must be modified to *I*=*I*_{0}e^{−(μ1+μ2+μ3+⋯+μn)x}. Moreover, the human organs are not only inhomogeneous tissues but also the tissue thicknesses are different, as illustrated in figure 1*c*. In this model, the numbers of photons are calculated by
1.1
1.2
1.3

For clinical CT scanners, an X-ray source produces an X-ray beam with a wide range of photon energies (polychromatic radiation). Therefore, equation (1.3) can be extended to
1.4
and
1.5
The summation of the attenuations in the particular X-ray path can be defined as
1.6
The natural logarithm of the ratio of the incident photons and the outgoing photons is known as ray sums. A set of ray sums at a particular angle is a projection as shown in figure 2*a*.

A back projection algorithm is typically used for calculating the linear attenuation coefficients *μ* in an image matrix of an object by projecting all projections backwards in the direction of each ray path. In the case of small numbers of projections, the reconstructed image would show a spiky artefact (figure 2*b*). This is also known as a star artefact. On the other hand, blurring around the edges of the images would be developed from spikiness with more projections used. In order to remove this artefact, a mathematical filtering process can be used. Digital filters are implemented to raw projection data and the results, called filtered projection, are back-projected to the image matrix. Next a brief overview of a number of iterative algorithms will be presented.

## 2. Algebraic reconstruction methods: ART, SART and OS-SART

An iterative reconstruction algorithm is a technique that uses the differences between the measured data and the calculated data to update an image. In this section, the approaches and definitions for SART and OS-SART in Jiang and Wang are described. For example, an image **x** is required which solves
2.1
where ** x**∈ℜ

^{N}is the image value for

*N*cells and

**b**∈ℜ

^{M}is measured data from

*M*measurements. Each element

*A*

_{mn}of the weight matrix makes

**∈ℜ**

*A*^{M×N}is the length of intersection of the

*m*th ray with the

*n*th cell, for example is the weight matrix of the intersection between the 1st ray and the image matrix 2-by-2 as illustrated in figure 3. For iterative reconstruction, equation (2.3) can be solved by 2.2

where *δ***x**^{(k)} is an update, λ_{n} is a relaxation parameter, *k* is the iteration index and *n* is the cycle index (one cycle is defined as all available datasets used in one iteration); let
2.3
Usually, the initial image values are defined as all zero, then **x**^{(0)}=**0**. Also, let
2.4
and
2.5
ART updates row-by-row so
2.6
where *i* is the row index, and the relaxation parameter λ_{n} is changed by the cycle index *n*. SART iteratively updates for the whole image as
2.7
where **V** is diagonal with *j*th diagonal element
2.8
and **W** is diagonal and with *i*th diagonal element
2.9
and **A**_{ij} is the (*i*,*j*)th element of **A**. However, **b** can also be divided into *B* blocks where
2.10
and . Similarly **A** can be divided into *B* blocks where
2.11
The reconstruction method of divided blocks is the so-called ‘OS-SART’ which updates an image block-by-block using
2.12
where **W**_{β} is diagonal with *i*th diagonal element
2.13
To counteract the potential divergence of reconstruction after some cycles, in this study a relaxation parameter λ_{n} is scheduled corresponding to the cycle index *n*(*n*≥1) as
2.14
and sf is the scaling factor. Denote that *α* is a rate of change of the relaxation parameter. This guarantees the convergence of the algorithm.

## 3. Conjugate gradient least-squares algorithm

The conjugate gradient least squares (CGLS) is an iterative solution used in large sparse least-squares problem of **Ax**=**b**. The weight matrix (** A**∈ℜ

^{M×N}) is created using a forward projection programme. The matrix

**A**is a combination of the sub-matrix for each projection. In the case where the matrix

**A**is positive-definite, the matrix

**A**

^{T}

**A**is positive-definite for any matrix

**A**. The iterative method is terminated at

*i*steps for

*k*=0,1,2,…,

*i*. The residual (

**r**

^{k}) at each step is computed by 3.1 Starting with an initial approximation (

**x**

^{0}), then,

**r**

^{0}=

**b**−

**Ax**

^{0},

**p**

^{0}=

**A**

^{T}

**r**

^{0}3.2 3.3 3.4 3.5 3.6 Normally, an initial approximation is

**x**

^{0}=0, then

**r**

^{0}=

**b**and

**p**

^{0}=

**A**

^{T}

**b**, where

**p**is conjugate vectors. In the case of the residual

**r**being zero, it implies that the problem is solved. Otherwise, if the residual

**r**is non-zero, the desired solution for

**r**

^{k}to be minimized, which can be monitored by the L2 norm of the

**r**

^{k}for each iteration. This can be achieved by updating the residual

**r**into the problem iteratively. Moreover, regularized CGLS was also studied here for ill-posed problems.

## 4. Conclusion

A brief overview of a number of iterative algorithm for X-ray tomography is presented in this article. The hope is that the next generation of X-ray tomography algorithms to be further developed in the next few years are better adapted in real application of X-ray tomography devices.

## Footnotes

One contribution of 11 to a theme issue ‘X-ray tomographic reconstruction for materials science’.

- Accepted March 11, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.