Ring artifacts correction in compressed sensing tomographic reconstruction

02/05/2015 ∙ by Pierre Paleo, et al. ∙ 0

We present a novel approach to handle ring artifacts correction in compressed sensing tomographic reconstruction. The correction is part of the reconstruction process, which differs from classical sinogram pre-processing and image post-processing techniques. The principle of compressed sensing tomographic reconstruction is presented. Then, we show that the ring artifacts correction can be integrated in the reconstruction problem formalism. We provide numerical results for both simulated and real data. This technique is included in the PyHST2 code which is used at the European Synchrotron Radiation Facility for tomographic reconstruction.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 16

page 17

page 18

page 20

page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

1.1 Rings artifacts in tomographic reconstruction

During a tomographic acquisition process, some flaws in the experimental setup can lead to unwanted artifacts. One instance is the presence of defective pixels in the detector, which appear as lines in the sinogram when the defects are independent of the projection angle. These spurious lines give raise to rings artifacts in the reconstructed object. Even after preprocessing steps like flat-field correction and median filtering, these artifacts can remain and are detrimental to the reconstruction quality. Therefore, multiple techniques have been developed to tackle this problem.

1.2 Related work

Various techniques have been proposed in the literature to reduce or suppress the rings artifacts. As reported in [ringscomparison]

, these techniques can be classified into two groups : sinogram preprocessing and reconstructed images post-processing. The preprocessing methods aims at detecting and correcting the spurious lines in the sinogram before applying the reconstruction process, thus, rings do not form if the method succeeds. On the other hand, post-processing techniques work directly on the reconstructed image, trying to extract the concentric circles and filter them. These methods often perform a transformation into polar coordinates to transform the concentric circles into straight lines

[prell09].

A comprehensive comparison of ring artifact removal methods can be found in [ringscomparison]. Although these methods certainly provide satisfactory results in their limited framework, the authors report that no existing method is really suitable for correcting different types of rings, since they always introduce other distortions. A recent work [miqueles] reports a compressed sensing approach for rings artifacts reduction using a total variation denoising of the sinogram before calling the reconstruction routine. It is a generalization of Titarenko’s Algorithm [titarenkoReg] which consists in a regularization of the sinogram. This can also be classified in the sinogram pre-processing techniques.

2 Preamble


In this section, we introduce the principle and the formalism of compressed sensing tomographic reconstruction. This formalism is extended in section 3 for ring artifacts correction.

2.1 Compressed sensing tomographic reconstruction

Computed tomography aims at reconstructing an image from a set of projections . Here denotes the acquired sinogram, is the slice to be reconstructed and is the projection operator (while its adjoint is the back-projection operator). The classical filtered backprojection algorithm enables to reconstruct the image, but the number of projections should be of the same order of the number of rows in the image to have an acceptable reconstruction according to the Shannon-Nyquist sampling theorem. This is often impracticable, and the subsampling leads to artifacts in the reconstructed image.

Compressed sensing techniques exploit a-priori knowledge on the image, like its sparsity, in order to bypass this limitation. Instead of computing a closed form solution like in the filtered backprojection technique, tomographic reconstruction by compressed sensing amounts to an optimization problem

(1)

where is a fidelity term of with respect to the acquired data (henceforth is denoted ), and contains a-priori knowledge on the image. In general, the regularization term makes the problem non-smooth, which precludes from using usual Gradient-like algorithms .

Advances in convex analysis provide adapted methods, based on proximal splitting methods [proxsplitting], which are a generalization of projected gradient. One instance is the Iterative Shrinkage-Thresholding Scheme (ISTA). For a functional split into a smooth term and a non-smooth term , the first order condition at an optimum reads

which suggests the use of a fixed point iterative scheme. The operator is called proximal operator :

(2)

so that one step of ISTA reads

where is the Lipschitz constant of the gradient of the fidelity term. Usually, the data fidelity term is a L2 norm, so the calculation of the proximity operator of is straightforward. On the other hand, the proximity operator of the regularization term does not always have a closed form expression.

2.2 Total variation regularization

