Suppression of ring artefacts when tomographing anisotropically attenuating samples

Titarenko, S.; Titarenko, V.; Kyrieleis, A.; Withers, P.; De Carlo, F.

doi:10.1107/S0909049511006005

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 18| Part 3| May 2011| Pages 427-435

doi:10.1107/S0909049511006005

Suppression of ring artefacts when tomographing anisotropically attenuating samples

Sofya Titarenko,^a Valeriy Titarenko,^b ^* Albrecht Kyrieleis,^b Philip J. Withers ^b and Francesco De Carlo ^c

^aDepartment of Mathematics, Faculty of Physics, Moscow State University, Leninskie Gory, Moscow 119991, Russia, ^bHenry Moseley X-ray Imaging Facility, School of Materials, University of Manchester, Grosvenor Street, Manchester M1 7HS, England, and ^cX-ray Science Division, Argonne National Laboratory, 9700 South Cass Avenue, Bldg 401, Argonne, IL 60439-4837, USA
^*Correspondence e-mail: valeriy.titarenko@manchester.ac.uk

(Received 16 June 2010; accepted 17 February 2011; online 2 April 2011)

There are many objects for which the attenuation varies significantly as they are rotated during computerized X-ray tomography, for example plate samples. This can lead to significant ring artefacts in the subsequent tomographic reconstructions. In this paper a new method is presented that can successfully suppress such ring artefacts and is applicable to both parallel and cone-beam geometries. Rapid correction is achieved via an analytical formula which involves only a matrix-vector multiplication, for which the matrix is known and depends on a regularization parameter. The efficacy of the method is demonstrated for a paleontological sample (calcified shark cartilage) and a carbon–carbon composite/Ti–SiC metal matrix composite test sample.

Keywords: X-ray absorption computerized tomography; laterally extended objects; ring artefacts; filtered back-projection; inverse problem; laminography.

1. Introduction

There are many objects which attenuate X-rays very differently in different directions. For example, this may be because geometrically they have a large aspect ratio, such as a plate, and so attenuate significantly more when viewed along the major axis, or because they comprise differently attenuating phases that are not homogeneously distributed. As a result when acquiring a series of radiographs (projections) as the sample is rotated in order to make a tomograph, the details in some projections, or regions within certain projections, will be only faintly represented, if at all. This can lead to significant ring artefacts in the subsequent reconstructions which vary in strength around the rings. In this paper a method is presented that can successfully suppress such ring artefacts.

Let us consider the following typical experimental set-up for X-ray tomography on a synchrotron X-ray beamline or laboratory X-ray imager. There is a source of X-rays which provides a beam of given geometry. The geometry of the beam may be arbitrary and is not restricted to classical parallel or cone-beam geometries. The only assumption is that the intensity of the beam is unchanged during the scan. The beam passes through a sample that is rotated about a vertical axis before falling onto a scintillator which absorbs the photon energy and re-emits it as visible light in turn. The light then passes through an optical system and is recorded by a CCD (charged couple device) camera. As a result a series of two-dimensional radiographs (projections) is obtained.

Given a set of projections one can apply various reconstruction algorithms in order to find the structure of the sample. For example, an FBP (filtered back-projection) algorithm for parallel geometry of the beam and an FDK (abbreviated from Feldkamp, Davis, Kress) algorithm for cone-beam geometry are currently the most popular algorithms (Natterer & Wübbeling, 2007 ; Kak & Slaney, 1988 ; Buzug, 2008 ). Let there be a monochromatic X-ray beam, having a beam path, L, through the sample, and I₀ be the intensity of the X-ray photons before the sample. Then the intensity of the X-rays after the sample can be written as I = I₀ ·A = I₀exp(-p), where the attenuation factor p is the line integral $[\int_L\mu({\bf r})\,{\rm d}{\bf r}]$ along L and $[\mu({\bf r})]$ is the linear attenuation coefficient at a point $[{\bf r}]$ . This means that p is the path through the sample along the beam trajectory. Let x and y denote the horizontal and vertical coordinates, which are parallel to the rows and columns of the pixels on the CCD sensor, respectively, and t be a moment of time. Then the reconstruction algorithms use the attenuation factor p(x,y,t) in order to find the linear attenuation coefficient $[\mu({\bf r})]$ at a point $[{\bf r}]$ .

In addition to electrons generated by the absorption of visible photons, the CCD sensor is also sensitive to electrons (so-called dark-noise electrons) generated by thermal processes within it. Therefore, to a first approximation the standard flat-field correction procedure can be defined (see, for example, Stock, 2008 ),

$[p = \ln{{W-D} \over {I-D}}, \eqno(1)]$

where W is a white (or flat) field (recorded with the sample removed from the beam) and D is a dark field (recorded with the beam switched off). Unfortunately, for various reasons mentioned below the measured attenuation factor $[\tilde p]$ differs from the real one p, so we can always write

$[\tilde p = p + q. \eqno(2)]$

