A priori information in a regularized sinogram-based method for removing ring artefacts in tomography

Titarenko, S.; Titarenko, V.; Kyrieleis, A.; Withers, P.J.

doi:10.1107/S0909049510010964

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 17| Part 4| July 2010| Pages 540-549

doi:10.1107/S0909049510010964

A priori information in a regularized sinogram-based method for removing ring artefacts in tomography

Sofya Titarenko,^a Valeriy Titarenko,^b ^* Albrecht Kyrieleis ^b and Philip J. Withers ^b

^aDepartment of Mathematics, Faculty of Physics, Moscow State University, Leninskie Gory, Moscow 119991, Russia, and ^bHenry Moseley X-ray Imaging Facility, School of Materials, University of Manchester, Grosvenor Street, Manchester M1 7HS, UK
^*Correspondence e-mail: valeriy.titarenko@manchester.ac.uk

(Received 18 September 2009; accepted 23 March 2010; online 14 May 2010)

Ring artefacts in X-ray computerized tomography reconstructions are considered. The authors propose a ring artefact removal method based on a priori information regarding the sinogram including smoothness along the horizontal coordinate, symmetry of the first and the final rows and consideration of small perturbations during acquisition. The method does not require prior reconstruction of the original or corrected sinograms. Its numerical implementation is based on quadratic programming. Its efficacy is examined with regard to experimental data sets collected on graphite and bone.

Keywords: inverse problems; image processing; flat-field correction; tomography; quadratic programming; artefact.

1. Introduction

We consider a typical arrangement for parallel-beam tomography, whereby a sample is rotated about an axis while illuminated by a parallel beam of X-rays (see Fig. 1). This does not restrict the method proposed in the paper, which can be also applied (with small modifications) to fan-beam and cone-beam geometries. A description of these geometries may be found by Natterer & Wübbeling (2007 ).

Figure 1
A simplified scheme of the experimental set-up.

Typically, a white or monochromatic beam of X-rays formed by a synchrotron or laboratory source passes through the sample and then falls onto a scintillator, which absorbs X-rays and then re-emits the absorbed energy in the form of visible light. A charge-coupled device (CCD) or a complementary metal-oxide-semiconductor (CMOS) camera then records the photons to form a radiographic image.

Unfortunately, the real system is always affected by errors; some of them are random (e.g. cosmic rays, dark noise on the pixels) while others have a deterministic nature. The following causes lead to ring artefacts in reconstructed slices (see, for example, Banhart, 2008 ; Vidal et al., 2005 ; Kowalski, 1977 ; Titarenko et al., 2010b ):

(i) dirt, dust or scratches on the scintillator, the optical system or the camera;

(ii) X-ray diffraction, reflection of the visible light from some part of the system;

(iii) defects in pixel elements such as dead pixels, drifts in the individual pixel response, or non-linear response;

(iv) defects or impurities on the scintillator crystals;

(v) beam variations either produced at the source or somewhere along its path (e.g. owing to defects in the monochromator).

Ring artefacts often appear as concentric circles (for 360° sample rotation) or arcs (half-circles with additional line segments on their ends, when the sample is rotated by 180°). The common means of reducing ring artefacts is the flat-field correction method [for an introduction, see Stock (2008 )]. Consider a row of pixels along x on the camera recording the projection of a two-dimensional horizontal slice (horizontal section across the sample) recorded with the sample rotated to an angle θ about a vertical axis. Then having measured flat-field intensity $[I_{\rm{flat}}(x,\theta)]$ (the intensity recorded without the sample), a dark-field intensity $[I_{\rm{dark}}(x,\theta)]$ (the intensity recorded without the beam, caused by the generation of thermal electrons on the camera pixels) and the sample illuminated intensity $[I(x,\theta)]$ , we calculate the attenuation through the slice,

$[p(x,\theta)=\ln\left[{{I_{\rm{flat}}(x,\theta) -I_{\rm{dark}}(x,\theta)}\over{I(x,\theta)-I_{\rm{dark}}(x,\theta)}}\right].\eqno(1)]$

The flat-field and the dark-field intensities are often taken only before and/or after the experiment and are interpolated for all other angles θ. To increase the accuracy several projections of the flat-field and the dark-field are sometimes taken and their averaged projections are used. This method is not entirely successful in suppressing the error (see, for example, Davidson et al., 2003 ), and so a more effective way to improve the image is needed. A number of approaches have already been proposed, and they can be divided into three main groups:

(i) Processing before reconstruction. Such methods often use the intensity $[I(x,\theta)]$ or the attenuation $[p(x,\theta)]$ sinogram and are based on smoothing the sinogram [either ordinary smoothing, filtering using a fast Fourier transform or the removal of high frequencies; see, for example, Boin & Haibel (2006 ), Raven (1998 ), Antoine et al. (2002 )].

(ii) Processing after reconstruction (for example, see Walls et al., 2005 ; Sijbers & Postnov, 2004 ; Yang et al., 2008 ). For instance, in Yang et al. (2008) the reconstructed slice is considered in polar coordinates, if a large variance occurs, an average value taken over a few pixels in the radial direction.

(iii) Modifying the experimental procedure. For example, using time delay integration whereby the detector is moved laterally during acquisition (Davis & Elliott, 1997 ).

Of course, ring artefact suppression methods can also use ideas from several groups. For example, the method proposed by Titarenko et al. (2009 ) is based on both pre- and post-processing ideas.

We propose an alternative method that can be applied before reconstruction. Consider the ideal case where the X-ray beam is monochromatic and constant over time, and the dark noise is zero; then the measured intensity is

$[I(x)=I_0(x)\exp\left[-p(x,\theta)\right],\eqno(2)]$

where I₀(x) is the initial (flat-field) intensity and $[p(x,\theta)]$ is the attenuation. If we assume that there is an unknown non-rotating object between the sample and the scintillator, then we could write

$[I(x)=I_0(x)\exp\left[-p(x,\theta)-q(x)\right],\eqno(3)]$

where q(x) is a function to be found and describes the overall contribution to the recorded intensity from all the deterministic errors listed above. In order to describe a range of physical effects we do not restrict q(x) to have only positive values. The effects of non-linearities on the scintillator or the detector pixels are not considered here, since they are not described by (3). This is also strictly true for white-beam optics; however, the method can also be used at least as a first approximation (see examples in §4).

There may be various random errors $[\varepsilon(x,\theta)]$ in the detected signal; for example, as caused by beam instability, dark noise in the pixels, analog-to-digital conversion errors associated with the camera. Hence, the measured attenuation $[\tilde p(x,\theta)]$ can be written as

$[\tilde p(x,\theta)=p(x,\theta)+q(x)+\varepsilon(x,\theta).\eqno(4)]$

Our aim is to find q(x); error reduction is not considered in this paper.

Note that it is not possible to solve (4) without employing some assumptions about the nature of $[p(x,\theta)]$ and q(x). In order to solve the equation we mainly invoke principles of the theory of inverse and ill-posed problems described below.

2. Ill-posed problems

Let there be a problem

$[Az=u,\qquad u\in U,\quad z\in Z,\eqno(5)]$

where U and Z are Hilbert spaces and $[A:Z\to U]$ is a linear operator. Instead of the exact operator A and the data u we are given approximate A_h and $[u_\delta]$ such that

$[\|A_h-A\|\leq h,\qquad \|u_\delta-u\|\leq \delta.\eqno(6)]$

Denote $[\eta]$ = $[\{h,\delta\}]$ . Our aim is then to find a solution $[z\in Z]$ for the given data $[\{A_h,u_\delta,\eta\}]$ .

In the earliest part of the 20th century Jacques Hadamard formulated the conditions which define the so-called well posed problems (see Hadamard, 1923 ), namely:

(i) the solution of the problem expressed in (5) exists and;

(ii) it is unique;

(iii) it is stable.