Depending on the regularization term used in (1), the reconstruction can yield very different results. If the regularization term is null, the problems amounts to a least-squares reconstruction. This approach tends to blur the edges in the reconstructed slice. A commonly used regularization term for image denoising is the Total Variation which was introduced in [ROF] as the Rudin-Osher-Fatemi model. This prior has the property to preserve the image edges, which is essential especially in tomographic reconstruction where the object in the sample should be distinguishable. In a discrete framework, the Total Variation is the L1 norm of the gradient :

(3)

Given an observed sinogram , the total variation tomography reconstruction problem aims at finding the regularized image satisfying

(4)

where weights the regularization with respect to the data fidelity term, and is the projection matrix. If

is equal to the identity matrix, the problem (

4) amounts to the so-called total variation denoising problem :

(5)

In [michel11], a dual approach enables to compute the proximity operator of the total variation :

(6)

Knowing the proximal map of the Total Variation regularization term, the denoising problem can be solved by an iterative shrinkage-thresholding scheme. In PyHST2, we use an accelerated version known as the Nesterov algorithm [nesterov] or FISTA [beckteboulleIEEE].

However, tomographic reconstruction is not a simple denoising problem, since a projection matrix appears in the problem (4). Constructing a smooth dual of (4) would entail to invert which is an ill-posed problem. The total variation tomographic reconstruction is tackled with two nested FISTA loops, each iteration being the solution of a denoising subproblem [beckteboulleIEEE]

(7)

Here is the Lipschitz constant of which is calculated with the power method.

2.3 Dictionary learning

Total Variation regularization performs well for piecewise-constant images since edges and uniform regions are reinforced. However, for non piecewise-constant images, the cartoon effect might be prejudicial for the reconstruction quality. Thus, another regularization technique has to be considered for such images.

Most natural images have an intrinsic sparsity which can be recovered by an adapted transform or by building the best sparsifying basis. The latter technique is called dictionary learning. Given a set of acquired signals , a number of

basis vectors (or atoms)

are built. Each signal is expressed as a linear combination of the atoms. Dictionary learning is a joint optimization of the atoms and the coefficients under the constraint of sparsity :

(8)

Here denotes the zero norm counting the number of non-zero component of a vector. The dictionary is learned with K-SVD [ksvd].

In image processing problems like tomographic reconstruction, the image is divided into square patches of size pixels. Every patch area of index is expressed as a linear combination of the atoms :

(9)

That is, for pixel of the image :

(10)

where denotes the patch containing the pixel , and is the center of this patch so that belongs to the patch support.

To overcome discontinuities in patches borders, the PyHST2 [pyhst2] approach allows the patches to overlap. Following formalism used in PyHST2, denotes the indicator function of patch while denotes the indicator function of the patch center111each patch center is a subset of the patch such that the patch centers are not overlapping. Every image pixel is then a linear combination of possibly overlapping atoms :

(11)

The tomographic reconstruction problem is

(12)

where is the matrix containing the patches coefficients, and

(13)

is a term promoting patches overlapping. Again, the optimization problem (12) is solved with the FISTA algorithm.

3 Rings correction in compressed sensing reconstruction

We now present how rings correction can be handled directly in the reconstruction process by integrating additional variables in the functional to minimize. This approach is independent of the regularization used and can therefore be applied in various frameworks like total variation and dictionary learning. Sinogram pre-processing techniques modify the acquired data to filter the unwanted lines. This filtering often introduces new artifacts. On the other hand, image correction techniques can also add new artifacts when circular features are detected as artifacts ; and the forward and backward Cartesian-polar coordinate transforms lead to a loss of precision even with a bilinear interpolation. When the rings correction is done in the reconstruction process, the data is not modified ; it is an initial guess of the solution for the iterative method. The solution minimizes the euclidean distance from the reference slice – the backprojected sinogram – under regularization constraints. The rings correction is included in both the constraints and the functional to minimize.

In our approach, the rings correction consists in splitting the sinogram into two components : the “genuine” sinogram and the spurious straight lines giving rings after back-projection. This splitting is done in the reconstruction process, so the two components are updated after each iteration. In the end, only the valid sinogram component is kept while the rings variables are discarded. We give two examples of frameworks using this approach : total variation regularization and dictionary learning reconstruction.

3.1 Rings correction in total variation framework