As a first approximation the term q often does not depend on time, i.e. $[\tilde p(x,y,t) = p(x,y,t)+q(x,y)]$ . This case corresponds, for example, to the case where there are some particles or pieces of dust/dirt on the surface of the scintillator which absorb X-rays. However, some other factors can lead to q being dependent on p or t. For example, the beam may not be monochromatic. Let the beam have an intensity $[I_0(\lambda)]$ dependent on the wavelength $[\lambda]$ , and p and q be also functions of $[\lambda]$ . The measured attenuation can be written as

$[\tilde p^* = -\ln\left\{{{\int I_0(\lambda)\exp[-p(\lambda)-q(\lambda)]\,{\rm d}\lambda} \over {\int I_0(\lambda)\,{\rm d}\lambda}}\right\} \eqno(3)]$

and in the general case it is a function of p^*, q^* and $[I_0(\lambda)]$ , where

$[\eqalignno{ p^* =\, & -\ln\left\{{{\int I_0(\lambda)\exp[-p(\lambda)]\,{\rm d}\lambda} \over {\int I_0(\lambda)\,{\rm d}\lambda}}\right\},\cr q^* = \, & -\ln\left\{{{\int I_0(\lambda)\exp[-q(\lambda)]\,{\rm d}\lambda} \over {\int I_0(\lambda)\,{\rm d}\lambda}}\right\}. &(4)}]$

This means that $[\tilde p^* = \tilde p^*[\,p^*, q^*, I_0(\lambda)]]$ , i.e. q^* is an implicit function of $[\tilde p^*]$ , p^* and $[I_0(\lambda)]$ , and $[\tilde p^* - p^*]$ is only the first term in the Taylor series of q^* for a given $[I_0(\lambda)]$ . Among the other causes of q being a function of p and t there are the following: the scintillator may have various impurities inside (the attenuation properties of these impurities/defects may be temperature dependent, so that intensity variation on passing through the sample as it is rotated may lead to temperature variations on the scintillator); scintillator performance may be non-linear; some pixels of the CCD sensor may have a non-linear response (see Bardelli et al., 2006 ; Moszyński, 2006 ; Krumm et al., 2008 ). In such cases the intensity of ring artefacts will depend on the shape of the sample: a ring artefact is more intense if the attenuation factor p has a greater value because the intensity of X-rays that pass through the sample is an exponential function of the attenuation factor p and the same errors in the measured intensity lead to greater errors in p if the intensity is small. Therefore for samples showing very different levels of attenuation in different directions, i.e. where p varies markedly with projection angle, the ring artefacts may have irregular shapes. See, for example, Fig. 1(b) which relates to a plate sample; the artefacts are most intense for arcs almost parallel to the longest axis of the sample.

Figure 1
The sinogram and reconstructed slices of the calcified shark cartilage (the right-hand column contains magnified regions of images shown in the left-hand column): (a) the sinogram (attenuation factor p) obtained after applying the standard flat-field correction technique; the slice reconstructed using FBP with (c) and without (b) the proposed ring artefact suppression method.

In this paper we try to find a priori information about the attenuation factor p(x,y,t) that will help us to suppress these irregular ring artefacts. I(x,y) is the intensity of the X-ray beam just before the scintillator after the beam has passed through the sample. Let I_v(x,y) be the intensity of the visible light just after the scintillator. Often it can be assumed that the transform $[I\to I_{\rm v}]$ is linear and may be written as a convolution,

$[\eqalignno{ I_{\rm v}(x,y) =\, & K(x,y)*I(x,y)\cr =\, & {\int}\!\!{\int} I(x-\xi,y-\eta)K(\xi,\eta)\,{\rm d} \xi \,{\rm d}\eta, &(5)}]$

where K(x,y) is the point-spread function (PSF). In many applications K(x,y) has an almost Gaussian shape (Banhart, 2008 ; Mahajan, 2001 $[Mahajan, V. N. (2001). Optical Imaging and Aberrations, Part II, Wave Diffraction Optics. Bellingham: SPIE Press.]$ ), i.e.

$[K(x,y) = K_0 \exp[-\gamma(x^2+y^2)], \eqno(6)]$

where $[K_0, \gamma\,\,\gt\,\,0]$ . Since

$[{\partial\over \partial x}K(x,y) = -2K_0\gamma \exp[-\gamma(x^2+y^2)]\eqno(7)]$

and $[0\leq I(x,y)\leq I_{\rm max}]$ , then, owing to the differentiation property of the convolution, we obtain

$[\eqalignno{\left|{\partial\over \partial x} I_{\rm v}(x,y)\right| =\, & \left|I(x,y) {\partial\over \partial x} K(x,y)\right|\cr &\leq 2K_0\gamma I_{\rm max} {\int}\!\!{\int} \left|\xi \exp[-\gamma(\xi^2+\zeta^2)]\right|\,{\rm d}\xi \,{\rm d}\zeta\cr & \leq 2K_0 I_{\rm max} (\pi/\gamma)^{1/2}. &(8)}]$

So it can be assumed that I_v(x,y) is a smooth function. Therefore, we can suppose that p(x,y,t) is also a smooth function of x and y.

