Generalized Titarenko's algorithm for ring artefacts reduction

A fast algorithm for ring artefact reduction in high-resolution micro-tomography with synchrotron radiation is presented. The new method is a generalization of the one proposed by Titarenko and collaborators, with a complete sinogram restoration prior to reconstruction with classical algorithms. The generalized algorithm can be performed in linear time and is easy to implement. Compared with the original approach, with an explicit solution, this approach is fast through the use of the conjugate gradient method. Also, low/high-resolution sinograms can be restored using higher/lower-order derivatives of the projections. Using different order for the derivative is an advantage over the classical Titarenko's approach. Several numerical results using the proposed method are provided, supporting our claims.

1. Introduction

Micro-computed tomography essentially relies on the inversion of the Radon transform, which is a well known operator in inverse problems, especially in image reconstruction. Typically, a three-dimensional volume must be restored from a stacking of two-dimensional reconstructed images, where the ith feature image satisfies the following principle,

The most established method in tomography consists of converting X-rays into visible light with a scintillator and projecting onto a CCD using standard microscope optics. These projections can be used as input for computerized tomographic reconstruction. A common source of degradation in tomographic reconstruction is through the superimposition of so-called `ring artifacts'. There are several known reasons for the generation of ring artifacts in reconstructed images. Primarily, these are due to defective pixels on the CCD detector or non-linear response of individual detector elements. Defective pixels typically result in a sharp single-pixel artifact unless several consecutive pixels are flawed. Larger artifacts are generally a result of imperfect or dusty scintillator crystals, which can easily exceed a single pixel. Other causes have recently been reported, such as monochromator vibration (see Titarenko et al., 2010a) or thermal processes generating additional electrons in the CCD (see Titarenko et al., 2009). For 180° rotation of the sample, semi-circle artifacts are generated; only for 360° scans will the full ring artifact be generated. As such, these imperfections are independent of the projection angle and will appear uniformly for all angles in ring format, superimposed on the sample image. Regardless of the source of these artifacts, the resulting degradation of the image quality, obscuring relevant details, is detrimental to the computed tomography process.

In order to illustrate the artefacts of the reconstructed image using the analytical inversion with filtered backprojection algorithm, we consider the example of Fig. 1, where a human tooth was introduced to a standard tomography using the Medipix detector. From the corrupted sinograms we clearly see the ring artifacts in the reconstructed image. Some authors in the literature used to filter the reconstructed image transforming to polar coordinates in order to identify the ring structures (see Sijbers & Postnov, 2004, and references therein). Indeed, this is because the polar transformation maps circles to straight lines, which can be detected using an appropriate filter. Converting from polar coordinates to Cartesian and vice versa can be done with a computational cost of O(n²) where n is the size of the feature image. Unfortunately, changing to polar/Cartesian coordinates introduces additional noise in the reconstructed image. Polar transformation is a fast algorithm compared with the average computational cost of a reconstructed slice, which is O(n²logn) (see Andersson, 2005). The main disadvantage of polar/Cartesian algorithms is the need to filter the high frequencies with several thresholds on the cut-off frequency, depending on some prior knowledge of the sample under investigation. Another technique for removal of stripe artefacts in the sinogram is the Fourier–Wavelet approach, presented by Münch et al. (2009). It is a fast technique, with linear computational cost O(n) based on a discrete wavelet decomposition of the sinogram. This is a more sophisticated technique which is able to remove horizontal and vertical stripes in many imaging problems.

Figure 1 Reconstruction of the feature image f, using a corrupted sinogram . (a) f in Cartesian coordinates. (b) f in polar coordinates. Transforming the reconstructed image to polar coordinates, in order to remove vertical stripes, was a popular approach for removing ring artefacts. See Sijbers & Postnov (2004).

A fast tomographic experiment, as carried out at the Brazilian Synchrotron Light Laboratory, goes through the following scheme,

Step (B) is usually done with a fast reconstruction algorithm, specially designed to handle a huge tomographic data set. In our case we use only analytical reconstruction through `pyhst', a fast hierarchical inversion algorithm implemented for graphics processing unit (gpu) (see Mirone et al., 2013). Step (A) is the main contribution of this manuscript, where a fast filtering of each sinogram has to be carried out, prior to reconstruction.

Titarenko et al. (2009, 2010a,b,c) introduced a fast algorithm for filtering the sinogram prior to reconstruction, based on the optimality conditions of a quadratic functional. Their method is direct and with a linear computational cost, i.e. O(n). This has a great impact on building fast three-dimensional volumes for the tomographic data set. Therefore, even with more sophisticated inversion schemes, such as the expectation maximization algorithm (Helou et al., 2014; De Pierro, 1995) or recent algorithms based on compressed sensing techniques (Candes et al., 2006), the stacking process remains almost unaffected.