Rings correction can be handled by additional variables in the fidelity term . These artifacts come from a spurious line along projection angles in the sinogram. Thus, rings variables are stacked in a vector and added to the sinogram for each projection. The fidelity term for one projection reads

(14)

here and do not have the same dimensionality ; means that a vector of rings variables is added to each line of the sinogram as illustrated on Figure 1. The functional is

(15)

being a parameter weighting the relative importance of spatial regularization, and being a penalization parameter for the rings.

While sinogram preprocessing techniques filter the lines parallel to the projection angle, this approach forces the sinogram to be decomposed as a sinogram and rings variables . The sparsity constraint forces the rings variables to have only a few not null components, since the L1 norm is a good approximation of the sparsity-inducing L0 norm.

We emphasize the fact that the sinogram decomposition into a genuine sinogram and spurious rings is not a pre-processing technique ; the rings removal is intrinsically part of the reconstruction process. At each iteration, the image and the rings variables are adapted to minimize the energy ).

Figure 1: Principle of the rings separation. is the projection angle and is the detector bin index. The vertical orange and green lines represent spurious lines giving raise to ring artifacts. The decomposition forces the ring values to be captured in the vector (independent of the projection angle). In the end, only the part without the rings is returned.

The minimum of is found with the total variation regularization solver presented in 2.2. This iterative algorithm has a step size where is the Lipschitz constant of the gradient

. This constant is an upper bound of the largest eigenvalue of the Hessian

. Since the gradient is now taken with respect to both image and rings variables, the Hessian is

(16)

its largest eigenvalue is calculated with the power method.

Once the Lipschitz constant is obtained, one iteration of the FISTA denoising algorithm needs to compute

(17)

Here denotes the augmented vector containing both image and rings variables, (resp. ) is the part of vector containing image (resp. rings) variables. Since the squared L2 norm is separable, the proximal operator (17) can be written

(18)

In short, the proximal operator is separable with respect to the image and rings variables. This is convenient because rings and image variables can be updated by solving two separate subproblems in a FISTA iteration.

The first term in (18) is the denoising problem (5). The second term is the proximity operator of the L1 norm, which has a closed-form expression called soft thresholding operator :

(19)

3.2 Rings correction in Dictionary Learning framework

Similarly to total variation (15), the proximal operator is separable with respect to patch variables and rings variables . The problem (12) becomes

(20)

On the other hand, in this case, the non-smooth term is only the L1 norm which has a closed-form proximal operator (19). Thus, the denoising problem is straightforward, and only one FISTA loop is required :

(21)

3.3 Preconditioning the optimization problem

In both cases of Total Variation (15) and Dictionary Learning (3.2), the optimization process can be sped up by a precondition process. It is well-known that in the back-projected slices, the low spatial frequencies are overrepresented with respect to high frequencies. In our case, this leads to an ill-conditioned reconstruction problem. Therefore, a ramp filter is applied in the Fourier domain to give the proper weight to low and high frequencies.

The filtering is done before back-projection in every iteration, using a discretized version of the high pass ramp filter [murrell]. If there is only one iteration, the optimization then reduces to the standard filtered back-projection. This multiplication by a ramp filter can be seen as a preconditioner transforming an ellipsoidal (ill-posed) problem into a spherical problem, thus fastening the rate of convergence.

The fidelity term then becomes

(22)

where is the preconditioner. The gradient is

(23)

Here the operator is the filtered back-projection ; so the preconditioner is identified with the high-pass filtering process. We notice that the gradient with respect to the rings variables should also be filtered.

On the other hand, adding rings variables in the functional modifies the condition number of the problem. In order to prevent the reconstruction problem to be ill-posed, we introduced a multiplicative coefficient doing a variable substitution : becomes . Choosing enables to pre-condition the problem with respect to the ring variables, at the expense of an additional parameter.

4 Results


We present here some results for both simulated and real data, and compare our method with two mainstream techniques of rings correction : sinogram pre-processing based on Fourier-Wavelet de-striping [munch09] and image correction using polar coordinates transformation [prell09].

4.1 Simulated data

We use the standard test image “Lena” containing both smooth components and texture details. The tests are divided in increasing levels of difficulty for rings removal methods. First, constant lines are added in the sinogram. These lines independent from the projection angle give raise to rings artifacts in the reconstructed slice. Since the spurious lines have a constant value, they should be well handled by pre-processing techniques.