Various assumptions on p and q have been made in a series of papers submitted by the authors. In the case of regular ring artefacts two methods have been proposed in Titarenko et al. (2009 , 2010a ). In the former, the correction of ring artefacts is based on a knowledge of the attenuation coefficients in some areas of a reconstructed slice. In this case the difference between the exact image and the image obtained from a measured sinogram is almost a constant over each half-ring (all full-rings are concentric with respect to the centre of rotation of the sample). In the latter, the smoothness of the sinogram along the horizontal coordinate, as well as equivalence between the first and the final rows of the sinogram (they should be the same in the case of 360° rotation or flipped with respect to the axis of rotation in the case of 180° rotation) are used to obtain a quadratic programming problem whose solution is a time-independent function, q(x). In contrast to the former, explicit values of the attenuation coefficients are not used. For some simple, but often encountered, cases an analytical formula has been proposed in Titarenko et al. (2010b ); this formula is also used in the method proposed in this paper. A case of irregular ring artefacts caused by small vibrations in experimental set-up has been considered in Titarenko et al. (2010c ). In the current paper we consider more general irregular ring artefacts.

We think it is impossible in principle to develop an algorithm that suppresses ring artefacts for any dataset. The reason lies in the nature of tomography. Image reconstruction from a given set of projections is an ill-posed problem. This means that at least one of the following three requirements is not satisfied: the solution exists, is unique and stable. Image reconstruction is always an unstable problem, i.e. small errors in input data lead to significant variations in the solution, see examples by Lavrent'ev et al. (2001 ) and Natterer & Wübbeling (2007). Also from the theory of ill-posed problems it is known that without a priori information about properties of a sample and noise levels in input data you may achieve a solution which varies significantly from the exact one. Note that many well known methods of image reconstruction assume some properties of the solution, which are often not met in real problems. In the case of filtered back-projection the solution is to be an infinite number multiplied by the differentiable function, i.e. it cannot be used to reconstruct a sample with any edges/boundaries, even a sphere or cylinder. In addition, the methods have some internal parameters allowing a user to `regularize' a solution in order to avoid high noise level in a solution. This `regularization' is often implemented as cutting high frequencies in the input data.

Similar ideas of `filtering' are used in all the methods we know to suppress ring artefacts, which have been developed over the last 20 years. Of course, this `filtering' is implemented in different ways. There seem to be three main approaches. For the first, a data acquisition is modified in order to `blur' ring artefacts. In Davis & Elliott (1997 ) time-delay integration is used, so a detector is moved laterally during acquisition. This approach allows users to avoid strong rings on a given slice; however, it may also increase the noise level on several neighbouring slices. The second approach is based on post-processing of the reconstructed images (see, for example, Sijbers & Postnov, 2004 ; Walls et al., 2005 ; Yang et al., 2008 ; Titarenko et al., 2009). A polar-to-Cartesian coordinates transform is often used to find the rings. The third approach is based on applying various filtering techniques to the input data (often sinograms) (see Antoine et al., 2002 ; Raven, 1998 ; Boin & Haibel, 2006 ; Titarenko et al., 2010b; Titarenko & Yagola, 2010 ). Of course, some methods may combine several approaches. The method proposed in this paper belongs to the last group.

In our opinion, owing to the ill-posedness of image reconstruction, data processing does not allow ring artefacts to be suppressed without real features being suppressed at the same time. A good example is a sample with a small dense particle near the centre of rotation. In this case a sinogram will have a column which is brighter than the neighbouring ones. The same sinogram may be formed in the case of a ring artefact and no dense particle. Other examples may involve samples with concentric features. It seems the best approach is to combine data processing and acquisition, so that if there are any doubts as to whether we see an artefact it is possible to acquire more data (moving/rotating/inclining a sample) to check this. Of course, it may be possible to develop a ring-suppression method specifically for a particular class of sample such that real features may be easily distinguished. However, in the case of artefacts developed during acquisition or from some non-linear effects, e.g. beam-hardening, it can be very hard to suppress ring artefacts.

2. Ring artefact suppression

Let us take a horizontal slice of data, i.e. the data recorded by a single row of pixels over all projections. This can be represented by the matrix $[\tilde p_{ij} = \tilde p(x_j,y,t_i)]$ , where t_i, $[i = 1,\ldots,m]$ , and x_j, $[j = 1,\ldots,n]$ , are grids (m and n are numbers), i is the projection number and j the pixel in the row. In many cases, e.g. when the sample rotates at constant speed, the variable t (time) may be replaced by the angle of rotation of the sample $[\theta]$ . However, in the case of variations in beam intensity, the temperature of the scintillator or the camera or the strength of the artefacts depends not only on $[\theta]$ so it is better to use t. The grid {x_j}₁ⁿ is uniform (x_j-x_j-1 is often the pixel size), {t_i}₁^m may be arbitrary, since we have not assumed that p(x,y,t) is a smooth function of t. Therefore

$[\tilde p_{ij} = p_{ij}+q_{ij},\eqno(9)]$

where p_ij must be found. Finding p_ij is equivalent to finding q_ij, which we will use thereafter. q denotes the matrix q_ij. Since p(x,y,t) is a smooth function of x, we require that

$[G(q) = \sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n-1} |p_{ij}-p_{i, j+1}|^2 \eqno(10)]$

attains its minimal value for a given error in measurement. Another assumption is that the error in measurements is relatively small, i.e.

$[H(q) = \sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n} |q_{ij}|^2 \eqno(11)]$

should tend to zero.

Let us introduce a smoothing (Tikhonov's) functional

$[M^\alpha[q] = G(q)+\alpha H(q), \eqno(12)]$

where $[\alpha\,\,\gt\,\,0]$ is a regularization parameter. Then the minimizer $[q^\alpha]$ of this equation tends to the exact solution as all errors in the measurements tend to zero. Note that the regularization parameter $[\alpha]$ depends on errors and cannot be chosen arbitrarily. However, for given non-zero errors one could vary $[\alpha]$ in order to satisfy some restrictions on the image structure. More information about methods of solving ill-posed problems can be found in Titarenko et al. (1995 , 1998 ), Engl et al. (1996 ), Bakushinsky & Kokurin (2004 ) and Ivanov et al. (2002 ).

Now we assume that q_ij can be written as a final sum

$[q_{ij} = \sum\limits_{s = 1}^S f_{is}c_{sj}, \eqno(13)]$

where f_is are orthonormal vectors, i.e.

$[\sum\limits_{i = 1}^m f_{is}\cdot f_{iv} = \left\{\matrix{1,& s = v,\cr 0,& s\ne v.}\right.\eqno(14)]$

The vector $[\tilde p_{ij}]$ is defined for $[1\leq j\leq n]$ . Let $[\tilde p_{i0} = \tilde p_{i1}]$ and $[\tilde p_{i,n+1} = \tilde p_{i,n}]$ . Then the vectors c_w and g_w are introduced such that

$[c_w = (c_{w1},c_{w2},\ldots,c_{wn}), \eqno(15)]$

$[g_w = (g_{w1},g_{w2},\ldots,g_{wn}), \eqno(16)]$

with elements defined as

$[g_{wj} = \sum\limits_{i = 1}^m(-\tilde p_{i,j-1}+2\tilde p_{ij}- \tilde p_{i,j+1})f_{iw}.\eqno(17)]$

Now our aim is to find a global minimizer c_sj of the smoothing functional (12). As shown in Appendix A the property (14) allows us to obtain S systems of linear equations (n equations with n variables)

$[Ac_w^T = g_w^T, \eqno(18)]$

where T is a transposed vector, i.e. a column. All these systems have the same matrix A.

To find the solution of (18) we use the analytical formula found in Titarenko et al. (2010b). Then

$[c_{wj} = \sum\limits_{k = 1}^{n}H_{jk}g_k, \eqno(19)]$

where

$[\eqalignno{H_{jk} ={}& {{\tanh(\tau/2)} \over {\alpha\sinh(n\tau)}}\Big\{\cosh\left[(n-|j-k|)\tau\right]\cr& +\cosh\left[(n+1-j-k)\tau\right]\Big\}&(20)}]$

and $[\tau = 2\;{\rm arcsinh}({\alpha^{1/2}}\,2/2)]$ . To avoid large numbers the last formula can be rewritten as

$[\eqalignno{H_{jk} ={}& {{\exp[-(\,j-k)\tau]} \over {[\alpha(\alpha+4)]_{\vphantom{\big{|}}}^{1/2}}}\cr&\times{{\big\{1+\exp[-(2n-2j+1)\tau]\big\}\big\{1+\exp[-(2k-1)\tau]\big\}} \over {1-\exp(-2n\tau)}}\cr&&(21)}]$

for $[j\geq k]$ . The matrix H_jk is symmetric, so $[j\,\,\lt\,\,k]$ is trivial.

For the practical applications described below we use the following set of orthonormal vectors (see Hsu, 1995 ),

$[\eqalignno{ f_{i1}& = {{1} \over m^{1/2}}\cdot\left(1,1,\ldots,1\right),\cr f_{i2}& = \left({2 \over m}\right)^{\!1/2}\cdot\left[\cos\left({{2\pi} \over {m}}\right), \cos\left({{4\pi} \over {m}}\right),\ldots, \cos\left({{2\pi m} \over {m}}\right)\right],\cr f_{i3}& = \left({2 \over m}\right)^{\!1/2}\cdot\left[\sin\left({{2\pi} \over {m}}\right), \sin\left({{4\pi} \over {m}}\right),\ldots, \sin\left({{2\pi m} \over {m}}\right)\right],\cr f_{i4}& = \left({2 \over m}\right)^{\!1/2}\cdot\left[\cos\left({{4\pi} \over {m}}\right), \cos\left({{8\pi} \over {m}}\right),\ldots, \cos\left({{4\pi m} \over {m}}\right)\right],\cr f_{i5}& = \left({2 \over m}\right)^{\!1/2}\cdot\left[\sin\left({{4\pi} \over {m}}\right), \sin\left({{8\pi} \over {m}}\right),\ldots, \sin\left({{4\pi m} \over {m}}\right)\right],\ldots . &(22)}]$

These vectors may be abbreviated to

$[\eqalignno{f_{i1} & = {{1} \over {{m^{1/2}}}},\cr f_{i,2s}& = \left({2 \over m}\right)^{\!1/2}\cos\left({{2\pi i s} \over {m}}\right),\cr f_{i,2s+1}& = \left({2 \over m}\right)^{\!1/2}\sin\left({{2\pi i s} \over {m}}\right).&(23)}]$

Note that there are only m independent vectors, so $[m\leq S]$ .

3. Practical applications

3.1. Test object 1: a high-aspect-ratio (laterally extended) sample

Test object 1 is calcified cartilage from the pectoral girdle of a fossilized symmoriid shark. A series of projections was acquired at the microtomography beamline 2-BM at the Advanced Photon Source (APS) at Argonne National Laboratory, USA. A double multilayer monochromator was used. The X-ray beam had a parallel geometry so a standard filtered back-projection (FBP) algorithm (see e.g. Natterer & Wübbeling, 2007; Kak & Slaney, 1988) was applied to the sinograms preprocessed by the algorithm proposed above. Although the sample comprises phases of similar attenuation coefficients it has a large aspect ratio ( [5.6&emsp14;:&emsp14;1.0] ) in the projections.

From Figs. 1(a) and 1(b) it can be seen that the standard flat-field correction technique is not effective owing to the large aspect ratio. As a result the strength of ring artefacts seen in the reconstructed horizontal slice in Fig. 1(b) varies with angle about the centre of the slice, being largest normal to the long axis of the sample, i.e. the arcs are most significant essentially where they are parallel to the long axis. After applying the proposed ring-suppression method [taking $[\alpha = 0.001]$ , S = 21 and the functions f_is defined in (23)] the artefacts have been suppressed, as demonstrated in Fig. 1(c). Note that weak ring artefacts still persist near the edge of the reconstructed slice. However, they are not intense and can be suppressed further if a larger value of S is used. The value S = 21 means that the minimal period of the orthogonal functions defined in (23) is one-tenth of the height of the sinogram (the height is measured along variable t). Roughly speaking one can expect that artefacts whose period is greater than the minimal period will be suppressed, whereas artefacts with a smaller period may still persist but with decreased strength. The parameter S controls how ring artefacts are suppressed along a radial coordinate (in a polar frame) and $[\alpha]$ along the angular coordinate. S = 21 is taken, since the minimal period of artefacts seems to be about one-tenth of the sinogram's height, and $[\alpha = 0.001]$ seems to be optimal for the suppression along the angular coordinate (for a larger $[\alpha]$ some rings remain while for a smaller $[\alpha]$ some genuine features become suppressed). Note that there is also a `spider' artefact near the centre of Fig. 1(c) which is discussed below.