In this manuscript we present a generalization of Titarenko's algorithm, still preserving the computational linear cost. Our method preserves the sinogram structure also introducing some radial smoothness on each projection, a similar approach also presented in their work (Titarenko et al., 2010c). For real tomographic data, measured in the Brazilian Synchrotron Light Laboratory, our method is competitive with Titarenko's original approach.

The manuscript is organized as follows. §2 presents the theoretical fundamentation of Titarenko's original work further generalized in §3 and §4. Some numerical experiments are presented in §5, discussion in §6, and §7 gives the final conclusions.

2. Generalized variational approach

Theoretically, the pair {s,f} from equation (1) is related through the Radon transform (Deans, 1983), given below,

where = stands for the direction of the X-ray path, corresponds to a pixel of the feature image , and  = corresponds to a pixel of the sinogram . Typically, U and V are inner product spaces of square integrable functions, e.g. L² or a Schwarz space. It can be shown that the linear operator is invertible for an appropriate choice of U and V. A complete discussion about the action of can be found by Natterer & Wubbeling (2001).

This means that, if we are given s = , the function f could be recovered using f = with an appropriate formula for . In practice, lies in a discrete mesh of points , where N and R are assumed to be large, e.g. N ≃ 1000 and R ≃ 3000. Therefore, a good estimate of the solution f could be found using the classical filtered backprojection algorithm, which is fast (mainly due to fast Fourier transforms) and approximates the inverse operator (Gabor, 2009).

The approach of Titarenko et al. (2009, 2010a,b,c) searches for a regularized sinogram that does not contain horizontal stripes. These horizontal stripes are the main source for ring artefacts in the image reconstruction. A regularized sinogram is obtained through the minimization of an appropriate Tikhonov functional (Tikhonov 1963a,b; Yagola et al., 2002). A brief mathematical description of Titarenko's method is given in Appendix A. In this case, we are looking for a regularized sinogram, the solution of

This was also explored in the work of Twomey (1963) for the solution of integral equations of the first kind. The solution is presented in Appendix A. For completeness, the pseudo-code of the original Titarenko's ring suppression algorithm is given below,

Titarenko's Algorithm (TA):

Input: Corrupted sinogram

[1] Set and

[2] Set matrix as in (42)