Then, lines with variable intensity are added to the sinogram. This makes the correction more difficult for pre-processing techniques, especially if the lines have sharp variations (i.e high frequencies components).

Finally, ring-shaped features are added in the phantom before adding spurious lines in the sinogram. This case is more challenging because correction methods should not remove any feature coming from the phantom (they belong to the “true” image), while actually removing rings coming from the sinogram (they come from a flaw in the experimental setup).

The reference sinogram pre-processing technique is the wavelet-Fourier filtering [munch09]. This methods first compute the wavelet decomposition at a level of the sinogram. The vertical detail coefficients at level

emphasize the spurious lines that give raise to rings artifacts. In these coefficients, a spurious line is nearly constant along the projection angle, thus, it has only low frequencies in the Fourier domain. Filtering these few low frequencies in the Fourier domain enable to suppress the line after taking the inverse Fourier transform. The filter used is a high-pass Gaussian filter whose standard deviation

tunes the bandwidth. Then, the sinogram is reconstructed from these filtered wavelet coefficients. The Matlab implementation of this method is provided in the author’s article.In the tests, denotes the standard deviation of the Gaussian filter and is the number of levels of the wavelet decomposition.

The image correction technique used here is Rings Correction in Polar Coordinates (RCP) [prell09]. It transforms the image into polar coordinates and performs a low-pass filtering in the radial direction. The filtered image is then subtracted from the original image, and a threshold is applied to ignore non-artifact structures. The result is filtered in the azimuthal direction. After a transformation into Cartesian coordinates, the image should only contain rings artifacts ; these are subtracted from the original image. A C++ implementation can be found at [blair14]. In the tests, the thresholding parameters are set so that all the image pixels can be considered as possibly part of an artifact. The important remaining parameters are the maximum ring width , and the maximum angular arc we expect the rings to have.

Figure 3 shows the results for the firsts test case. The rings are reduced by the sinogram pre-processing technique (3

.c), but they do not totally disappear. Besides, additional artifacts appear after the correction. The parameters of this method – wavelet name, decomposition levels and filter bandwidth – depend on the image type and the ring width/intensity. In particular, choosing a wavelet with many vanishing moments like in

[munch09] does not always yield to better results. The RCP performs slightly better (3.e) ; a strong artifact is added to the right but the result is qualitatively better. The Total Variation regularization entirely removes the rings (3.g).

Figure 4 shows the results for the second test case. The sinogram filtering adds many spurious rings (4.c). The RCP technique removes most of the rings, but a small rings details remain. Again, the difference image shows that strong artifacts are added (4.f) even if they are not obviously perceptible on the reconstructed image. The Total Variation regularization (4.g) entirely removes the rings artifacts ; the result is qualitatively very close to the original phantom.

A plot of the rings variables during the total variation reconstruction shows that the rings are actually captured by the ring vector . The six peaks representing the detected lines in the sinogram are clearly visible in Figure 2.

Figure 2: Vector of rings variables for second test case (Figure 4.b). The horizontal axis goes from zero to the number of bins of the detector – that is, in this simulated case, for the test image.