3.2. Test object 2: compound sample

In order to more seriously challenge the correction algorithm a sample was devised showing more extreme heterogeneity of the attenuation coefficient. In essence the sample comprises a (5.5 mm thick) carbon–carbon composite (CX-270G grade material supplied by Toyo Tanso Ltd, fully graphitized with a nominal density of 1.67 g cm^-3; the material has a flat two-dimensional woven architecture of high modulus fibres) sandwiched between two (2.5 ×1.0 and 6.5 ×1.2 mm) blocks of Ti matrix/unidirectional 140 µm-diameter SiC monofilament composite all wrapped up in (30 µm) aluminium foil. The dataset was acquired at the B-16 beamline of the Diamond Light Source, UK. The X-ray beam had a white spectrum and parallel geometry. FBP was used to reconstruct slices. While the aspect ratio of the sample is relatively small (1.5 : 1.0) the attenuation coefficient is highly anisotropic because the carbon composite at the centre attenuates only slightly whereas the metal composite attenuates the X-ray beam significantly. In addition, the presence of the SiC fibres in the metal composite makes this region locally heterogeneous. Beam hardening arising from the white beam is also apparent as black lines parallel to the metal composite blocks and diagonal lines in Fig. 2, making the treatment of the ring artefacts more complicated.

Figure 2
The compound sample (Ti–SiC composite/CX-270G carbon–carbon composite sandwich wrapped in aluminium foil): (a) the full sinogram; (b), (c) enlarged regions (marked by white boxes) of the sinogram; (d), (e) the standard flat-field correction is applied; (f), (g) the proposed method is applied ( $[\alpha = 0.0001]$ , S = 20).