All problems which do not satisfy the above conditions are called ill-posed. Hadamard believed that solving practical problems one had to consider only well posed ones. However, since then many important ill-posed problems have been identified.

In 1963, A. N. Tikhonov formulated a famous definition of the regularizing algorithm that has become a basic tenet in the modern theory of ill-posed problems (see Tikhonov, 1963 ). By definition, a regularizing algorithm (regularizing operator) is called an operator/function $[R(A_h,u_\delta,\eta)]$ possessing two properties:

(i) $[R(A_h,u_\delta,\eta)]$ is defined for any $[\delta]$ > 0, $[h\geq 0]$ , $[u_\delta\in U]$ , $[A_h \in L(Z,U)]$ ;

(ii) for any $[z\in Z]$ , for any $[u_\delta \in U]$ such that Az = u, $[\|u-u_\delta \|\leq \delta]$ , $[\delta]$ > 0 and for any $[A_h \in L(Z,U)]$ such that $[\|A_h -A\|\leq h]$ , $[h\geq 0]$ , $[z_\eta]$ = $[R(A_h,u_\delta,\eta) \to z]$ as $[\eta \to 0]$ .

Here L(Z,U) is a space of linear bounded operators acting from Z into U.

Tikhonov proposed the following approach to construct regularizing algorithms. He introduced the smoothing functional (now often called Tikhonov's functional). The simplest form of this functional is

$[M^{\alpha}[z]=\|A_hz-u_\delta\|^2+\alpha\|z\|^2,\eqno(7)]$

where $[\alpha]$ > 0 is a regularization parameter. The minimiser $[z^\alpha_\eta]$ of the functional is considered as an approximate solution of (5). Note that if there are several solutions of (5) then $[z^\alpha_\eta]$ tends to the normal solution (i.e. the solution with the minimal norm in Z). And if (5) has no solution, then $[z^\alpha_\eta]$ tends to a normal pseudosolution (pseudosolutions are defined as minimisers of $[\|Az- u\|]$ ).

The main problem is the choice of the regularization parameter $[\alpha]$ . There are two approaches: a priori and a posteriori. The a priori choice is demonstrated by the following theorem (Tikhonov & Arsenin, 1977 ).

Theorem 1. Let A be a one-to-one operator, $[z \in Z]$ . Then $[z^{\alpha(\eta)}_\eta\to z]$ as $[\eta\to 0]$ and $[\alpha(\eta)\to 0]$ so that $[(h+\delta)^2/\alpha(\eta)\to 0]$ .

As one of the a posteriori choices we consider the simplest form of the generalized discrepancy principle (Tikhonov et al., 1995 ). Let a solution $[z^\alpha_\eta]$ be found for a given $[\alpha]$ > 0 and define a generalized discrepancy as

$[\rho_\eta(\alpha)=\|A_hz^\alpha_\eta-u_\delta\|^2-\left(\delta+h\|z^\alpha_\eta\|\right)^2\,\,-\,\,\sigma^2,\eqno(8)]$

where $[\sigma\equiv \min_{z\in Z} \|A_h z- u_\delta\|]$ .

(1) If the condition $[\|u_{\delta}\|^2\,\,\gt\,\,\delta]$ does not hold, then as an approximate solution of (5) we take $[z_{\eta}]$ = 0.

(2) Otherwise, the generalized discrepancy $[\rho_\eta(\alpha)]$ has a root $[\alpha^*\,\,\gt\,\,0]$ and we take $[z_{\eta}]$ = $[z_{\eta}^{\alpha^*}]$ , and if $[\rho_{\eta}(\alpha)\,\,\gt\,\,0]$ for all $[\alpha\,\,\gt\,\,0]$ then $[z_{\eta}]$ = $[\lim_{\alpha\to 0+0}z^{\alpha}_{\eta}]$ .

To find the root one can exploit the properties of $[\rho_\eta(\alpha)]$ :

(1) $[\rho_\eta (\alpha)]$ is continuous and monotonically non-decreasing for $[\alpha\,\,\gt\,\,0]$ .

(2) $[\lim_{\alpha \to +\infty} \rho_\eta (\alpha)]$ = $[\|u_\delta\|^2 -\delta^2-\sigma^2]$ .

(3) $[\lim_{\alpha \to 0+0} \rho_\eta (\alpha) \leq -\delta^2]$ .

It is important to mention that for ill-posed problems the regularization parameter $[\alpha]$ should depend explicitly on the errors. If there is no such dependency then only well posed problems can be solved (Bakushinskii, 1984 ). Therefore the best known `error-free methods', namely the `L-curve method' (Hansen, 1992 ) and `the generalized cross-validation method' (Wahba, 1977 ), cannot be applied to ill-posed problems [see examples by Yagola et al. (2002 )].

It has been recognized that the problem (5) could be generalized for non-linear problems (Tikhonov et al., 1998 ). In such a case one wants to minimize a functional J(z). If the solution is not unique, one chooses the solution which minimizes a functional $[\Omega(z)]$ . Then for the simplest case we obtain the modified smoothing functional

$[M^\alpha[z]=J_\eta(z)+\alpha\Omega(z),\eqno(9)]$

where $[J_\eta(z)]$ is an approximate functional [see Tikhonov et al. (1998) for more details]. For (5), $[\Omega(z)]$ = $[\|z\|^2]$ , J(z) = $[\|Az-u\|^2]$ , $[J_\eta(z)]$ = $[\|A_hz-u_\delta\|^2]$ .

More information about ill-posed problems can be found in books by Tikhonov et al. (1995), Bakushinsky & Kokurin (2004 ), Engl et al. (1996 ), Groetsch (1993 ), Ivanov et al. (2002 ), Kaltenbacher et al. (2008 ), Lavrentiev et al. (2003 ) and Tikhonov et al. (1998).

3. Ring artefact removal algorithm

3.1. A priori information

3.1.1. Smoothness

Here we discuss the solution of (4). The first step is to find additional (so-called a priori) information about the exact solution $[p(x,\theta)]$ . The main assumption to be made is that $[p(x,\theta)]$ is a smooth function of x for the following reason.

Let us denote the intensity of the X-ray beam just before the scintillator as I(x,z), where z is the vertical coordinate. The intensities of the visible light just after the scintillator and after the optical system are denoted as I_scin(x,z) and I_opt(x,z), respectively. It may be assumed that transformations $[I(x,z)\to]$ $[I_{\rm{scin}}(x,z)\to]$ I_opt(x,z) are linear and could be written as convolutions, i.e.

$[\eqalignno{I_{\rm{scin}}(x,z)&=I(x,z)*K_{\rm{scin}}(x,z) \cr&\equiv \textstyle\int\!\!\int I(x-\xi,z-\zeta)K_{\rm{scin}}(\xi,\zeta)\,{\rm{d}}\xi\,{\rm{d}}\zeta&(10)}]$

and I_opt(x,z) = I_scin(x,z) * K_opt(x,z) = I(x,z) * K(x,z), K(x,z) = K_opt(x,z)*K_scin(x,z), where K_scin(x,z) and K_opt(x,z) are point-spread functions (PSFs) of the scintillator and the optical system. It is often assumed that the PSFs have Gaussian shape (Banhart, 2008; Mahajan, 2001 $[Mahajan, V. N. (2001). Optical Imaging and Aberrations, Part II, Wave Diffraction Optics. SPIE Press.]$ ). Thus we can write K(x,z) = $[k\exp[-\beta(x^2+z^2)]]$ , where $[k, \beta\,\,\gt\,\,0]$ . Owing to the differentiation characteristics of the convolution we obtain

$[{{\partial}\over{\partial x}} I_{\rm{opt}}(x,z) = I(x,z)*{{\partial}\over{\partial x}} K(x,z).\eqno(11)]$

Since $[{{\partial}}K(x,z)/{\partial x}]$ = $[-2k\beta x \exp[-\beta(x^2+z^2)]]$ and $[I(x,z)\leq]$ I_max, then

$[\eqalignno{\left|{{\partial}\over{\partial x}}I_{\rm{opt}}(x,z)\right|&\leq 2k\beta I_{\rm{max}} \int\!\!\int \left|\xi \exp\left[-\beta(\xi^2+\zeta^2)\right]\right|\,{\rm{d}}\xi\,{\rm{d}}\zeta\cr& \leq 2k\beta I_{\rm{max}} \, \left(\pi^{1/2}/\beta\beta^{1/2}\right) \cr& = 2k I_{\rm{max}} (\pi/\beta)^{1/2}.&(12)}]$

Thus the attenuation $[p(x,\theta)]$ is also a smooth function of x. For example, let us consider an object with a flat edge, so the intensity of X-rays after they have passed through the sample is

$[I(x,z) = \Bigg\{\matrix{I_{\rm{max}}\hfill & x \,\,\lt\,\,\, 0,\hfill\cr \sigma I_{\rm{max}}\hfill & x\,\geq\, 0,\hfill}\eqno(13)]$

where $[\sigma\in[0,1]]$ , i.e. the sample is in the left half-plane ( $[x\geq 0]$ ) and the optical path length of the X-rays passing through the sample is $[-\ln\sigma]$ . Then the intensity of the visible light recorded by the camera sensor is

$[\eqalignno{I_{\rm{opt}}(x,z)&= K(x,z)\,I(x,z) \cr&= I_{\rm{max}}k\textstyle\int\limits_{-\infty}^\infty \exp\left(-\beta\eta^2\right)\,{\rm{d}}\eta\Big[\textstyle\int\limits_x^\infty1\exp\left(-\beta\xi^2\right)\,{\rm{d}}\xi \cr&\quad + \textstyle\int\limits_{-\infty}^x\sigma\exp\left(-\beta\xi^2\right)\,{\rm{d}}\xi\Big] \cr& = I_{\rm{max}}k(\pi/\beta)\left(\pi^{1/2}/2\beta^{1/2}\right) \Big\{{\rm{erf}}(+\infty)-{\rm{erf}}\left(\beta^{1/2}x\right)\cr&\quad+\sigma\Big[{\rm{erf}}\left(\beta^{1/2}x\right)-{\rm{erf}}(-\infty)\Big]\Big\} \cr& = I_{\rm{max}}k(\pi/2\beta) \left[1+\sigma-{\rm{erf}}\left(\beta^{1/2}x\right)(1-\sigma)\right],&(14)}]$

where $[{\rm{erf}}(x)\equiv ({{2}/{{\pi^{1/2}}}})\!\textstyle\int_0^x\!\exp({-t^{\,2}})\,{\rm{d}}t]$ . If $[I_1\equiv]$ $[\lim_{x\to-\infty} I_{\rm{opt}}(x,z)]$ = $[I_{\rm{max}}k\pi/\beta]$ , then I_opt(x,z) = $[I_1 [1+\sigma-{\rm{erf}}({\beta^{1/2}}x)(1-\sigma)]/2]$ ,

$[{{\partial}I_{\rm{opt}}(x,z)/{\partial x}}=-I_1\left({{{\beta^{1/2}}}/{{\pi^{1/2}}}}\right)(1-\sigma) \exp(-\beta x^2)\eqno(15)]$

and $[{{\partial I_{\rm{opt}}(x,z)}/{\partial z}}]$ = 0, therefore the maximum absolute value of the first derivative is $[I_1{\beta^{1/2}}(1-\sigma)/{\pi^{1/2}}]$ . So instead of the discontinuous function I(x,z) we record the continuous function I_opt(x,z). The values of I_opt(x,z) at two points (x₁,z₁) and (x₂,z₂) can differ by at most

$[I_1{\beta^{1/2}}(1-\sigma)\left[\left(x_1-x_2\right)^2+\left(z_1-z_2\right)^2\right]^{1/2}/{\pi^{1/2}}.\eqno(16)]$

However, if the PSFs tend to two-dimensional delta functions (e.g. when the scintillator is very thin), then $[\beta]$ increases and a possible difference of the intensities between two points may be very large. Therefore intensities of light recorded at neighbouring pixels of the camera sensor may be uncorrelated and the approach described in the paper cannot be used. However, for most data sets collected at modern synchrotrons the assumption of smoothness is valid. Of course, one may expect that scintillators will become thinner in the near future; however, it is likely that the sizes of camera pixels will also decrease, so maximal changes of intensities at neighbouring pixels should be similar.

We also suppose that the optical path lengths found using images recorded by the camera and real X-ray images virtually recorded before the scintillator are the same. Therefore the real optical path length is assumed to be a smooth function on (x,z).

3.1.2. Small perturbations

Let us introduce a Cartesian frame $[(\xi,\eta)]$ on a given slice. Suppose (0,0) is at the point of intersection of the vertical axis of rotation of the sample and the given slice, and the coordinate axis $[\xi]$ is parallel to the axis x on the plane of the camera sensor. For simplicity and not restricting generality we suppose that (0,0) projects into on the axis x. Denote the attenuation coefficient across the given slice by $[\mu(\xi,\eta)]$ . Then the attenuation $[p(x,\theta)]$ is the line integral

$[p(x,\theta)=\textstyle\int\limits_L\mu(\xi,\eta)\,{\rm{d}}\xi\,{\rm{d}}\eta,\eqno(17)]$

where L is defined by the equation $[\xi\cos\theta]$ + $[\eta\sin\theta]$ = x.

Let there be a radially symmetric function $[\tilde\mu(\xi,\eta)]$ = $[f(\xi^2+\eta^2)]$ . The corresponding attenuation $[\tilde p(x,\theta)]$ defined by (17) does not depend on the angle θ, i.e. $[\tilde p(x,\theta)]$ = $[\tilde p(x)]$ . If there is no information about the exact $[p(x,\theta)]$ and q(x) in (4), then there may be several solutions. For example, we may define $[\mu_\tau(\xi,\eta)]$ = $[\tau\tilde\mu(\xi,\eta)]$ and $[q_\tau(x)]$ = $[(1-\tau)\tilde p(x)]$ , then we get the corresponding $[p_\tau(x)]$ = $[\tau\tilde p(x)]$ and $[p_\tau(x)+q_\tau(x)]$ = $[\tilde p(x)]$ . Even if we restrict $[\mu(\xi,\eta)]$ to non-negative functions only, then $[\tau\geq 0]$ and we get infinitely many pairs $[(p_\tau,q_\tau)]$ that give us the same $[\tilde p(x)]$ . This example shows that the solution of (4) may be non-unique. We need an additional assumption which allows us to find a unique solution among all possible solutions satisfying (4). If there is no a priori information about $[\mu(\xi,\eta)]$ , then the most natural requirement is that q(x) is small. In the case of large dust/dirt particles with high attenuation we cannot guarantee that there is no other `virtual' solution $[\mu_v(\xi,\eta)]$ that gives $[p_v(x,\theta)]$ and this $[p_v(x,\theta)]$ is closer to the recorded $[\tilde p(x,\theta)]$ than the real one $[p(x,\theta)]$ . Therefore if there are no other reasons we choose a solution $[\mu(\xi,\eta)]$ such that the corresponding $[p(x,\theta)]$ is close to $[\tilde p(x,\theta)]$ . So if several solutions $[p(x,\theta)]$ exist, then we choose the solution that has the minimal norm of $[\|q(x)\|]$ . Of course, the term `small' may be defined in different ways. For instance, if one needs only small absolute values of q(x), then one can use the L² norm, where $[\|q(x)\|^2]$ = $[\textstyle\int q^2(x)\,{\rm{d}}x]$ ; if the first derivative of q(x) should also be `small', then one can use the W²₁ norm $[\|q(x)\|^2]$ = $[\textstyle\int[q^2(x)+\{q'(x)\}^2]\,{\rm{d}}x]$ .

3.1.3. Evenness

As we have supposed that the axis of rotation intersects the given slice at (0,0) and this point projects into x = x_c, then

$[p\left[-\left(x-x_c\right),\theta+180^\circ\right]=p\left(x-x_c,\theta\right),\eqno(18)]$

(see, for example, Natterer & Wübbeling, 2007). Not restricting generality, we set x_c = 0 and

$[p\left(-x,\theta+180^\circ\right)=p(x,\theta).\eqno(19)]$

Physically, it simply means that an X-ray beam attenuates by the same amount whichever direction it travels along a given path through the sample, so that the same attenuation is recorded for the beam path when it is rotated by 180°. Thus if the sample is rotated by $[180^\circ]$ , then there will be the following requirement:

$[\textstyle\int\left\{\,p(x,0)-p\left(-x,180^\circ\right)\right\}^2\,{\rm{d}}x=0.\eqno(20)]$

If the sample is rotated by $[\Theta\in(180^\circ\semi360^\circ],]$ then

$[\textstyle\int\limits_0^{\Theta-180^\circ}\,{\rm{d}}\theta\textstyle\int{\rm{d}}x\left\{\,p(x,\theta)-p\left(-x,\theta+180^\circ\right)\right\}^2=0.\eqno(21)]$

3.1.4. Other possible requirements

Sometimes it is known that for some regions the attenuation has a known value, e.g. zero. For instance, a sample may have several holes such that for some angles θ certain pixels of the camera measure only a flat-field. Then one can find q(x) for some intervals of x.

As a consequence if it is known that for some projection angle θ part of the illuminating beam does not intersect with the sample, then this can be used to help find the exact solution. Let the sample be rotated by 180° and $[p(x,0)\equiv 0]$ for $[x\leq 0]$ , then $[p(x,180^\circ)\equiv 0]$ for $[x\geq 0]$ and

$[q(x)=\Bigg\{\matrix{\tilde p(x,0),\hfill & x\,\geq\,0,\hfill\cr \tilde p(x,180^\circ),\hfill & x\,\,\lt\,\,\,0.\hfill}\eqno(22)]$

Let the sample have a known radius R. By the radius we mean the maximum distance between the origin (0,0) and any point that belongs to both the sample and the given horizontal slice. This does not restrict the sample to have a specific shape, any real sample is permitted; however, the result below will explicitly depend on R. For $[[x_1,x_2]\subset[-R,R]]$ we define a strip $[D \equiv \{(\xi,\eta): \xi^2+\eta^2\leq R, \xi\in[x_1,x_2]\}]$ (see Fig. 2). The integral of $[\mu(\xi,\eta)]$ on D is invariant with respect to a rotation of the $[(\xi,\eta)]$ coordinate system by an angle $[\triangle \theta]$ . This means that if we virtually exclude all parts of the sample that are outside D, then the integral path length is the same whether the sample is rotated by any angle $[\triangle \theta]$ . Since, for a given θ, $[p(x,\theta)]$ = $[\textstyle\int_{-\infty}^{+\infty}\mu(\xi,\eta)|_{\xi = x}\,{\rm{d}}\eta]$ , i.e. the integral over the line $[\xi = x]$ , then the integral path length is

$[\eqalignno{\textstyle\int\limits_{x_1}^{x_2}p(x,\theta)\,{\rm{d}}x&= \textstyle\int\limits_{x_1}^{x_2}\textstyle\int\limits_{-\infty}^{+\infty}\mu(\xi,\eta)|_{\xi=x}\,{\rm{d}}\eta\,{\rm{d}}x\cr& = \textstyle\int\limits_{x_1}^{x_2}\textstyle\int\limits_{-\infty}^{+\infty}\mu(\xi,\eta)\,{\rm{d}}\eta\,{\rm{d}}\xi\cr&= \textstyle\int\!\!\textstyle\int\limits_D\mu(\xi,\eta)\,{\rm{d}}\eta\,{\rm{d}}\xi,&(23)}]$

i.e. the area integral, and the area integral does not depend on the angle of rotation of the coordinate system, i.e. on $[\triangle \theta]$ . In this case D becomes $[D_{\triangle \theta}]$ in new coordinates. Define x₁^* = $[\min_{D_{\triangle \theta}} \xi]$ and x₂^* = $[\max_{D_{\triangle \theta}} \xi]$ . Since the area $[D_{\triangle \theta}]$ is inside the area $[D^*\equiv \{(\xi,\eta): \xi^2+\eta^2\leq R^2, \xi\in[x_1^*,x_2^*]\}]$ , then $[\textstyle\int_{D^*}\mu(\xi,\eta)\,{\rm{d}}\xi\,{\rm{d}}\eta \geq \textstyle\int_{D_{\triangle \theta}}\mu(\xi,\eta)\,{\rm{d}}\xi\,{\rm{d}}\eta]$ and

$[\textstyle\int\limits_{x_1}^{x_2} p(x,\theta)\,{\rm{d}}x\leq \textstyle\int\limits_{x_1^*}^{x_2^*} p(x,\theta+\triangle \theta)\,{\rm{d}}x.\eqno(24)]$

We have the inequality sign since the area $[D_{\triangle\theta}]$ is larger than D^*, i.e. we also integrate $[\mu(\xi,\eta)]$ over some areas that have been virtually excluded before.

Figure 2
The slice is rotated by $[\triangle \theta]$ . In the new coordinates D becomes $[D_{\triangle \theta}]$ . D^* projects on the same segment [x₁^*,x₂^*] of [-R,R] as $[D_{\triangle \theta}]$ does.

In this paper we consider only the first three requirements, i.e. smoothness, small errors and evenness.

3.2. Finite dimensional approximation

Above we have considered the case of continuous functions, i.e. when we can measure $[p(x,\theta)]$ for any x and θ. This is close to the case when pixels of the sensor have very small size, there is no space between the pixels, we acquire a lot of projections and the acquisition angles are uniformly separated. Here we consider a discrete problem. Let n_x be the number of pixels on a row of the sensor and $[n_\theta]$ be the number of acquisition angles. There are two grids: a grid {x_j}₁^n_x, such that x_j = x₁+\triangle x( j-1), j = $[1,\ldots,n_x]$ , and a grid $[\{\theta_i\}_1^{n_\theta}]$ , $[\theta_i]$ = $[\triangle \theta(i-1)]$ , i = $[1,\ldots,n_\theta]$ , $[\triangle \theta]$ = $[\Theta/(n_\theta-1)]$ . Define matrices p_ij = $[p(x_j,\theta_i)]$ , $[\tilde p_{ij}]$ = $[\tilde p(x_j,\theta_i)]$ , $[\varepsilon_{ij}]$ = $[\varepsilon(x_j,\theta_i)]$ and a vector q_j = q(x_j). Then equation (4) becomes

$[\tilde p_{ij} = p_{ij}\,+\,q_j\,+\,\varepsilon_{ij}, \qquad i = 1,2,\ldots,n_\theta,\quad j = 1,2,\ldots,n_x.\eqno(25)]$

3.2.1. Smoothness

Let there be a function f(x) and an arbitrary point x = a. Then the Taylor polynomial of f(x) about x = a of degree n be denoted by

$[\eqalignno{P_{n,a}(x)={}&f(a)+f'(a)(x-a)+f''(a){{(x-a)^2}\over{2!}}+\ldots \cr& + f^{(n)}(a) {{(x - a)^n}\over{n!}}.&(26)}]$

If the (n+1)-derivative of f(x) exists on the segment (a,x), then according to Taylor's theorem (see, for example, Hirst, 2006 ),

$[f(x)=P_{n,a}(x)+E_n(x),\eqno(27)]$

the error (or remainder) term E_n(x) is given by

$[E_n(x)={{f^{(n+1)}(c)}\over{(n + 1)!}}(x - a)^{(n+1)},\eqno(28)]$

where $[c\in (a,x)]$ .

Let us consider the ith row. We omit its index in p_ij, i.e. by p_j we mean p_ij. Using the Taylor series we can write down

$[p_{j+\triangle j} = p_j + p'_j\triangle j \triangle x + p''_j{{(\triangle j \triangle x)^2}\over{2}} + p'''_j{{(\triangle j \triangle x)^3}\over{6}} + \ldots,\eqno(29)]$

where [p'_j] , [p''_j] , [p'''_j] are the first, second and third derivatives of $[p(x,\theta_i)]$ at x_j. Since the function $[p(x,\theta)]$ is assumed to be smooth, then using formulas similar to (11) we could estimate the remainder (28) for $[p(x,\theta_i)]$ . Therefore p_j could be approximated by a linear combination of $[p_{j\pm1}]$ , $[p_{j\pm2}]$ , $[p_{j\pm3}\ldots]$ . Let 2m be a number of points outside the point x_j (m from the left and m from the right). We want to find coefficients a₁, $[a_{2},\ldots a_{m}]$ such that

$[p_j\simeq\sum_{l=1}^ma_l\left(\,p_{j-l}+p_{j+l}\right).\eqno(30)]$

For each $[p_{j\pm l}]$ we use the Taylor series (29) and find a_l in order to set to zero the maximal number of coefficients of derivatives [p'_j] , $[p''_j, \ldots]$ . For m = 1, we get a₁ = 1/2; for m = 2, a₁ = 2/3, a₂ = −1/6; for m = 3, a₁ = 3/4, a₂ = −3/10, a₃ = 1/20; for m = 4, a₁ = 4/5, a₂ = −2/5, a₃ = 4/35, a₄ = −1/70.

Then the condition of smoothness of $[p(x,\theta)]$ leads to the minimization of

$[\textstyle\sum\limits_{j=m+1}^{n_x-m}\left[p_j-\textstyle\sum\limits_{l=1}^ma_l\left(\,p_{j-l}+p_{j+l}\right)\right]^2.\eqno(31)]$

If we introduce a vector b = $[(a_m,\ldots,a_1,-1,a_1,\ldots,a_{m})]$ , then we can rewrite the last formula as

$[\textstyle\sum\limits_{j=0}^{n_x-2m-1}\left[\sum\limits_{l=1}^{2m+1}b_l\,p_{j+l}\right]^2.\eqno(32)]$

Since the remainder defined in (28) could be estimated, then the constant $[\tilde h\,\,\gt\,\,0]$ in

$[\left|\textstyle\sum\limits_{l=1}^{2m+1}b_l\,p_{j+l}\right|\leq\tilde h \Big(\textstyle\sum\limits_{j=1}^{n_x}p_j^2\Big)^{1/2}\eqno(33)]$

could also be found. Instead of p_j we should substitute $[\tilde p_j-q_j]$ , then find the sum of similar functions for all rows. As the result we get

$[q^TH_1q+2c_1^T q+w_1,\eqno(34)]$

where q = $[(q_1,\ldots,q_{n_x})]$ , w₁ is a constant and T means the transpose of a vector. Equation (34) can be considered as $[\|A_hz-u_\delta\|^2]$ in (7) and errors similar to h and $[\delta]$ in (6) could also be found owing to (33). Note that the matrix H₁ is symmetric, positive-semidefinite (i.e. $[q^TH_1q\geq 0]$ for any q), (4m+1)-diagonal. For instance, for m = 1, H₁ equals

$[{{n_\theta}\over{4}}\left(\matrix{ \quad\!1 & -2 & \quad\!1 & \quad\!0 & \quad\!\ldots & \quad\!0 & \quad\!0 & \quad\!0\cr -2 & \quad\!5 & -4 & \quad\!1 & \quad\!\ldots & \quad\!0 & \quad\!0 & \quad\!0\cr \quad\!1 & -4 & \quad\!6 & -4 & \quad\!\ldots & \quad\!0 & \quad\!0 & \quad\!0\cr \quad\!0 & \quad\!1 & -4 & \quad\!6 & \quad\!\ldots & \quad\!0 & \quad\!0 & \quad\!0\cr \quad\!\vdots & \quad\!\vdots & \quad\!\vdots & \quad\!\vdots & \quad\!\ddots & \quad\!\vdots & \quad\!\vdots & \vdots\cr \quad\!0 & \quad\!0 & \quad\!0 & \quad\!0 & \quad\!\ldots & \quad\!6 & -4 & \quad\!1\cr \quad\!0 & \quad\!0 & \quad\!0 & \quad\!0 & \quad\!\ldots & -4 & \quad\!5 & -2\cr \quad\!0 & \quad\!0 & \quad\!0 & \quad\!0 & \quad\!\ldots & \quad\!1 & -2 & \quad\!1 }\right).\eqno(35)]$

3.2.2. Small perturbations

If we suppose that among all possible solutions $[p(x,\theta)]$ we should find such a solution that the corresponding function q(x) has small absolute values, then on the set of the solutions we have to find an element that minimizes the norm $[\|q(x)\|^2_{L^2}]$ , which may be approximated by

$[\gamma_0\textstyle\sum\limits_{j=1}^{n_x}q_j^2=q^TH_{20}q,\eqno(36)]$

where $[\gamma_0\,\,\gt\,\,0]$ and H₂₀ = $[\rm{diag}\,\{\gamma_0,\ldots,\gamma_0\}]$ is a diagonal matrix. The meaning of $[\gamma_0]$ will be discussed in §3.3. Note the matrix is positive-definite, i.e. $[ q^TH_{20}q\,\,\gt\,\,0]$ for any vector $[q\not\equiv0]$ .

If one also requires small values of [q'(x)] or higher derivatives, then one has to minimize

$[q^TH_{2}q,\eqno(37)]$

where H₂ = $[\textstyle\sum_{l=0}^rH_{2l}]$ , r is the maximal order of the derivatives. Note that matrices H_2l for $[l\,\,\gt\,\,0]$ are only positive-semidefinite. For instance, for l = 1 the matrix H₂₁ may be written as

$[\gamma_1 \left(\matrix{\quad\!1&-1&\quad\!0&\quad\!\ldots&\quad\!0&\quad\!0&\quad\!0\cr -1&\quad\!2&-1&\quad\!\ldots&\quad\!0&\quad\!0&\quad\!0\cr \quad\!0&-1&\quad\!2&\quad\!\ldots&\quad\!0&\quad\!0&\quad\!0\cr \quad\!\vdots&\quad\!\vdots&\quad\!\vdots&\quad\!\ddots&\quad\!\vdots&\quad\!\vdots&\quad\!\vdots\cr \quad\!0&\quad\!0&\quad\!0&\quad\!\ldots&\quad\!2&-1&\quad\!0\cr \quad\!0&\quad\!0&\quad\!0&\quad\!\ldots&-1&\quad\!2&-1\cr \quad\!0&\quad\!0&\quad\!0&\quad\!\ldots&\quad\!0&-1&\quad\!1}\right)\eqno(38)]$

or as (35). If the coefficient $[\gamma_0\,\,\gt\,\,0]$ , then the matrix H₂ is positive-definite as a sum of the positive-definite H₂₀ and positive-semidefinite H_2l, $[l\,\,\gt\,\,0]$ .

3.2.3. Evenness

Since real problems always have errors in data, then requirements (20) and (21) should be changed to the minimization of the corresponding functions in the equations. Let us consider (20). If there is $[j^*\in\overline{1,n_x}]$ such that x_j^* = 0, i.e. the axis of rotation is projected in this point, then define n^* = min( j^*-1,n_x-j^*) and obtain the finite dimensional approximation of (20),

$[\textstyle\sum\limits_{j=j^*-n^*}^{j^*+n^*}\left(\,p_{0j}-p_{n_\theta,2j^*-j}\right)^2.\eqno(39)]$

Now consider the case when there are $[j^*\in\overline{1,n_x-1}]$ and $[\Delta^*\in(0,\triangle x)]$ such that $[x_{j^*}+\Delta^*]$ = 0. Define $[\tau]$ = $[\Delta^*/\triangle x]$ and n^* as the maximal number such that $[x_1\leq -n^*\triangle x]$ and $[x_{n_x}\geq (n^*+1)\triangle x]$ . Then using linear interpolation on a segment [x_j,x_j+1] we get

$[\textstyle\sum\limits_{j = j^*-n^*}^{j^*+n^*}\left[(1-\tau)\left(\,p_{0j}-p_{n_\theta,2j^*-j}\right)+\tau\left(\,p_{0,j+1}-p_{n_\theta,2j^*-j+1}\right)\right]^2.\eqno(40)]$

In a similar way equation (21) may be rewritten. In both cases we get a quadratic function

$[q^TH_3 q + 2 c_3^T q + w_3,\eqno(41)]$

where H₃ is a positive-semidefinite matrix.

3.2.4. Other possible requirements

Let there be some areas where the sinogram should have known values. The simplest way is to find an average value of all elements $[\tilde p_{ij}]$ of a jth column that are inside the area, then subtract the derived value and add the known value to the whole column. However, because of the noise in real data the number of these elements should be large enough in order to decrease the influence of the noise.

The presence of the noise should also be taken into account for equation (24), which for instance may be approximated as

$[\textstyle\sum\limits_{j=j_1}^{j_2}p_{ij}\leq\sum\limits_{j=j_1^*}^{j_2^*} p_{{i^*}j}+C,\eqno(42)]$

where C is a constant and j₁^*,j₂^* depend on j₁, j₂ and $[\triangle \theta]$ .

3.3. The functional to be minimized

Assuming smoothness and evenness of the function $[p(\theta,x)]$ as a function of x we get the following functional to be minimized,

$[\eqalignno{J(q)&\equiv q^THq+2cq+w\cr&\equiv \left(q^TH_1q+2c_1^Tq+w_1\right)+\gamma\left(q^TH_3 q + 2 c_3^T q + w_3\right),&(43)}]$

where $[H\equiv H_1+H_3]$ , $[c\equiv c_1+c_3]$ , $[w\equiv w_1+w_3]$ , $[\gamma\geq 0]$ . Note that H is a positive-semidefinite matrix, thus the solution may be non-unique. To choose a unique solution we suppose that the errors are small and introduce

$[\Omega(q)\equiv q^TH_{2}q.\eqno(44)]$

Then the modified smoothing functional can be written as

$[\eqalignno{M^\alpha[q]&\equiv J(q)+\alpha \Omega(q)\cr&= q^T(H+\alpha H_2)q+2c q+w,&(45)}]$

where $[\alpha\,\,\gt\,\,0]$ is the regularization parameter. Then matrix $[H+\alpha H_2]$ is positive-definite. So if there are no other restrictions on q [such as is defined by (42)], then the functional $[M^\alpha[q]]$ always has a unique solution. In the case of other restrictions, a unique solution exists if the restrictions define a non-empty set.

Note that matrices H₁ and H₂ are band matrices, i.e. their non-zero entries are confined to a diagonal band, comprising the main diagonal and zero or more diagonals on either side. Thus if $[\gamma]$ = 0 one can use special mathematical algorithms to find the minimiser of $[M^\alpha[q]]$ . The matrix H₃ is a sparse matrix but almost all its diagonals have at least one non-zero element. Hence for $[\gamma\,\,\gt\,\,0]$ we are restricted to use general numerical methods to minimize $[M^\alpha[q]]$ .

To find the solution we solve a quadratic programming problem using the C-library implemented by Numerical Algorithms Group (NAG). The solution is based on an inertia-controlling method that maintains a Cholesky factorization of the reduced Hessian and is described by Gill et al. (1991 ).

4. Applications

4.1. Experimental set-up

To illustrate the new method we apply it to two data sets acquired at the Daresbury synchrotron X-ray source (station 16.3; second-generation synchrotron) using a PCO-4000 CCD camera (resolution: 4000×2672; pixel size: 9 µm × 9 µm; dynamic range: 14 bit, Kodak image sensor KAI-11002). The following two samples were used:

(i) Graphite, a cylindrical piece of graphite (diameter 5 mm) with small metallic particles inside, exposure time 150 ms, 3600 projections, CdWO₄ scintillator (500 µm thick), a pixel corresponds to a 2.25 µm × 2.25 µm area.

(ii) Bone, a piece of human bone with tissue all embedded in epoxy (size about 4 mm), exposure time 16 s, 1800 projections, YAG:Ce scintillator (35 µm thick), a pixel corresponds to 0.9 µm × 0.9 µm area.

For the two data sets, 10 dark and 20 flat/white-field images were taken (10 dark and 10 flat-field images before and 10 flat-field images after the acquisition of ordinary projections). Their averaged images were obtained for the standard flat-field correction procedure. The range of acquisition angles for both samples was $[180^\circ]$ .

As mentioned above, the method is strictly applicable only for a monochromatic spectrum of X-rays. However, if attenuation coefficients of elements in a sample vary smoothly over the range of wavelengths representative of the white beam then the method can also be applied. Of course some possible beam hardening effects cannot be fully excluded. We have scanned the two samples in white beam.

After all images (dark/flat fields and ordinary projections) were taken, we applied the standard flat-field correction technique (1) to each projection. Then instead of projections, i.e. vertical slices, we recovered horizontal slices. Each horizontal slice was preprocessed by the method proposed in this paper and then reconstructed using a filtered backprojection algorithm (FBP) for a parallel-beam geometry [see the description by Natterer & Wübbeling (2007)]. Each slice was reconstructed independently of the other slices. Therefore the preprocessing and reconstruction can be done in parallel on a cluster of workstations/desktops. In order to process a 4000×4000 sinogram it usually takes 1–2 s [Intel processor E7200 (dual core), NAG C-library]. Since an optimization procedure is used, the preprocessing time is a function of the initial estimate of the solution, which in our case we set to zero.

4.2. Graphite

We have applied the standard flat-field correction technique to the graphite sample; see the sinogram and the reconstructed slice in Figs. 3(a) and 3(c). Note that the choice of the regularization parameter $[\alpha\,\,\gt\,\,0]$ is very important. If $[\alpha]$ is too large, then the suppression of ring artefacts will be minimal, since $[\|q(x)\|]$ may have only small values. If $[\alpha]$ is too small, then q(x) may have very large values and therefore may change the whole sinogram; see Figs. 3(e) and 3(f), where additional `waves' are seen on both the sinogram and the reconstructed slice (the sinogram tends to flatten in order to avoid possible jumps near the edges).

Figure 3
Sinograms and the reconstructed slices for graphite: (a) and (c) based on the standard flat-field correction; (b) and (d), (g) the proposed method is applied and the `optimal' $[\alpha]$ = 0.00025 is used (only smoothness is assumed); (e) and (f) the proposed method is applied, $[\alpha]$ = 10^-9, $[\gamma]$ = 0.2. The dashed line shows the position of the 150th row, i.e. $[\theta_{150}]$ = 7.5°.

In Fig. 4 several profiles are shown for the 150th row, i.e. for the row that corresponds to $[\theta]$ = 7.5°. One can see that when there is no suppression the profile is close to zero outside the sample's shadow, i.e. for x < 1.6 and x > 6.7 mm. To compare profiles found after the proposed method has been applied we also have applied the mean filter in order to avoid fluctuations. The mean filter has not changed the main properties of the profiles and the profiles are easier to see in Fig. 4(b).

Figure 4
The attenuation p(x) for the 150th row ( $[\theta_{150}]$ = 7.5°) of the sinogram: (a) corresponds to the row shown by the dashed line in Fig. 3(a)

, the standard flat-field correction; (b) the proposed method is applied, then the mean filter (radius is 10 pixels) is applied to exclude fluctuations, the optimal regularization parameter $[\alpha]$ = 0.00025, $[\gamma]$ = 0.0, corresponds to Fig. 3(b)

(black solid line), $[\alpha]$ = 10^-9, $[\gamma]$ = 0.2, corresponds to Fig. 3(e)

(grey solid line), $[\alpha]$ = 10^-9, $[\gamma]$ = 0.0 (black dotted line).

If $[\alpha]$ is small (e.g. $[\alpha]$ = 10^-9), then the profile flattens. However, the shape depends on other parameters. If evenness is used, then the profile tends to have a symmetrical slope outside the sample's shadow; see the grey line in Fig. 4(b) for the segments $[x\in[0,1.6]]$ and $[x\in[6.7,7.8]]$ . Note that evenness has not been used for $[x\,\,\gt\,\,7.8]$ , since the centre of rotation (x = 3.90 mm) is not at the centre (x = 4.27 mm) of the sinogram, so there are constant values for p(x) for x > 7.8 mm.

If $[\alpha]$ is large, then $[\|q(x)\|]$ may have only small values, so that while large jumps for q(x) are possible the number of these jumps is very small. For example, after the mean filter has been applied there is no visible difference between the profiles found after applying the standard flat-field correction or the proposed method ( $[\alpha]$ = 0.00025); the profile is shown Fig. 4(b) by the black solid line.

The best way to find $[\alpha]$ is to find a root of the generalized discrepancy $[\rho(\alpha)]$ . However, it is often difficult to make correct estimations of the errors h and $[\delta]$ . For some cases we may use a priori information about the sinogram. For instance, let us consider a segment of a row where p(x) should be zero [e.g. the left hand 30% of the central row ( $[\theta]$ = 90°) in Fig. 3(a)]. For a given $[\alpha]$ we define $[\chi]$ as the standard deviation of the values of p(x) in this segment (see Fig. 5). For small values of $[\alpha]$ , $[\chi]$ is large as p(x) is not constant on average as already pointed out; see region x < 1.6 mm in Fig. 4 (grey and dotted lines). As the slope of p(x) gets larger towards smaller $[\alpha]$ , $[\chi]$ is a decreasing function for small values of $[\alpha]$ . Note that very small values of $[\alpha]$ should be avoided, since the optimization problem becomes very unstable and the solution may be non-unique. As $[\alpha]$ becomes large the suppression of ring artefacts reduces and $[\chi]$ becomes an increasing function of $[\alpha]$ . It is clear that $[\chi(\alpha)]$ tends to a constant as $[\alpha\to +\infty]$ (see Fig. 5), since there is no suppression at all. Consequently there is an `optimal' value of $[\alpha]$ , where $[\chi(\alpha)]$ has a minimal value. For our case we get $[\alpha]$ = 0.00025 (ln0.00025 ≃ −8.29). The corresponding sinogram and the reconstructed slice are shown in Figs. 3(b) and 3(d).

Figure 5
The standard deviation $[\chi]$ of the values p(x) as a function of $[\ln\alpha]$ (graphite); see text.

4.3. Bone

The method described above can also be a first step in the suppression of more complicated ring artefacts. For instance, we have scanned a piece of a human bone. Several regions on the projections change their values over time [i.e. q depends on time, so we have q(x,t)]. From a physical point of view it can be explained in the following way. Absorption properties of some dust/dirt particles and defects inside the scintillator depend on its temperature. X-rays heat the scintillator, therefore there is a cooling system connected to the scintillator. However, the system is turned on only when the temperature of the scintillator reaches some level and it takes some time to cool the scintillator to the initial temperature. Therefore one can see almost periodical changes of q.

The original sinogram and the corresponding reconstructed image are shown in Figs. 6(a) and 6(c). Applying the above method we obtain Figs. 6(b) and 6(d). There are still some ring artefacts on the image; however, these artefacts are caused by the time-dependent effects explained above. Unfortunately these artefacts do not allow us to estimate errors in input data correctly and therefore find an `optimal' regularization parameter $[\alpha]$ , so we simply tried $[\alpha]$ = 0.00001.

Figure 6
A piece of a human bone: (a) and (b) sinograms, (c) and (d) reconstructed images; (a) and (c) the standard flat-field correction is applied, (b) and (d) the proposed method is used ( $[\alpha]$ = 0.00001).

The correction of the remaining irregular ring artefacts requires changes in the experimental procedure to avoid thermal effects by (i) changing a scintillator, (ii) moving a sensor as in Davis & Elliott (1997), (iii) acquiring more flat-field images during the acquisition. An example of a method correcting irregular ring artefacts caused by vibrations in a set-up can be found by Titarenko et al. (2010b).

5. Discussion

The primary aim of the paper is to propose the new method. Here we compare it with the method of Boin & Haibel (2006) as implemented within the High Speed Tomography Reconstruction (PyHST) reconstruction software package (https://www.esrf.eu/UsersAndScience/Experiments/TBS/SciSoft ), a very popular software tool used at the ESRF (European Synchrotron Radiation Facility, Grenoble, France). To get the best images we have tuned the parameter N, the so-called span factor (see Boin & Haibel, 2006). The results of the comparison are shown in Fig. 7. One can see that for these data sets at least a better suppression of ring artefacts is achieved by the newly proposed method.

Figure 7
The pieces of graphite and human bone reconstructed by the proposed method (a) and (c); and using the PyHST software (b) (N = 51) and (d) (N = 171).

Note that a simplified version of this algorithm can be found by Titarenko et al. (2010a ), where only the smoothness of $[p(x,\theta)]$ and the small values of q(x) were used. There it was required that $[\textstyle\sum_{j=1}^{n_x-1}(\,p_j-p_{j+1})^2]$ and $[\textstyle\sum_{j=1}^{n_x}q_j^2]$ should be small in order that $[H+\alpha H_2]$ from (45) was a tridiagonal symmetric matrix, so a fast (analytical) method of finding q could be applied. The method proposed here is more general, since higher derivatives of $[p(x,\theta)]$ and q(x) can be used and a condition of evenness has also been introduced.

In the future more sophisticated a priori restrictions are expected to be proposed, some of which will apply to particular sample types of acquisition geometries. Nevertheless, simpler methods such as that proposed by Titarenko et al. (2010a) are also needed because they can be used in cases where fast algorithms are needed to suppress artefacts.

A priori information has been used previously by Titarenko et al. (2009), but it requires knowledge of the attenuation coefficients in some areas in order to apply a correction to a wider area. In the current paper no information about these values is used, so that the method proposed here is applicable to a wider range of complex samples which gives the better suppression in cases where Titarenko et al. (2009) can be applied.

It was mentioned in §4.3 that in some cases ring artefacts can have irregular shapes such that those described by Titarenko et al. (2010b) caused by various time-dependent vibrations in the experimental set-up. In such cases it may be better to transform these irregular artefacts to regular ones first by shifting a reference flat field and then only correct by the current method. Unfortunately, this approach is not applicable to all possible irregular ring artefacts.

Now let us discuss possible bottlenecks when applying the method proposed here to some real data set. For the first group of samples there is an area near the centre of the rotation where the attenuation coefficient differs sufficiently compared with other areas of the slice. For instance it may be a metal particle, a fibre crossing the slice, a small pore in a foam/bone or a cell in a plant. As a result the corresponding sinogram will often have a jump near the central vertical line. For the second group there is an axial symmetric feature, e.g. a thin glass tube or container protecting a biological sample. In both cases it is difficult to decide whether we suppress an artefact or a real feature. However, the following two approaches may help to find an answer. Firstly, try to consider neighbouring slices as dependent, i.e. their attenuation coefficients should not vary very rapidly from one slice to another and therefore $[p(x,z,\theta)]$ or $[\textstyle\int_{0}^{180^\circ}\!\!p(x,z,\theta)\,{\rm{d}}\theta]$ should have similar values if points (x₁,z₁) and (x₂,z₂) are close. As a result you may get an algorithm that will suppress the ring artefacts for the whole sample at once. One example of such algorithms can be found by Titarenko & Yagola (2010 ). Secondly, modify the data acquisition procedure. For example, take at least one additional projection at the end of acquisition just shifting a sample perpendicular to the beam. In this case if your vector or matrix q is found correctly, then the two projections acquired before and after the shift will be the same after correction; otherwise you may find the points where the projections differ, and correct the corresponding elements of vector q. Ideally, it would be better to have an adaptive method to correct ring artefacts, i.e. shift, rotate or incline a sample depending on results of the preprocessing. For instance, it can be done in the following way: acquire only odd projections, then use them to find the vector q, at the same time continue acquiring even projections, so at the end of acquisition you may check whether your vector q is found correctly; if not, acquire additional projections.

6. Conclusions

A new method for suppressing ring artefacts has been proposed which is applied directly to sinograms prior to image reconstruction. Generally this is quicker than approaches that must be applied to reconstructions. It is based on ideas and methods from the theory of inverse and ill-posed problems. The new technique is applied to two data sets obtained using a parallel white-beam synchrotron X-ray illumination and is found to work very well. This can be a first step towards the suppression of more complicated ring artefacts when one removes regular artefacts first and then applies another technique to remove remaining irregular artefacts. In principle the method can also be extended to fan-beam reconstructions.

Acknowledgements

The authors would like to thank STFC for funding this project and Daresbury Laboratory (Warrington, UK) for the use of station 16.3 to image the graphite and the bone samples. The authors are thankful to Alessandro Mirone for his help with the PyHST software. ST would like to thank the EPSRC for funding her visit to Manchester through a `Collaborating for Success' grant.

References

Antoine, C., Nygård, P., Gregersen, Ø. W., Holmstad, R., Weitkamp, T. & Rau, C. (2002). Nucl. Instrum. Methods Phys. Res. A, 490, 392–402. Web of Science CrossRef CAS Google Scholar
Bakushinskii, A. B. (1984). Comput. Math. Math. Phys. 24, 1258–1259. CrossRef Web of Science Google Scholar
Bakushinsky, A. B. & Kokurin, M. Yu. (2004). Iterative Methods for Approximate Solution of Inverse Problems, Mathematics and its Applications, Vol. 577. Dordrecht: Springer. Google Scholar
Banhart, J. (2008). Advanced Tomographic Methods in Materials Research and Engineering, Monographs on the Physics and Chemistry of Materials, No. 66. Oxford University Press. Google Scholar
Boin, M. & Haibel, A. (2006). Opt. Express, 14, 12071–12075. Web of Science CrossRef PubMed Google Scholar
Davis, G. R. & Elliott, J. C. (1997). Nucl. Instrum. Methods Phys. Res. A, 394, 157–162. CrossRef CAS Web of Science Google Scholar
Davidson, D. W., Fröjdh, C., O'Shea, V., Nilsson, H.-E. & Rahman, M. (2003). Nucl. Instrum. Methods Phys. Res. A, 509, 146–150. Web of Science CrossRef CAS Google Scholar
Engl, H. W., Hanke, M. & Neubauer, A. (1996). Regularization of Inverse Problems, Mathematics and its Applications, Vol. 375. Dordrecht: Kluwer. Google Scholar
Gill, P. E., Murray, W., Saunders, M. A. & Wright, M. H. (1991). SIAM Rev. 33, 1–36. CrossRef Web of Science Google Scholar
Groetsch, C. W. (1993). Inverse Problems in the Mathematical Sciences. Braunschweig: Vieweg. Google Scholar
Hadamard, J. (1923). Lectures on Cauchy's Problem in Linear Partial Differential Equations. New Haven: Yale University Press. Google Scholar
Hansen, P. C. (1992). Inv. Probl. 8, 849–872. CrossRef Web of Science Google Scholar
Hirst, K. E. (2006). Calculus of One Variable. London: Springer. Google Scholar
Ivanov, V. K., Vasin, V. V. & Tanana, V. P. (2002). Theory of Linear Ill-Posed Problems and its Applications. Utrecht: VSP. Google Scholar
Kaltenbacher, B., Neubauer, A. & Scherzer, O. (2008). Iterative Regularization Methods for Nonlinear Ill-Posed Problems, Radon Series on Computational and Applied Mathematics. Berlin: Walter de Gruyter. Google Scholar
Kowalski, G. (1977). IEEE Trans. Nucl. Sci. 24, 2006–2016. CrossRef Web of Science Google Scholar
Lavrentiev, M. M., Avdeev, A. V. & Lavrentiev-jun, M. M. (2003). Inverse Problems of Mathematical Physics. Utrecht: VSP. Google Scholar
Mahajan, V. N. (2001). Optical Imaging and Aberrations, Part II, Wave Diffraction Optics. SPIE Press. Google Scholar
Natterer F. & Wübbeling, F. (2007). Mathematical Methods in Image Reconstruction, SIAM Monographs on Mathematical Modeling and Computation, Vol. 5. Philadelphia: SIAM. Google Scholar
Raven, C. (1998). Rev. Sci. Instrum. 69, 2978–2980. Web of Science CrossRef CAS Google Scholar
Sijbers, J. & Postnov, A. (2004). Phys. Med. Biol. 49, N247–N253. Web of Science CrossRef PubMed Google Scholar
Stock, S. R. (2008). Micro-Computed Tomography: Methodology and Applications. CRC Press. Google Scholar
Tikhonov, A. N. (1963). Sov. Math. Dokl. 5, 1035–1038. Google Scholar
Tikhonov, A. N. & Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. New York: John Wiley and Sons. Google Scholar
Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V. & Yagola, A. G. (1995). Numerical Methods for the Solution of Ill-Posed Problems, Mathematics and its Applications, Vol. 328. Dordrecht: Kluwer. Google Scholar
Tikhonov, A. N., Leonov, A. S. & Yagola, A. G. (1998). Nonlinear Ill-Posed Problems, Vols. 1 and 2, Applied Mathematics and Mathematical Computation, Vol. 14. London: Chapman and Hall. Google Scholar
Titarenko, S. S. & Yagola, A. G. (2010). Moscow Univ. Phys. Bull. 65, 65–67. Web of Science CrossRef Google Scholar
Titarenko, S., Withers, P. J. & Yagola, A. (2010a). Appl. Math. Lett. Submitted. Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. J. (2009). Appl. Phys. Lett. 95, 071113. Web of Science CrossRef Google Scholar
Titarenko, V., Titarenko, S., Withers, P. J., De Carlo, F. & Xiao, X. (2010b). J. Synchrotron Rad. Submitted. Google Scholar
Vidal, F. P., Létang, J. M., Peix, G. & Cloetens, P. (2005). Nucl. Instrum. Methods Phys. Res. B, 234, 333–348. Web of Science CrossRef CAS Google Scholar
Wahba, G. (1977). SIAM J. Numer. Anal. 14, 651–667. CrossRef Web of Science Google Scholar
Walls, J. R., Sled, J. G., Sharpe, J. & Henkelman, R. M. (2005). Phys. Med. Biol. 50, 4645–4665. Web of Science CrossRef PubMed Google Scholar
Yagola, A. G., Leonov, A. S. & Titarenko, V. N. (2002). Inv. Probl. Eng. 10, 117–129. Web of Science CrossRef Google Scholar
Yang, J., Zhen, X., Zhou, L., Zhang, S., Wang, Z., Zhu, L. & Lu, W. (2008). The 2nd International Conference on Bioinformatics and Biomedical Engineering, pp. 2386–2389. Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 17| Part 4| July 2010| Pages 540-549

doi:10.1107/S0909049510010964