Figure 5 shows the results for the third test case. Here two features are added to the original phantom (5.a) : a black disk and a semi-circular “ring”. These features are part of the phantom, they should not be filtered by rings correction techniques. Lines added to the sinogram are not constant along the projection angle, and their width can be several pixels. This leads to a back-projected image (5.b) with large rings. The RCP technique (5.e) is more efficient than the sinogram preprocessing (5.c). The Total Variation rings removal slightly impacts the semi-circular white feature while preserving the black disk.

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 3: First test case.
(a) Original phantom. (b) Result of filtered backprojection after adding constant lines in the sinogram. (c) Image back-projected after applying the Munch et al. de-striper algorithm with , and the “Daubechies 15” wavelet. (d) Difference between the phantom and the corrected image. The PSNR is 26.6. (e) Result of the correction with the RCP technique with and . (f) Difference between the phantom and the corrected image. The PSNR is 29.6. (g) Result of the reconstruction using the Total Variation regularization in PyHST2, with parameters , . (h) Difference between the phantom and the corrected image. The PSNR is 39.0.

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 4: Second test case.
(a) Original phantom. (b) Result of filtered backprojection after adding lines of variable width and intensity in the sinogram. (c) Image back-projected after applying the Munch et al. de-striper algorithm, with , and the “Daubechies 15” wavelet. (d) Difference between the phantom and the corrected image. The PSNR is 29.2. (e) Result of the correction using the RCP technique with and . (f) Difference between the phantom and the corrected image. The PSNR is 25.1. (g) Result of the reconstruction using the Total Variation regularization in PyHST2, with parameters , . (h) Difference between the phantom and the corrected image. The PSNR is 39.4.

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 5: Third test case.
(a) Original phantom. (b) Result of filtered backprojection after adding lines of variable width and intensity in the sinogram. (c) Image back-projected after applying the Munch et al. de-striper algorithm, with , and the “Daubechies 20” wavelet. (d) Difference between the phantom and the corrected image. The PSNR is 23.2. (e) Result of the correction using the RCP technique, with and . (f) Difference between the phantom and the corrected image. The PSNR is 21.2 (g) Result of the reconstruction using the Total Variation regularization in PyHST2, with parameters , . (h) Difference between the phantom and the corrected image. The PSNR is 29.9.

Beside the visual aspect of the corrected image, we use the Peak Signal to Noise Ratio (PSNR) as a measure of the correction quality :

(24)

where is the mean square error between the two images and of dimensions and bit depth of bits. Although PSNR gives a score to the overall similarity between the corrected image and the original phantom, it is inconsistent with the eye perception of quality. For example, RCP performed better than sinogram filtering in these tests, but got a lowest PSNR for cases 2 and 3. The structural similarity index gives the same kind of results.

We were surprised that the sinogram pre-processing technique performed poorly on these simulated cases. The reason is probably that the intensity of the lines added to the sinogram is too large for this method. The technique presented in

[munch09] only filters the wavelet detail coefficients : if the sinogram line features are too large, they are captured by the approximation coefficients rather than the detail coefficients. By trying with lines of smallest amplitude, the Fourier-Wavelet method actually worked without adding big artifacts in the reconstructed image. For our particular experimental cases, however, the sinogram unwanted lines have non-negligible amplitudes, making this method ineffective.

4.2 Experimental data

We give here some results for reconstructions performed on real data. The samples were kindly provided by the ESRF beamline ID19.

4.2.1 Syntactic foam

The reconstruction technique was used on a syntactic foam sample. In this case, the rings are “large” in the extent that the radius difference between the exterior and the interior of the ring is several pixels. This means that the spurious lines in the sinogram have several pixel of width along the detector bins axis, forming “bands”. However, the intensities of the lines forming a band has too many variations to be efficiently filtered by sinogram pre-processing techniques. This case is also difficult for slice-correction algorithms which detect circular features, since the sample itself have circular features we do not want to remove.

(a) (b) (c)

Figure 6: Syntactic foam sample tomography acquired at ESRF ID19, with energy of 19 keV and pixel size of . (a) Filtered back-projection. (b) Correction with RCP technique, using the parameters and . (c) Reconstruction with Dictionary Learning technique in PyHST2, using the parameters , and .

4.2.2 Rhynie chert

We apply the reconstruction on a rhynie chert sample. This situation is almost the opposite of the previous case : the rings artifacts have a small intensity in the reconstructed slice, and the sample borders form a nearly circular polygonal shape. These border have a huge amplitude with respect to the rest of the sample, and the transition between the border and the interior/exterior is very sharp. Thus, slice correction techniques would try to remove the borders before any other feature in the slice depending on the thresholding parameters.

We realized that the rings correction was difficult for the total variation reconstruction : the procedure added rings tangent to one of the slice borders. It turned out that the problem was due to the rotation center for (back)projection improperly set, leading to accumulating errors in the iterative reconstruction. Indeed, total variation and dictionary learning reconstruction require to compute the projection for the functional, and the back-projection for the functional gradient. If the rotation center for these operations is not the same that the one used for actually rotate the sample, slight errors appear in the (back)projection ; these error accumulate with the number of iterations and take the form of circular features (Figure 7.c).