4. Discussion

In applying our approach a key issue is the choice of $[\alpha]$ and S. First let us discuss how $[\alpha]$ should depend on s. One may use the same $[\alpha]$ for all f_is. However, the function q(x_i,y,t) is assumed to be smooth along t. Note that Fourier coefficients of a smooth function decrease as 1/s (see Beerends et al., 2003 ). Since the functions f_is defined in (23) are discrete Fourier coefficients, the use of $[\alpha_s]$ = $[\alpha_1\cdot s^2]$ would seem to be more reasonable than $[\alpha_s]$ = $[\alpha_1]$ . This strategy can be assessed by comparing images from the right- and left-hand columns of Fig. 3. For small values of S the quality of the reconstruction of the titanium composite blocks is almost the same for both cases, while for large S the case of $[\alpha]$ being constant for all s provides more wave artefacts. For example, the carbon cores of the 140 µm SiC fibres are clearer in many of the fibres in Fig. 3(b), while there are additional wave features around the blocks in Fig. 3(c) which decrease the quality of the image.

Figure 3
Reconstructed slices of the compound sample: left-hand column (a) and (c) $[\alpha_s = 0.0001]$ , right-hand column (b) and (d) $[\alpha_s = 0.0001s^2]$ ; with top row (a) and (b) S = 5 and bottom row (c) and (d) S = 51.