[3] Define explicitly ; (see Titarenko et al. (2010a)

[4] Define

[5] Compute

Output: Restored sinogram

3. Generalized Titarenko's algorithm without angle dependency

Assuming that the deviation from s to m does not depend on the angle, i.e.

and changing to the discrete case, we are now searching for such that

where is a second-order finite difference operator given by

Since V is quadratic, solving = 0 guarantees a solution of (6). To make the calculations easier, let = + , with D related to the finite difference operator. Denoting as the ith column of , a straightforward calculation gives us

where is a pentadiagonal matrix,

Therefore, the optimality criteria gives us the noise vector in the same fashion as in (40) (see Appendix A),

with the average projection, as defined in (41).

3.1. Further generalizations

It is a well known fact, especially in discrete differential equations, that finite difference operators are closely related to tridiagonal or pentadiagonal matrices. The structure of the matrix depends on the type of approximation used for the associated derivatives. In the present case, we are also approximating derivatives and we have exactly the same problem. The finite difference operator used by Titarenko et al. (2010a,c) (see Appendix A) was easily expanded to a finite difference operator when minimizing (3). Of course, this can be expanded even more.

Considering only the discrete case, it is easy to realise that

is the general form of the finite difference operator applied to the ith column of matrix , and is the vector associated with the linear functional , i.e. is the jth row of a matrix . In the case of and we have the following matrices,

and

The number of rows of matrix strongly depends on the operator ; as a rule, the offset from the main diagonal indicates the number of rows of . As a consequence, one can easily verify that matrices and satisfy = and = .

Inspired by (3), (37) and (11), we want to solve the following generalized optimization problem,

Using the same notation as in (8), and letting , it is not difficult to show that

with the jth row of . From (18) and = , we must have

We refer to equation (19) as a generalization of Titarenko's algorithm, clearly including case (40) and (10). We can now choose a whole family of finite difference operators in order to build a ring suppression algorithm. The generalized approach is given below, where stands for the Kronecker product (Horn & Johnson, 1985).

 

Generalized Titarenko's Algorithm (GTA):

Input: Corrupted sinogram

[1] Set and

[2] Set finite difference vector/matrix

[3] Define matrix

[4] Define

[5] Solve

Output: Restored sinogram as in (4)

3.2. Finite difference coefficients

We remark that matrix acting on the operator is a convolution matrix based on the kernel h. In fact, from (41) we have h = and h = in (7). From Bengt (1988) and Mahesh (1998) we have a wide family of backward and forward finite difference coefficients, with different order of accuracy. Some coefficients are depicted in Table 1. The kernel vector is denoted h_kj, with k standing for the order of derivative and j for the corresponding accuracy.

Table 1
Finite difference coefficients to compute matrix ; see text for details

Derivative	Accuracy	Kernel
	1	h1,1 = (−1, 1)
	2	h1,2 = (−3/2, 2, −1/2)
	3	h1,3 = (−11/6, 3, −3/2, 1/3)
	6	h1,6 = (−49/20, 6, −15/2, 20/3, −15/4, 6/5, −1/6)
	1	h2,1 = (1, −2, 1)
	2	h2,2 = (2, −5, 4, −1)
	6	h2,6 = (469/90, −223/10, 879/20, −949/18, 41, −201/10, 1019/180, −7/10)
	1	h3,1 = (−1, 3, −3, 1)
	5	h3,5 = (−967/120, 638/15, −3929/40, 389/3, −2545/24, 268/5, −1849/120, 29/15)

3.3. Solving the optimality conditions (19)

We want to solve (12) in order to complete step [5] in the Generalized algorithm GTA. Taking as in (12), linear system (19) reduces to (40), which has analytical solution

where was explicitly given by Titarenko et al. (2010a). In the general case, a similar approach is not trivial. Huckle (1994) presents a Fourier-strategy for solving tridiagonal linear systems arising from elliptic partial differential equations. Although we are not restricted to partial differential equations, the matrix inversion presented in the original work of Titarenko et al. is highly connected to the work of Huckle.

It is easy to observe that matrix is rank deficient, i.e. = R-r, where r+1 is the length of the Kernel vector h, as presented in Table 1. For instance, in (40), = R-1 whereas = R-2 in (10). The relaxation parameter λ ensures that = R so that attains a minimum.

Let us denote = and = . We propose a fast solution of = through the use of the conjugate gradient method (CGM). This is possible since is a symmetric and positive definite matrix.¹ The CGM iterations (see Bazaraa et al., 2007) are particularly interesting for this linear system since is a convolution matrix. Indeed, the CGM strongly depends on matrix-vector product , for some , which is computed using fast Fourier transforms. In fact,

where h is the convolution kernel and is a symbol for the circular convolution.

3.4. Computational complexity

The solution of (19) for Titarenko's case, i.e. h = (1,-1) (see Table 1), could be computed really fast using a tridiagonal matrix algorithm (TDMA) (see Thomas, 1949). The computational complexity of TDMA is O(R). If the analytical inverse of is computed previously, we have to perform the matrix-vector product for each slice, which has complexity O(R²). This is not the best programming choice, and it is preferable to use TDMA per iteration. A generic pseudocode for the three-dimensional reconstruction, using Titarenko's method, is given by

For the generalized case, we still have complexity O(R) for step (a), using the conjugate gradient method. Therefore, the fast reconstruction strategy remains unnaltered. Operator stands for the analytical inversion with filtered backprojection, or an iterative strategy.

Fig. 2(a) depicts the CPU time for each slice, using algorithms TA and GTA in a high-level programming language like Python. The reconstructions for this particular example are presented in more details in §5. Even with an analytical formulation for the inverse matrix = in Titarenko's original approach, the CPU time is higher if compared with the solution of = through the CGM. Indeed, A is a symmetric sparse band matrix, and the iterations of CGM converge in at most R steps, where R is the number of rays in the sinogram. In fact, for a particular example with R = 2048, Fig. 2(b) shows the number of iterations for the convergence of the CGM using vector kernal h_1,2 and h_2,2 for the generalized Titarenko's algorithm. Even with TA computing just once the dense matrix , a matrix-product has to be performed, for each slice, which has a quadratic computational cost. Using CGM is a more efficient strategy since the matrix-product is computed through convolution.

Figure 2 Complexity of algorithms TA and GTA. (a) Comparison of the CPU time/slice used by Titarenko's algorithm and his generalization using the conjugate gradient method. (b) Number of iterations used by the conjugate gradient method with order 1 and 2.

3.5. Reduction by blocks

The generalized algorithm (19) could be easily implemented by blocks of columns. This is important whenever the blank scan (the flat field) is corrected in between projections. Also, such a strategy could be used to remove artefacts in the sinogram that are not eliminated using (19).

Let the sinogram be given by

where is such that N = cb. Each is restored according to the generalized algorithm, i.e.

where is the solution of (19), and is the average projection within the block, i.e. = . In this case, matrix = remains the same for each block.

3.6. Geometric mean sinogram

Solving (19) through the conjugate gradient method could be made fast using convolution and the fast Fourier transform. Therefore, obtaining two solutions is a simple task, say and . In this case, we have two approximate restored sinograms and , respectively. Each one is given by (20). If we take as the solution of (19) with kernel h_1,j, we are minimizing (40). On the other hand, if we take as the solution of the same linear system, with kernel h_2,j, we are minimizing (41).

Our aim is to combine these two sinograms in order to have a better one. A convex combination of and does not fit to our needs since we have to choose a proper weight to favour first or second derivatives. Instead, we set a new sinogram, denoted by , whose entries ( j,i) are defined by

for some previously defined tolerance > 0. If we take 0 we have the geometric mean.² The choice of ∊ is such that is close to the arithmetic mean.

4. Generalized Titarenko's algorithm with angle dependency

This case was already presented by Titarenko et al. (2011), where the corrupted sinogram also depends on the angle of rotation. In this case, hypothesis (5) is replaced by

taking into account physical properties that depend on the angle of rotation θ. After a similar discussion given in Appendix A, a functional is minimized and a linear equation is obtained. Without going into further detail, we claim that a similar analysis of §3 is easily extended for this case, using a finite difference operator . Here, we are looking for a correction matrix such that = is the updated sinogram, with the measured one.

In this case, the rows of the matrix are written as a linear combination of a previously known basis, i.e. with the unknown coefficients and the basis. After computing the optimality conditions, the following linear system is obtained,

where are the kth columns of matrices and , respectively. Here, matrix is given by = with = . In a matricial form, the correction matrix is the one that satisfies = . Generalization comes from the fact that we can choose a kernel vector h, see Table 1, in order to define the finite difference matrix . Taking h = , matrix becomes exactly the one proposed by Titarenko et al. (2010a), with an analytic inverse. Further details of this algorithm can be found in the work of Titarenko et al. (2011).

 

Generalized Titarenko's algorithm with angle dependency [GTA(θ)]

Input: Corrupted sinogram ,

[1] Define orthonormal by rows

[2] Set =

[3] Define matrix =

[4] Define matrix =

[5] Solve = , for eack k =

Output: Restored sinogram =

 

Compared with the GTA algorithm, GTA(θ) demands a high computational cost. Loop [5] could be eliminated if the inverse of matrix is computed a priori, although this is only true for kernel vector h = (1, −1); in this case, = . For different kernels, the explicit inverse matrix is no longer available, and the conjugate gradient method is used.

5. Numerical experiments

In 2011, work began in earnest on the design of an imaging beamline at the Brazilian Synchrotron Light Source (https://lnls.br ). IMX, which has an electron source size of 391 µm × 97 µm and beam divergence of 808 µrad × 26 µrad, was built on a 2 T bending magnet of the 1.37 GeV storage. This beamline can operate in either white beam or monochromatic beam. The white-beam energy spectrum ranges from 4 keV to 15 keV, with a photon flux at the sample position of approximate 10¹⁵ photons s⁻¹. Several scintillators are available at the beamline, the most utilized being the YAG:Ce (scintillator crystal), which has an emission wavelength of 550 nm and a photon output per 8 keV. The scintillator is mounted on a commercially available microscope system from Optique Peter³ (Douissard et al., 2012) coupled with a PCO.2000 CCD camera (Germany). The PCO.2000 is a high-resolution 14-bit CCD cooled camera with a quantum efficiency of 55% at 500 nm and has a full well capacity of 40000. The resulting array size of this detector system is 2048 × 2048, with pixel size from 0.37 µm to 3.7 µm. For a typical tomographic scan, 1000 images were acquired over 180°. A schematic of the experimental system for the acquisition of projected data is shown in Fig. 3.

Figure 3 Schematic of the imaging beamline IMX at LNLS.

5.1. Medipix detector

Hybrid silicon photon-counting detectors can meet the requirements of very high dynamic range and high detection efficiency specially needed for tomographic applications using synchrotron radiation. Compared with other families of hybrid pixel detectors, such as Pilatus (see Kraft et al., 2009) and XPAD (see Pangaud et al., 2008), the Medipix3RX readout chip enables a higher spatial resolution with a pixel size of 55 µm by 55 µm. Furthermore, this chip takes advantage of the advanced CMOS (complementary metal oxide semiconductor) technology to allow a high level of functionality in each pixel, enabling reconstruction of the charge generated by a photon detected on several neighbour pixels (charge summing mode) or allowing up to eight energy thresholds per cluster of four pixels (spectroscopic mode). Unfortunately, some characteristics of analog and digital circuits of the ASICS (application specific integrated circuit) are responsible for a non-linear relationship between measurement and incident flux. One of these characteristics is the shaping time directly related to the dead time of the system, which is responsible for pile-up effects (see Rinkel et al., 2011). The pixel-to-pixel dispersion of these non-linear parameters cannot be corrected by simple flat-field normalization. In tomography, these dispersions are responsible for ring artefacts, in addition to global false estimation of the reconstructed attenuation coefficients.

To assess the robustness of the proposed algorithms for ring artefact correction, tomographic measurements were performed on the tomographic beamline of the Brazilian Synchrotron (LNLS) with a Medipix complementary detector. A Medipix3RX ASIC bump bonded to a 300 µm silicon sensor was used in conventional single-pixel mode configured with single 24-bit counters and high-gain mode (see Gimenez et al., 2011). The detector was read out using Medipix3 USB interfaces and Pixelman software (see Vykydal et al., 2006).

5.2. Examples

In the absence of noise, the theoretical sinogram satisfies the following mathematical relation (Helgason, 1980),

We refer to constant z as the mass of the sinogram, i.e. the mean remains constant for all possible values of θ. For noisy data, (28) is no longer true. For the LNLS synchrotron data, as the beam current decays in time, function is far from being constant, with a monotonic behaviour.

Since in (28) is not constant for real measured data , the standard deviation (std) of each projection is a function , say

Function w = is shown in Fig. 4 for some measured sinogram at the angle set . We take the regularization parameter λ [see (14)] as the standard deviation of , i.e.

As is clear from (30), the value of λ depends on the measured sinogram . As = is closer to a constant vector, λ is lower as we are dealing with a well behaved sinogram (at least in theory). The opposite remains true, i.e. the more distant is from a constant vector then the higher will be the penalty since we have an ill-behaved sinogram.

Figure 4 Function w = , as defined in equation (29), is the standard deviation of each measured projection, for a fixed sinogram. In this particular example the sinogram is presented in Fig. 5(a).

We use the following notation:

is the restored sinogram = + ; where is the solution of (19) with kernel h_k,j given in Table 1;

= is the geometric mean sinogram, with parameter λ defined in (30), using and as references [see (25)];

is the measured sinogram, usually corrupted as in (5).

A measured sinogram at the LNLS, , is shown in Fig. 5(a). The region within the rectangle is presented in Fig. 5(b), where several horizontal stripes are visible. Each stripe generates a circle in the reconstruction technique.

Figure 5 Sinogram I (see text for details). (a) Original sinogram from the LNLS beamline. (b) Zoomed region from the original sinogram. (c) Sinogram . (d) Smoothed sinogram . (e) Geometric mean sinogram .

The generalized Titarenko's suppression algorithm, with kernel h_1,3, see Table 1, gives the result in Fig. 5(c). Using a high-order finite difference, the horizontal stripes are clearly reduced. However, as can be noted in Fig. 5(c), there are two severe horizontal stripes remaining in the restored sinogram. The one on the top is due to the average projection [see vector in (41)]. The second is probably due to dead/damaged pixels in the CCD camera.

Fig. 5(d) is the restored sinogram using the generalized algorithm with kernel h_2,2. The artefact corresponding to the first stripe is strongly reduced, although the second stripe still remains, as it is a major corruption in the data. Fig. 5(e) depicts the geometric sinogram using the combination of Figs. 5(d) and 5(c).

The projection at the angle = for the restored sinogram is presented in Fig. 6(a), compared with the same projection of the geometric sinogram . As expected, the geometric map overcomes the sinogram , as shown in the rectangle, zoomed in Fig. 6(b).

Figure 6 Projections at = . (a) Projection versus projection . (b) Zoomed region of part (a).

To remove the second horizontal stripe (see Fig. 5e), we adopt the block strategy presented in (24), using only two blocks. The restored sinograms [only the critical region of Fig. 5(a)], are shown in Figs. 7(a) and 7(b). The first was obtained using kernel h_2,2 at each block, whereas the second was obtained using the geometric mean sinogram of each block, taking h_2,2 and h_1,3 as appropriate convolution kernels.

Figure 7 (a) Restoration using two blocks [see equatiion (24)] with kernel h2,2 at each block. (b) Geometric mean sinogram using two blocks. (c) Reconstruction with the sinogram of (b). (d) Magnified part of (c).

The first horizontal stripe present in the restored sinogram (see Fig. 5c) leads to the circular shadow outside the sample. The second horizontal stripe, uncorrected by first-order finite differences, leads to only one strong ring artefact in the feature image. Using second-order finite differences, i.e. using , we obtain the reconstruction in Fig. 8(c). Since the first horizontal stripe was reduced, the shadow outside the sample, as in Fig. 8(b), was dramatically reduced. Nevertheless, since the second stripe still remains in the restored sinogram, the reconstructed feature image presents one strong ring artefact. Through the geometric mean sinogram (see Fig. 5e) we obtain the reconstruction in Fig. 8(d).

Figure 8 FOV of reconstruction from sinograms presented in Fig. 5. See text for details.

To reduce the leading ring artefact, still present in the reconstructed images [see Figs. 8(b), 8(c) and 8(d)], we use the filtered backprojection on the sinogram of Fig. 7(b), which is the geometric combination of (by blocks) and (by blocks). As depicted in Fig. 7(c), the leading ring artefact was completely eliminated, while still preserving the aspects of the image; as seen in Fig. 7(d) [zoomed part of Fig. 7(c)].

In order to illustrate the action of the derivative kernel h_j,k, see Table 1, we present a low-resolution simulated example. A 327×327 image, representing a simulated micro-gear, is shown in Fig. 9(a). After adding noise to the 527×180 sinogram (i.e. horizontal stripes), the filtered backprojection is applied to the corrupted sinogram, originating the reconstructed Fig. 9(b), where the ring artefacts are clearly visible.

Figure 9 Simulated example. (a) Image with resolution 327×327, representing an ideal micro-gear. (b) Reconstructed image (with strong ring artefacts) after adding horizontal stripes to the 527×180 sinogram.

The generalized Titarenko's algorithm (see §2) with different orders for the derivative was applied to the corrupted simulated gear data. The results are depicted in Fig. 10. Part (a) presents the result obtained using the original Titarenko's algorithm (see §2 for an algorithmic description of TA). In this case the strong ring artefacts were reduced, but new artefacts are introduced, mainly because the sinogram has low resolution in the pixel axis. We remark that Titarenko's original algorithm is obtained with the generalized algorithm and kernel vector h_1,1. Using the second-order derivative, we reduce the stripe effects through smoothness in the radial direction of the sinogram. The reconstructed image using second derivatives, i.e. using kernel vector h_2,1 and h_2,2, are presented in part (b) and (c), respectively. Although the rings still remain in the reconstructed images, they were reduced if compared with (a). Using the block strategy and different kernel vectors, we obtain Figs. 10(d), 10(e) and 10(f).

Figure 10 Simulated example, for TA and GTA. (a) TA for corrupted data. (b) GTA with h1,2. (c) GTA with h2,2. (d) GTA with h1,2 and 6 blocks. (e) GTA h2,2 and 60 blocks. (f) GTA with h3,1 and 6 blocks. See text for details.

The Wavelet-FFT approach for ring removal, presented by Münch et al. (2009), is based on a discrete wavelet decomposition of the sinogram. Without going into further details, the algorithm to recover the sinogram is given by

where is the measured sinogram, L is the highest decomposition level using the wavelet type ω, and σ is a damping factor, fixed for each decomposition level. For example, ω could represent a signature to Daubechies wavelets (e.g. DB25 or DB42) or to a Haar wavelet (DB1). Fig. 11 presents some reconstructions using DB25 with three different choices for . In comparison with the generalized Titarenko's algorithm, we note that that the method is highly sensitive to the choice of .

Figure 11 Reconstruction of the simulated micro-gear using the Wavelet-FFT correction and Daubechies wavelet DB25 (Münch et al., 2009). (a) Parameters L = 7, = 2. (b) Parameters L = 3, = 3. (c) Parameters L = 2, = 7.

The block strategy was applied to the corrupted gear data, originating the reconstructions in Figs. 10(d), 10(e) and 10(f). The restored sinograms were recovered using the generalized Titarenko's approach using b = 6 blocks [see (23)], i.e. using GTA in between 30 angles.

The sinogram of the rock, presented in Fig. 5, was restored using the generalized Titarenko's algorithm with angle dependency [see §4]. Using kernel vector h = (1,-2,1) from Table 1 we obtain the results in Fig. 12. The restored sinogram presented in Fig. 12(b) is similar to the one obtained using the classical Titarenko's approach. The same strong stripe artefact remains, even with a smoothing kernel. In this example, GTA(θ) is competitive with GTA, but with a higher computational cost. Fig. 12(a) shows the correction matrix , which is non-constant through columns. Fig. 13 depicts GTA(θ) in the simulated micro-gear example, with different values of parameter S, and with fixed parameter λ, as in (30).

Figure 12 Sinogram restoration using GTA(θ) with kernel h = (1,-2,1) and S = 5. (a) Correction matrix . (b) Updated sinogram . See §4 for details.

Figure 13 Sinogram restoration using GTA(θ) with kernel h = (1,-2,1) and different values of . (a), (b), (c) S = 2, S = 5 and S = 20, respectively. (d), (e), (f) Reconstruction of sinograms (a), (b) and (c), respectively.

Finally, some image reconstructions, with corrupted sinograms obtained with the Medipix detector, are shown in Figs. 14 and 15. As Medipix are low-resolution images, the same discussion of the micro-gear example of Fig. 9 applies here. Noisy and low-resolution sinograms are better restored with a higher-order derivative, and using a block strategy.

Figure 14 Reconstruction of the human tooth of Fig. 1 with the Medipix detector. (a) Using sinogram . (b) Using the geometric sinogram. (c) Enlarged image of (a). (d) Enlarged image of (b).

Figure 15 Reconstruction of a rock sample with the Medipix detector. (a) Using the corrupted sinogram. (b) Using the geometric sinogram. (c) Enlarged image of (a). (d) Enlarged image of (b).

6. Discussion

As mentioned in the Introduction, the generation of ring artifacts in reconstructed images are usually related to non-linearity of the pixels' response. The model of sinogram corruption given by (5) accounts for the presence of strips. It assumes that the deviation from s to m does not depend on the angle. In other words, this deviation is supposed to be independent of the sample inhomogeneity.

Consider a pixel receiving a signal (intensity or number of photons) denoted by I₀ without sample and I in the presence of the sample. The associated measurements are denoted Im₀ and Im, respectively. The measured sinogram m and the original one s of (5) are obtained by calculating the measured and real X-ray attenuations of the object, respectively,

We can deduce the value of the deviation term,

The second term of the last equation, -log(Im/I), depends on the object composition and in this way on the projection angle for each pixel. Most non-linear effects, such as the pile-up phenomenon, tend to increase with the incident flux. This term is expected to become smaller with higher flux, resulting in better artefact correction for high-attenuation objects. This interpretation could explain the behaviour of the ring suppression algorithm on the images obtained with Medipix (see Figs. 14, 15 and 16). The ring suppression within the samples is really efficient for the tooth (Fig. 14) and for the rock (Fig. 15). This could be related to their high attenuations with maximum attenuation values over the whole sinograms of 8.3 and 6.2, respectively, for these samples. For the last sample, the ring suppression within the sample is only partial, which could be explained by this low attenuation (maximum attenuation value of 1.8 over the whole sinogram).

Figure 16 Reconstruction of the human tooth of Fig. 1 with the Medipix detector, using the Wavelet-FFT correction and Daubechies wavelet DB25 (Münch et al., 2009). (a) Parameters L = 7, = 2. (b) Parameters L = 3, = 3. (c) Parameters L = 2, = 7.

We are following the approach of Prince & Willsky (1990), where a quadratic functional is minimized. This was also the approach adopted by the classical Titarenko's approach. In fact, using a quadratic norm, we obtain an easy linear system to solve. This is one of the main goals of Titarenko's algorithm. Using a total variation norm (TV) penalty should be more appropriate, but this would introduce another type of algorithm, more difficult and not as fast as the one proposed. One of the most important aspects of imaging beamlines at a synchrotron is computational time, where GTA is appropriate.

7. Conclusions

We have presented a fast algorithm for ring artefact reduction in tomography by generalizing the approach proposed by Titarenko and collaborators. We provided several numerical and simulated results using the proposed method. The algorithm was tested on experimental data acquired with two different detectors: a scintillator mounted on a microscope system coupled with a CCD camera, and a counting detector based on a 300 µm Si sensor connected to a Medipix3RX readout ASIC. For a Medipix detector (with low resolution), we obtain good results using the generalized approach with high accuracy on the finite difference operator.

The main limitation of the method relies on the precision of the model of data corruption. Further work is focused on Titarenko's algorithm allowing the dependency of the correction term with the projection angle, which has already been carried out in their work of 2011 (Titarenko et al., 2011), i.e. by considering = instead of (5). We have applied such a technique in this manuscript, but the computational cost is high. We have also compared our strategy, without the angle dependency, with the Wavelet-FFT algorithm presented by Münch et al. (2009). The wavelet approach is fast, with a linear computational cost O(N). Although it is a very competitive method, it depends on three parameters, , making it difficult to choose the best configuration. In contrast with the present method, which depends on only one parameter, this is clearly a disadvantage of the wavelet algorithm.

Applying a TV penalty to the corrupted sinogram could be a good choice to restore the sinogram, penalizing smaller perturbation. This is an excellent starting point for a new study on rings artefacts reduction.

APPENDIX A

Titarenko's Algorithm

In real experiments, the measured sinogram is a corrupt version of the original one, that we denote by , i.e.

where is a sinogram to be determined, obeying some mathematical criteria (see Titarenko et al., 2011). Following the classical approach of Rudin et al. (1992), the restored sinogram is the one that minimizes the TV, subject to constraints of mean, say , and standard deviation, given by ,

Here, the total variation norm is defined by TV(s) = . The equivalence of the constrained and unconstrained optimization problems in (36) was previously shown by Chambolle (2004). The choice of the regularization parameter λ depends on the nature of the problem and could be estimated with an appropriate method (see Yagola et al., 2002). A TV approach will penalize smaller perturbations on the corrupted sinogram, but in this case the optimality conditions are more difficult to solve.

The approach adopted by Titarenko et al. (2010a,b,c) searches for a solution that minimizes the quadratic functional

This approach was also explored in the work of Prince & Willsky (1990, 2002), where the sinogram is restored through the mimization of the above functional, subject to the Ludwig–Helgason consistency conditions (Ludwig, 1966; Helgason, 1980).

Assuming that the deviation from s to m does not depend on the angle, i.e.

we arrive at the minimization, in terms of the radial (or noise) function n = n(t), i.e.

Assumption (38) is motivated by the stripes along the θ-axis on each sinogram.

In a discrete sense, the sinogram m is represented by matrix , and we are looking for a noise vector such that V = is minimum, with the integral replaced by the sum and the derivative replaced by finite differences. Since we are minimizing a quadratic, the function attains a minimum where the gradient is zero, i.e. = 0. The first-order conditions and a tridiagonal system must be solved to find . Titarenko's algorithm now follows,

with being a constant tridiagonal matrix, a first-order finite difference and the average of the projections, i.e.

Titarenko et al. (2010a) presented the exact analytical inverse of using properties of hyperbolic trigonometric functions. For completeness, matrix is a tridiagonal matrix, whose diagonal entries are given by

APPENDIX B

Octave codes

The codes for ring suppression, as presented in this manuscript, are presented below, in Octave open-source programming language. Function ring() performs the correction by columns, while function ring_blocks() corrects by blocks of columns; see (20), (23) and (24). Auxiliary function performs the matrix-vector product (22) using the convolution function conv() implemented by Octave. Finally, function ringCGM() iterates using the Conjugate gradient method in order to solve linear system (21). Functions kron() (Kronecker product), ones() and std() are also implemented in Octave (see Eaton et al., 2009).

By changing vector h at line 6 of functions ring() and ring_blocks(), we modify the order of the derivative, according to Table 1. Transcription of these codes to a high-level language such as Python could be easily done. We emphasize that the original Titarenko's algorithm is obtained using h = [1-1] at line 6 of function ring(). The same applies to function ring_blocks().

Function ring_a() stands for the algorithm GTA(θ) (see §4) with S as the input parameter; see Titarenko et al. (2011) for further details on matrices and value .

 

function  new = ring(old);

 

alpha = 2*std(std(old));

N = size(old,2);

h = [-1   2  -1];

n = ringCGM(h,alpha,f);

new = old + kron(n,ones(1,N));

 

 

function  new = ring_blocks(old,step);

 

alpha = 2*std(std(old));

R = size(old,1);

N = size(old,2);

h = [-1   2  -1];

blocks = fix(N/step);

 

for  k = 1:blocks,

   col = (k-1)*step+1:k*step;

   sino_block = old(:,col);

   q = ringCGM(h,alpha,f);

   new(:,col) = sino_block + diag(q)*ones(R,step);

end

 

 

 

s = conv(x,fliplr(h));

u = s(length(h):length(x));

y = conv(u,h);

 

 

function  x = ringCGM(h,alpha,f)

 

x0 = zeros(size(f));

w = -r;

x = x0 + a*w;

B = 0;

 

r = r - a*z;

      break;

   endif

   w = -r + B*w;

   x = x + a*w;

end

 

 

function   new = ring_a(old,S);

 

R = size(old,1);

N = size(old,2);

F = zeros(S,N);

i = 1:N;

 

for  k = 1:fix(S/2),

   D(2*k,:) = sqrt(2/N)*cos(pi*i*2*k/N);

   D(2*k-1,:) = sqrt(2/N)*sin(pi*i*2*k/N);

end

D(1,:) = sqrt(1/N) * ones(1,N);

h = [-1   2 -1];

 

 

alpha = 2*std(std(old));

 

for  w = 1:N,

   new(:,w) = old(:,w)+ringCGM(h,alpha,G(:,w));

end

Footnotes

¹Indeed, for all we have = > 0.

²This is always lower than or equal to the arithmetic mean.

³https://www.optiquepeter.com/en/applications-synchrotron.php .