After setting the correct rotation center, we were able to remove the ring artifacts (Figure 7.d), especially the one near the center of Figure 7.a. In this case, the RCP technique did quite well (Figure 7.b).

(a) (b) (c) (d)

Figure 7: Rhynie chert sample tomography acquired at ESRF ID19, with energy of keV and pixel size of . (a) Filtered back-projection with the correct rotation axis. (b) Correction with the RCP technique, using the parameters and . (c) Reconstruction with the Total Variation regularization using the incorrect rotation axis. (d) Reconstruction with the Total Variation regularization using the correct rotation axis. The parameters were and .

4.3 Execution time and convergence rate

In this section we measure the execution time required to obtain an acceptable reconstruction. All the tests are performed on a machine with an Intel Xeon CPU E5-1607 v2 @ 3.00GHz processor and a GeForce GTX 750 Ti graphic card.

Figure 8: Execution time as a function of the number of projection. The image used is the test image “Lena” corrupted with the rings presented in the second test case.

We measured that the execution time is the same with rings correction and without rings correction, which is coherent since the ring correction is part of the functional. Hence, the execution time gives an insight of the time needed by Compressed Sensing reconstruction, but it is not the best choice to measure the influence of rings correction on the overall reconstruction time. We now focus on the number of iterations required to converge to an acceptable solution. The cost of a single iteration depends on the number of projections, as it can be guessed with Figure 8.

To measure the convergence rate, we use the values of the objective function. For Total Variation reconstruction, the objective function is given by (15), it includes both the fidelity term (Euclidean distance) and the regularization term (L1 norm of the image gradient). One can notice that the energy is not a monotonic function of the number of iterations, which can be seen as a drawback of FISTA with respect to ISTA.

(a) (b)

Figure 9: Energy as a function of the number of iterations for the Total Variation tomographic reconstruction. The energy is normalized by the energy of the last iteration in order to have the same scale in the two cases. The image used is the test image “Lena” corrupted with the rings presented in the second test case. (a) Evolution of energy with 800 projections. (b) Evolution of energy with 200 projections.

With rings correction, the reconstruction process takes more iterations to converge, as illustrated in Figure 9. For simulated and real data, it turned out that a satisfactory reconstruction can be achieved with less than iterations without rings correction. When the rings correction is activated, it takes about iteration to correctly remove the ring artifacts. Thus, while the rings correction has no additional cost per iteration, it takes nevertheless more iterations to converge to an image with removed ring artifacts. The “energy transfer” between the fidelity term and the L1 norm of the rings is actually quite slow.

The convergence rate also depends on the number of projections. Figure 9 shows that the reconstruction process converges in 500 iterations (400 with rings correction) for 200 projections, when it takes about 2000 iterations (1500 without rings correction) for 800 projections. It is due to the fact that the weight of fidelity term virtually increases as the number of projection increases. To counter-weight this, one has to increase the penalty term of the regularization, which makes the energy transfer between the fidelity term and the ring variables a little faster.

The reconstruction quality is still excellent with four times less projections – which is a crucial interest of Compressed Sensing tomographic reconstruction. Therefore, this method makes sense rather for low-dose tomography when there are few projections ; and for now it is quite expensive to use it only for rings correction when a reconstructed volume is already available.

4.4 Conclusions

We presented a new way to correct the rings artifacts that appear in tomographic reconstruction. Compressed sensing tomographic reconstruction is a promising technique which enables to obtain a slice of good quality with fewer projections. Including the rings artifacts correction in the iterative reconstruction process has shown to be efficient while requiring no extra pre or post-processing steps. Besides, additional artifacts are less likely to appear thanks to the regularization. This method can be adapted to any compressed sensing approach, since the only things to do are modifying the functional and the iterative correction step accordingly.

In a further work, we would like to improve the convergence rate of the rings correction in order to make it more attractive. We also would like to extend this method to lines that are not constant along the projection angle in the sinogram, in order to cover more general and difficult cases.

Acknowledgements We thank Elodie Boller and Paul Tafforeau from the ESRF beamline ID19 for having kindly provided experimental samples. We also thank Michel Desvignes from Gipsa-Lab, Grenoble, for his reviews and advices.

References