Another important question is the most appropriate value of S. From the arguments of the previous section one might suppose that increasing S indefinitely might give better and better results. However, increasing S can have important implications with regard to periodic structures. Compare Figs. 3(b) (S = 5) and (d) (S = 51) where there are new dark arcs near the C cores of the SiC fibres. In Figs. 4(a)–4(c) it can be seen how genuine features (thin curves) in the sinogram are suppressed and how the sinogram flattens as S increases. Note the horizontal lines in Figs. 4(a)–4(c) are due to intensity profile variations (see Titarenko et al., 2010c). Therefore, if the value of S is increased then the genuine features in the reconstructed image may be suppressed.

Figure 4
Reconstructed slices of the shark cartilage preprocessed by the proposed algorithm: a region of the sinogram (top row) and the central region of the slice (bottom row). Left-hand column (a) and (d) S = 5, central column (b) and (e) S = 20, right-hand column (c) and (f) S = 51; in all cases $[\alpha_s = 0.0001s^2]$ .

An interesting feature is the `spider' feature that is absent for S = 5, but becomes increasingly prominent as S increases (see Fig. 4). In order to consider its origin let us undertake a gedanken experiment. Consider imaging a simple ball for which we vary the distance between the ball and the centre of rotation, and for which there are no errors in the simulations (see Fig. 5). The proposed correction is applied to both sinograms using the same parameters (S = 100, $[\alpha]$ = 0.0001s²). Comparing Figs. 5(a) and 5(b) one can see that the footprint of the ball over the sinogram is affected more drastically when it is closer to the axis of rotation. It would seem that features which are closer to the centre of rotation are suppressed earlier when S increases and/or $[\alpha]$ decreases. Clearly this is unsatisfactory and leads to flattening of the detail very near to the centre of rotation, giving rise to the characteristic `spider' noted above [see Figs. 4(d)–4(f)].

Figure 5
Sinograms of a ball after the proposed method has been applied. The width of the sinogram is 2000, the radius of the ball is 20, and the distances of the ball from the centre of rotation are (a) 500 and (b) 50; the attenuation coefficient is the same for both cases; $[\alpha_s = 0.0001s^2]$ , S = 100.

In order to circumvent this problem it is preferable to increase S with the distance from the centre of the sinogram. For example, see Fig. 6 where S = 5 inside the dashed circle and S = 30 outside. Of course, one can also vary S continuously from the centre to the edges; however, even the stepwise approach gives good results. In the general case when varying S a compromise between suppression of the ring artefacts towards the periphery of the slice and the suppression of genuine features towards the centre of the slice should be found.

Figure 6
The central region of a slice of the shark cartilage reconstructed using the proposed algorithm: (a) S = 30 everywhere, (b) S = 5 inside the dashed circle, otherwise S = 30; $[\alpha_s = 0.0001 s^2]$ .

Now we provide some estimates for the times required to preprocess a sinogram with the method proposed in this paper. We have used a desktop computer with a dual core processor (E7200, Intel Core 2 Duo, 2.53 GHz) and Intel performance libraries. To form the n×n matrix H, see (21), it takes 0.11 s for n = 2048 (the case of the shark cartilage) and 0.44 s for n = 3880 (the case of the compound sample). Once the matrix is found the solution can be obtained in less than 0.01 s in both cases. Therefore, in order to process a three-dimensional volume it is better to find S n×n matrices first and store them to memory.

5. Conclusions

We have developed a new method for the suppression of ring artefacts suited to anisotropically attenuating objects which can be applied to both parallel and cone-beam geometries. This anisotropy may be because the object is laterally extended such as a plate, or because the attenuation across the sample varies heterogeneously. The series of projections can be processed as a set of parallel two-dimensional slices. The slices can be chosen by fixing either rows or columns of the CCD sensor. A solution can be quickly realised because an analytical solution exists that is a product of a known matrix and a vector. Critical to the method are appropriate choices of $[\alpha]$ and S. We have found that $[\alpha_s = \alpha_1 s^2]$ works well. As to the number of terms, S, we found that increasing S ends to preferentially suppress genuine features towards the centre of the slice, giving rise to a `spider' artefact. As a result a better solution is to increase S with distance from the centre of the sinogram.

Two datasets preprocessed by the proposed method were acquired at different synchrotron sources, i.e. using a parallel geometry of the beam, and reconstructed by FBP. Of course, the method can also be applied to data obtained by laboratory tomography, i.e. in cone geometry of the beam, if one preprocesses the data taken by a given row of the CCD sensor. The method is based on preprocessing two-dimensional data, i.e. data taken by one row of the sensor. This means that neighbouring slices/sinograms are preprocessed independently. However, in forthcoming papers the authors intend to extend the proposed method to three-dimensional data, i.e. use the assumption that two neighbouring slices should vary slightly. The authors believe that the studies presented here suggest that the method will be widely applicable to biological and medical specimens taken using white X-ray beams.

APPENDIX A

Formula derivation

Let there be a smoothing functional

$[M^\alpha[q] = \sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n-1} |p_{ij}-p_{i, j+1}|^2 + \alpha \sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n} |q_{ij}|^2, \eqno(24)]$

with the regularization parameter $[\alpha \,\,\gt\,\,0]$ . We know that an element q_ij can be written as a final sum

$[q_{ij} = \sum\limits_{s = 1}^S f_{is}c_{sj}, \eqno(25)]$

where vectors f_is are orthonormal, i.e.

$[\sum\limits_{i = 1}^m f_{is}\cdot f_{iv} = \left\{\matrix {1,& s = v,\cr 0,& s\ne v.}\right. \eqno(26)]$

Function (24) is continuously differentiable at any point c with elements c_sj, $[1\leq s\leq S]$ , $[1\leq j \leq n]$ . A point c can be a local minimizer of the smoothing functional only if

$[{\partial\over \partial c_{wu}} M^\alpha[q] = 0,\qquad w = 1,\ldots,S,\quad u = 1,\ldots,n, \eqno(27)]$

see, for example, Dennis & Schnabel (1996 ).

Taking (26) into account we obtain

$[\eqalignno{ {\partial\over \partial c_{wu}} \left(\sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n} |q_{ij}|^2\right)& = {\partial\over \partial c_{wu}}\left[\sum_{i = 1}^m\sum_{j = 1}^n\left(\sum_{s = 1}^Sf_{is}c_{sj}\right)^2\right]\cr & = \sum_{i = 1}^m2\left(\sum_{s = 1}^S f_{is}c_{su}\right)f_{iw}\cr & = 2\sum_{s = 1}^Sc_{su}\left(\sum_{i = 1}^mf_{is}\,f_{iw}\right) = 2c_{wu}. &(28)}]$

Since we have defined p_i0 = p_i1 and p_i,n+1 = p_in, then

$[\eqalignno{\sum\limits_{j = 1}^{n-1} |p_{ij}-p_{i, j+1}|^2 & = \sum\limits_{j = 0}^{n} |p_{ij}-p_{i, j+1}|^2\cr & = \sum\limits_{j = 0}^{n} |(q_{ij}-q_{i, j+1}) - (\tilde p_{ij}-\tilde p_{i,j+1})|^2. &(29)}]$

Therefore

$[\eqalignno{{\partial\over \partial c_{wu}}& \left(\sum\limits_{i = 1}^{m}\sum\limits_{j = 1}^{n-1} |p_{ij}-p_{i, j+1}|^2\right)\cr & = 2\sum\limits_{i = 1}^m\left\{\left[\sum\limits_{s = 1}^Sf_{is}(c_{sj}-c_{s,j+1})\right]- (\tilde p_{ij}-\tilde p_{i,j+1})\right\}\,f_{iw}\cr & \quad+ 2\sum\limits_{i = 1}^m\left\{\left[\sum\limits_{s = 1}^Sf_{is}(c_{sj}-c_{s,j-1})\right]-(\tilde p_{ij}-\tilde p_{i,j-1})\right\}\,f_{iw}\cr & = 2\sum\limits_{s = 1}^S(2c_{su}-c_{s,u-1}-c_{s,u+1})\sum\limits_{i = 1}^m f_{is}\, f_{iw}\cr & \quad - 2 \sum\limits_{i = 1}^m(2\tilde p_{iu}-\tilde p_{i,u-1}-\tilde p_{i, u+1})f_{iw}\cr & = 2(2c_{wu}-c_{w,u-1}-c_{w,u+1})\cr & \quad - 2 \sum\limits_{i = 1}^m(2\tilde p_{iu}-\tilde p_{i,u-1}-\tilde p_{i, u+1})f_{iw}. &(30)}]$

Using vectors c_w, g_w defined in (15) and (16) the local minimizer of the smoothing functional (24) should satisfy systems

$[Ac_w^T = g_w^T, \eqno(31)]$

where the matrix is the same for all w,

$[A = \left(\matrix{ 1+\alpha & -1 & 0 &\ldots &0&0\cr -1 & 2+\alpha & -1 &\ldots&0&0\cr 0 & -1 & 2+\alpha & \ldots &0&0\cr \vdots & \vdots &\vdots &\ddots&\vdots&\vdots\cr 0& 0& 0& \ldots & 2+\alpha & -1\cr 0& 0& 0& \ldots &-1 & 1+\alpha }\right).\eqno(32)]$

Note that the Hessian matrix of $[M^\alpha[q]]$ with respect to variables c_wu is a diagonal matrix with diagonal elements $[4+2\alpha\,\, \gt\,\,0]$ , i.e. it is positive definite. Therefore, the solution found is a local minimizer (see Dennis & Schnabel, 1996). Since the systems have only unique solutions, the element found, c, is the global minimum of Tikhonov's functional $[M^\alpha[q]]$ .

Acknowledgements

The authors would like to thank the Diamond Light Source for use of the X-ray light source. The use of the Advanced Photon Source was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. The piece of shark cartilage was provided by Professor Michael Coates (University of Chicago) and scanned by Mason Dean (Max Planck Institute). ST is grateful to EPSRC for funds through the `Collaborating for Success' grant. VT and AK are grateful to STFC for funding the project.

References

Antoine, C., Nygåard, P., Gregersen, Ø. W., Holmstad, R., Weitkamp, T. & Rau, C. (2002). Nucl. Instrum. Methods Phys. Res. A, 490, 392–402. CrossRef CAS Google Scholar
Bakushinsky, A. B. & Kokurin, M. Y. (2004). Iterative Methods for Approximate Solution of Inverse Problems, Vol. 577, Mathematics and its Applications. Dordrecht: Springer. Google Scholar
Banhart, J. (2008). Editor. Advanced Tomographic Methods in Materials Research and Engineering. No. 66, Monographs on the Physics and Chemistry of Materials. Oxford University Press. Google Scholar
Bardelli, L. et al. (2006). Nucl. Instrum. Methods Phys. Res. A, 569, 743–753. CrossRef CAS Google Scholar
Beerends, R. J., ter Morsche, H. G., van den Berg, J. C. & van de Vrie, E. M. (2003). Fourier and Laplace Transforms. Cambridge University Press. Google Scholar
Boin, M. & Haibel, A. (2006). Opt. Express, 14, 12071–12075. Web of Science CrossRef PubMed Google Scholar
Buzug, T. M. (2008). Computed Tomography: From Photon Statistics to Modern Cone-Beam CT. Berlin: Springer. Google Scholar
Davis, G. R. & Elliott, J. S. (1997). Nucl. Instrum. Methods Phys. Res. A, 394, 157–162. CrossRef CAS Web of Science Google Scholar
Dennis, J. E. & Schnabel, R. B. (1996). Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Philadelphia: SIAM. Google Scholar
Engl, H. W., Hanke, M. & Neubauer, A. (1996). Regularization of Inverse Problems, Vol. 375, Mathematics and its Applications. Dordrecht: Kluwer Academic Publishers. Google Scholar
Hsu, H. P. (1995). Schaum's Outline of Signals and Systems. New York: McGraw-Hill. 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
Kak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographics Imaging, Vol. 33, Classics in Applied Mathematics. Philadelphia: SIAM. Google Scholar
Krumm, M., Kasperl, S. & Franz, M. (2008). NDT E Int. 41, 242–251. CrossRef Google Scholar
Lavrent'ev, M. M., Zerkal, S. M. & Trofimov, O. E. (2001). Computer Modelling in Tomography and Ill-Posed Problems. Utrecht: VSP. Google Scholar
Mahajan, V. N. (2001). Optical Imaging and Aberrations, Part II, Wave Diffraction Optics. Bellingham: SPIE Press. Google Scholar
Moszyński, M. (2006). Radiation Detectors for Medical Applications, pp. 293–315. NATO Security through Science Series. Berlin: Springer. Google Scholar
Natterer, F. & Wübbeling, F. (2007). Mathematical Methods in Image Reconstruction, Vol. 5, SIAM Monographs on Mathematical Modeling and Computation. 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). MicroComputed Tomography: Methodology and Applications. Boca Raton: CRC Press. Google Scholar
Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V. & Yagola, A. G. (1995). Numerical Methods for the Solution of Ill-posed Problems, Vol. 328, Mathematics and its Applications. Dordrecht: Kluwer Academic Publishers. Google Scholar
Tikhonov, A. N., Leonov, A. S. & Yagola, A. G. (1998). Nonlinear Ill-Posed Problems, Vols. 1, 2 and 14, Applied Mathematics and Mathematical Computation. London: Chapman and Hall. Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. J. (2009). Appl. Phys. Lett. 95, 071113. Web of Science CrossRef Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. J. (2010a). J. Synchrotron Rad. 17, 540–549. Web of Science CrossRef IUCr Journals Google Scholar
Titarenko, S., Withers, P. J. & Yagola, A. (2010b). Appl. Math. Lett. 23, 1489–1495. Web of Science CrossRef Google Scholar
Titarenko, S. S. & Yagola, A. G. (2010). Moscow University Phys. Bull. 65, 65–67. Web of Science CrossRef Google Scholar
Titarenko, V., Titarenko, S., Withers, P. J., De Carlo, F. & Xiao, X. (2010c). J. Synchrotron Rad. 17, 689–699. Web of Science CrossRef CAS IUCr Journals 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
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. Piscataway: IEEE. 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 18| Part 3| May 2011| Pages 427-435

doi:10.1107/S0909049511006005

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Suppression of ring artefacts when tomographing anisotropically attenuating samples

1. Introduction

2. Ring artefact suppression

3. Practical applications

3.1. Test object 1: a high-aspect-ratio (laterally extended) sample

3.2. Test object 2: compound sample

4. Discussion

5. Conclusions

APPENDIX A

Formula derivation

Acknowledgements

References

research papers