Controlled Wavelet Domain Sparsity in X-ray Tomography

TL;DR

This paper introduces an automatic method for selecting the regularization parameter in wavelet-based X-ray tomography, improving the sparsity control of the reconstructed image through a novel control algorithm.

Contribution

A new automatic parameter selection method for wavelet sparsity regularization in X-ray tomography using a control algorithm.

Findings

Effective control of wavelet coefficient sparsity.

Improved reconstruction quality with automatic parameter tuning.

Demonstrated robustness across different datasets.

Abstract

Tomographic reconstruction is an ill-posed inverse problem that calls for regularization. One possibility is to require sparsity of the unknown in an orthonormal wavelet basis. This in turn can be achieved by variational regularization where the penalty term is the sum of absolute values of wavelet coefficients. Daubechies, Defrise and De Mol (Comm. Pure Appl. Math. 57) showed that the minimizer of the variational regularization functional can be computed iteratively using a soft thresholding operation. Choosing the soft threshold parameter $\mu>0$ is analogous to the notoriously difficult problem of picking the optimal regularization parameter in Tikhonov regularization. Here a novel automatic method is introduced for choosing $\mu$, based on a control algorithm driving the sparsity of the reconstruction to an {\it a priori} known ratio of nonzero versus zero wavelet coefficients in the unknown function.

Tables

Table 1. Table 1: Initial values μ ( 0 ) superscript 𝜇 0 \mu^{(0)} of the sof-thresholding parameter.

		$120$ projections	$30$ projections
	Shepp-Logan	$0.0202$	$0.0195$
	Walnut	$0.0019$	$0.0021$

Table 2. Table 2: The relative error of the Shepp-Logan phantom reconstructions for FBP and CWDS.

	$120$ projections	$30$ projections
FBP	$0.15$	$0.27$
CWDS	$0.04$	$0.08$

Table 3. Table 3: Computation times for FBP reconstructions and CWDS reconstructions in seconds.

		$120$	$30$
walnut	FBP	$0.45$	$0.09$
	CWDS	$17.40$	$16.30$
Shepp-Logan	FBP	$0.02$	$0.01$
	CWDS	$98.90$	$29.50$

Equations

$${\boldsymbol{f}}_{\tiny{\text{S}}}=\operatornamewithlimits{argmin}_{{\boldsymbol{f}}\in{\mathbb{R}}^{N^{2}}}\Bigg{\{}\frac{1}{2}\|\mathbf{A}{\boldsymbol{f}}-{\boldsymbol{m}}\|_{2}^{2}+\mu\sum_{\gamma\in\Gamma}|\langle f,\psi_{\gamma}\rangle|\Bigg{\}}.$$

$${\boldsymbol{f}}_{\tiny{\text{S}}}=\operatornamewithlimits{argmin}_{{\boldsymbol{f}}\in{\mathbb{R}}^{N^{2}},\,{\boldsymbol{f}}>0}\Bigg{\{}\frac{1}{2}\|\mathbf{A}{\boldsymbol{f}}-{\boldsymbol{m}}\|_{2}^{2}+\mu\sum_{\gamma\in\Gamma}|\langle f,\psi_{\gamma}\rangle|\Bigg{\}}.$$

$$\begin{split}{\boldsymbol{y}}^{(i+1)}&=\mathbb{P}_{C}\Big{(}{\boldsymbol{f}}^{(i)}-\tau\nabla g({\boldsymbol{f}}^{(i)})-\lambda\mathbf{W}^{T}{\boldsymbol{v}}^{(i)}\Big{)}\\ {\boldsymbol{v}}^{(i+1)}&=\Big{(}I-\mathcal{T}_{\mu}\Big{)}\Big{(}\mathbf{W}{\boldsymbol{y}}^{(i+1)}+{\boldsymbol{v}}^{(i)}\Big{)}\\ {\boldsymbol{f}}^{(i+1)}&=\mathbb{P}_{C}\Big{(}{\boldsymbol{f}}^{(i)}-\tau\nabla g({\boldsymbol{f}}^{(i)})-\lambda\mathbf{W}^{T}{\boldsymbol{v}}^{(i+1)}\Big{)}\end{split}$$

$${\mathcal{T}}_{\mu }(c)=\{\begin{array}{ll}c+\frac{\mu }{2}& \text{if }x\le -\frac{\mu }{2}\\ 0& \text{if }|x|<\frac{\mu }{2}\\ c-\frac{\mu }{2}& \text{if }x\ge -\frac{\mu }{2}.\end{array}$$

$${\mu }^{(i+1)}={\mu }^{(i)}+\beta ({\mathcal{C}}^{(i)}-{\mathcal{C}}_{pr}),$$

$${\int }_{L_X}f(x)ds=\sum _{i=1}^N\sum _{j=1}^N{a}_{ij}{\tilde{f}}_{ij},$$

$$\mathbf{m}=\mathbf{A}\mathbf{f},$$

$$\varphi(x)\equiv 1,\qquad\psi(x)=\left\{\begin{array}[]{rl}1&\;\text{if}\quad 0\leq x<1/2,\\ -1&\;\text{if}\quad 1/2\leq x\leq 1.\end{array}\right.$$

$${\psi }_{jk}(x):=2^{j/2}\psi (2^{j}x-k)\text{mbox}forj\le 0,0\le k\le 2^j-1,$$

$${\phi }_{jk}(x):=2^{j/2}\phi (2^{j}x-k)\text{mbox}forj\le 0,0\le k\le 2^j-1,$$

$${\phi }_{jk}(x){\psi }_{jk}(y),{\psi }_{jk}(x){\phi }_{jk}(y),{\psi }_{jk}(x){\psi }_{jk}(y),$$

$$f⟶\mathcal{W}f=⟨f,{\psi }_{\gamma }⟩,\gamma \in \Gamma$$

$$\mathbf{W}\mathbf{f}\in {\mathbb{R}}^{N^2},$$

$${\boldsymbol{f}}_{\tiny{\text{S}}}=\operatornamewithlimits{argmin}_{{\boldsymbol{f}}\in{\mathbb{R}}_{+}^{N^{2}}}\Bigg{\{}\frac{1}{2}\|\mathbf{A}{\boldsymbol{f}}-{\boldsymbol{m}}\|_{2}^{2}+\mu\|\mathbf{W}{\boldsymbol{f}}\|_{1}\Bigg{\}}.$$

$${\#}_{\kappa }w:=\#\{i|1\le i\le N^2,|w_i|>\kappa \}.$$

$${\mathcal{C}}_{pr}=\frac{{\#}_{\kappa }\{\mathbf{W}{\mathbf{f}}_{pr}\}}{N^2},$$

$${\mu }^{(i+1)}:={\mu }^{(i)}+\beta ({\mathcal{C}}^{(i)}-{\mathcal{C}}_{pr}),$$

Full text

Keywords : tomography, wavelet, sparsity, regularization, control, limited data tomography, X-ray

1 Introduction

Tomographic imaging is based on recording projection images of an object along several directions of view. The resulting data can be interpreted as a collection of line integrals of an unknown attenuation coefficient function $f(x)$. In this work, we discretize the problem by approximating $f$ as a vectorized pixel image ${\boldsymbol{f}}\in{\mathbb{R}}^{N^{2}}$ and using the pencil-beam model for X-rays, so the indirect measurement is modelled by a matrix equation $\mathbf{A}{\boldsymbol{f}}={\boldsymbol{m}}$. The inverse problem of reconstructing ${\boldsymbol{f}}$ from tomographic data is highly sensitive to noise and modelling errors, or in other words ill-posed.

We focus on overcoming ill-posedness by enforcing sparsity of ${\boldsymbol{f}}$ in an orthonormal wavelet basis $\{\psi_{\gamma}\}_{\gamma\in\Gamma}$.

In practice, the sparse reconstruction ${\boldsymbol{f}}_{\tiny{\text{S}}}\in{\mathbb{R}}^{N^{2}}$ is defined as the minimizer of this variational regularization functional:

[TABLE]

The parameter $\mu$ in (1) describes a trade-off between emphasizing more the data fidelity term or the regularizing penalty term. In general, the larger the noise amplitude in the data, the larger $\mu$ needs to be.

One popular method to solve problem (1) is the so-called iterative soft-thresholding algorithm (ISTA). Such algorithm has been studied already in [1]; the adaptation to sparsity-promoting inversion was introduced in [2] and further developed in [3]. Nevertheless, convergence rate for a constrained problem, such as non-negativity constraints, is not taken into account in [2, 3]. However, in tomographic problems, enforcing non-negativity on the attenuation coefficients is highly desired. This is based on the physical fact that the X-ray radiation can only attenuate inside the target, not strengthen. Thus, the problem we need to solve reads as:

[TABLE]

where the inequality ${\boldsymbol{f}}>0$ is meant component-wise. In their seminal paper [4], Peijun Chen, Jianguo Huang, and Xiaoqun Zhang show that the minimizer of (1) can be computed using the primal-dual fixed point (PDFP) algorithm:

[TABLE]

where $\tau$ and $\lambda$ are positive parameters, $g({\boldsymbol{f}})=\frac{1}{2}\|\mathbf{A}{\boldsymbol{f}}-{\boldsymbol{m}}\|_{2}^{2}$, the matrix $\mathbf{W}$ is a digital implementation of the wavelet transform and $\mathcal{T}$ is the soft-thresholding operator defined by

[TABLE]

Here $\mu>0$ represents the thresholding parameter, while $\tau$ and $\lambda$ are parameters that needs to be suitably chosen to guarantee convergence. In detail, $0<\lambda<1/\lambda_{\max}(\mathbf{W}\mathbf{W}^{T})$, where $\lambda_{\max}$ denotes the maximum eigenvalue, and $0<\tau<2/\tau_{\text{lip}}$, being $\tau_{\text{lip}}$ the Lipschitz constant for $g({\boldsymbol{f}})$. Furthermore, in (3) the non-negative “quadrant” is denoted by $C={\mathbb{R}}_{+}^{N^{2}}$ and $\mathbb{P}_{C}$ is the euclidian projection. In other words, $\mathbb{P}_{C}$ replaces any negative elements in the input vector by zero.

Choosing the soft-thresholding parameter $\mu$ is analogous to the notoriously difficult problem of picking the optimal regularization parameter in Tikhonov regularization. Many approaches for the regularization parameter selection have been proposed. For a selection of methods designed for total variation (TV) regularization see the following studies: [5, 6, 7, 8, 9, 10, 11, 12]. In this paper we introduce a novel automatic method for choosing $\mu$ based on a control algorithm driving the sparsity of the reconstruction to an a priori known ratio $0\leq{\mathcal{C}_{pr}}\leq 1$ of nonzero wavelet coefficients in ${\boldsymbol{f}}$. Our approach is based on the following idea: in sparsity-promoting regularization, it is natural to assume that the a priori information is given as the percentage of nonzero coefficients in the unknown. The idea of using the a priori known level of sparsity has been used previously [13, 14], however the idea of using feedback control to achieve this is new.

We think of the iteration (3) as a plant which takes the current threshold parameter $\mu^{(i)}$ as an input and returns $\mathcal{C}^{(i)}$, the level of sparsity in the iterate ${\boldsymbol{f}}^{(i)}$, as an output. Then, we apply a simple incremental feedback control to $\mu^{(i)}$. The feedback loop we propose is inspired by the proportional-integral-derivative (PID) controllers, which are widely used to control industrial processes [15, 16, 17]. If $\mathcal{C}_{pr}$ is the expected degree of sparsity, and $\mathcal{C}^{(i)}$ is the degree of sparsity at the current iterate, we change $\mu^{(i)}$ adaptively as follows:

[TABLE]

where $\beta>0$ is a parameter used to tune the controller. We propose a simple method for choosing $\beta$ based on the wavelet coefficients of the backprojection reconstruction, which is quick and easy to compute. If $\beta$ is chosen too large, then the controller results in an oscillating behavior for the sequence $(\mu^{(i)})_{i}$. On the other hand, if $\beta$ is chosen too small, reaching the expected sparsity level may take a long time. Therefore we also account for an additional fine-tuning of the controller by exploiting the zero-crossings of the controller error $e^{(i)}=\mathcal{C}^{(i)}-\mathcal{C}_{pr}$.

We test our fully automatic controlled wavelet domain sparsity (CWDS) method on both simulated and real tomographic data. The results suggest that the method produces robust and accurate reconstructions, when the suitable degree of sparsity is available.

CWDS has a connection to the following studies, which also use a parameter changing adaptively during the iterations: [18, 19, 20, 21, 22, 23, 24]. However, our approach is different from all of them as it promotes an a priori known level of sparsity. Also, this is not the first study which uses the wavelet transform as a regularization tool in limited data tomography. A non-exhaustive list includes [25, 26, 27, 28, 14, 29, 30]. However, the proposed approach is different from the previous works, since it promotes a fully automatic choice for the regularization parameter.

2 Materials and Methods

2.1 Tomography setup

Consider a physical domain $\Omega\subset{\mathbb{R}}^{2}$ and a non-negative attenuation function $f:\Omega\subset{\mathbb{R}}^{2}\rightarrow{\mathbb{R}}_{+}$. As outlined in the Introduction, we represent $f$ by a matrix ${\boldsymbol{f}}=[\tilde{f}_{ij}]\in{\mathbb{R}}^{N\times N}$ that is later on intended as a vector belonging to ${\mathbb{R}}^{N^{2}}$, obtained by stacking the entries of the matrix column by column. In X-ray tomography, the detector measures the incoming photons and the measurement data are collected from the intensity losses of X-rays from different directions or angles of view. After calibration, the measurements can be modeled as

[TABLE]

where $a_{ij}$ is the distance that a X-ray line $L_{X}$ travels through the pixel $(i,j)$. This results in the following matrix model:

[TABLE]

where the measurement matrix $\mathbf{A}=[a_{ij}]\in{\mathbb{R}}^{P\times N^{2}}$ contains the information about the measurement geometry, and ${\boldsymbol{m}}\in{\mathbb{R}}^{P}$ is the vector representing the measured data (also called sinogram), $P$ being the number of angles of view multiplied by the number of detector cells.

Notice that, in the following, we assume both the measurement matrix $\mathbf{A}$ and the measured data ${\boldsymbol{m}}$ to be normalized by the norm $\|\mathbf{A}\|$ of the matrix $\mathbf{A}$.

2.2 2D Haar wavelets

For the readers’ sake of convenience, we briefly recall here the main ideas about Haar wavelets.

Consider the two real-valued functions $\varphi(x)$ and $\psi(x)$ defined on the interval $[0,1]$. Generally, $\varphi(x)$ is referred to as scaling function and $\psi(x)$ as mother wavelet. They are defined as follows:

[TABLE]

A Haar wavelet system is built by appropriately scaling and translating the mother wavelet $\psi(x)$:

[TABLE]

and the scaling function $\varphi(x)$:

[TABLE]

where $\varphi(x)=0$ for $x<0$ and $x>1$. Here, $j,\,k\in\mathbb{Z}$.

It is well known that the above 1D construction leads to an orthonormal system. In 2D, we consider the standard tensor-product extension of the 1D Haar wavelet transform. In detail, a 2D Haar system is spanned by four types of functions. Three of these types have the following form:

[TABLE]

and the fourth type is given by $\varphi_{j_{0}k}(x)\varphi_{j_{0}k}(y)$. Notice that the fourth type describes the coarsest scale $j_{0}$. The associated wavelet transform of a function $f$ is given by

[TABLE]

where $\langle f,\psi_{\gamma}\rangle$ denotes the so-called wavelet coefficients. Here, for notational convenience, we use the index $\gamma\in\Gamma$ to combine together the three types (7) at several scales $j$ and locations $k$, and the fourth type at several locations.

In the following, we are interested in the digital setting, i.e., we consider the matrix underlying the wavelet transform, which we shall denote by $\mathbf{W}\in{\mathbb{R}}^{N^{2}\times N^{2}}$. If ${\boldsymbol{f}}\in{\mathbb{R}}^{N^{2}}$, the vector collecting all the wavelet coefficients is given by:

[TABLE]

where it is clear that the matrix product $\mathbf{W}{\boldsymbol{f}}$ is the digital counterpart of (8). With the above notation, the minimization problem (2) reads as

[TABLE]

One of the main benefit of wavelets is that the transform coefficients are easy to compute and many fast algorithmic implementation are available.

For more information about the Haar wavelet transform, and its implementation, we refer to the classic text [31].

2.3 Sparsity promoting-regularization

We consider the functional in (2) with $\{\psi_{\gamma}\}_{\gamma\in\Gamma}$ being the Haar wavelet basis as described in Subsection 2.2. To solve the minimization problem (2), we implement the PDFP algorithm (3).

2.4 Sparsity selection

We assume that we have available an object ${\boldsymbol{f}}_{pr}$ similar to the one we are imaging.

Given $\kappa\geq 0$, for a vector $w\in{\mathbb{R}}^{N^{2}}$ we define the number of elements larger than $\kappa$ in absolute value as follows:

[TABLE]

Now, the prior sparsity level is defined by

[TABLE]

where $N^{2}$ is the total number of coefficients. In practical computations the value of $\kappa$ is set to be small but positive.

2.5 Automatic selection of the soft-thresholding parameter $\mu$

Assume that we know a priori the expected degree of sparsity $0<\mathcal{C}_{pr}\leq 1$ in the reconstruction. We introduce a simple feedback loop to drive the soft-thresholding parameter $\mu$ to the desired ratio $\mathcal{C}_{pr}$ of nonzero wavelet coefficients.

The core idea is to allow $\mu=\mu^{(i)}$ to vary during the iterations by adaptively tuning it at each iteration by the following updating rule:

[TABLE]

where $0\leq{\mathcal{C}}^{(i)}\leq 1$ is the sparsity level of the reconstruction ${\boldsymbol{f}}^{(i)}$ at the $i$-th iteration. The above controller is a special case of an incremental PID-controller, where only integral control is performed.

2.6 The tuning parameter $\beta$

Selecting the tuning parameter $\beta$ is easier than selecting the soft-thresholding parameter $\mu$. Indeed, $\beta$ has to be small enough to avoid oscillations in the sparsity $\mathcal{C}^{(i)}$ of the iterates as a function of $i$. If $\beta$ is chosen too small, this only result in a slower convergence of the algorithm.

To this purpose, we choose $\beta$ by making a suitable guess for the initial $\mu^{(0)}$. First, we compute the back-projection of the measured data to get a rough reconstruction. Back-projection is quick to compute and shows the dominant features of the target, but noise and artefacts are still predominant. As a result, the back-projection reconstruction is only good enough for estimating an initial guess for $\mu^{(0)}$, which is done by computing its wavelet coefficients. The initial value of the thresholding parameter $\mu^{(0)}$ is set equal to the mean of the absolute values of the $M$ smallest wavelet coefficients. In our case, we choose $M=n\,(1-\mathcal{C}_{pr})$, where $n$ is the total number of wavelet coefficients. Lastly, the tuning parameter is set to be $\beta=\omega\mu^{(0)}$, where $\omega$ is a positive parameter. To start with a small value of $\beta$, $\omega$ is required to be small, and vice versa.

In addition, the controller is fine tuned by detecting when the sign of difference $e^{(i)}=\mathcal{C}^{(i)}-\mathcal{C}_{pr}$ changes. When this happens, $\beta$ is updated by $\beta|e^{(i)}-e^{(i-1)}|$. The underlying idea is that, if the desired sparsity level is crossed, that is, $e$ changes sign, either $\beta$ is far too large and oscillations have emerged, or we are already reasonably close to the optimal $\mu$ and $\beta$ can be decreased without affecting the performance too much.

2.7 Pseudo-algorithm

A step-by-step description of the proposed CWDS algorithm is summarized in Algorithm 1.

3 Data Acquisition

In this paper, we consider both simulated data (see Section 3.1) and real data (see Section 3.2).

3.1 Simulated data

We use the Shepp-Logan phantom, available, for instance, in the Matlab Image Processing toolbox (see Figure 1). The phantom is sized $N\times N$, with $N=328$. The projection data (i.e., sinogram) of the simulated phantom is corrupted by a white Gaussian process with zero mean and $0.1\%$ variance.

3.2 Real data

We use the tomographic X-ray real data of a walnut, consisting of a 2D cross-section of a real 3D walnut measured with a custom-built CT device available at the University of Helsinki (Finland). The dataset is available and freely downloadable at http://fips.fi/dataset.php. For a detailed documentation of the acquiring setup, see [32]. Here we only mention that the sinogram is sized $328\times 120$. Sinograms with different resolutions for the angle of view can be obtained by further downsampling.

4 Numerical Experiments

In this Section, we present preliminary numerical results in the framework of 2D fan-beam geometry.

4.1 Algorithm parameters

In all the experiments, we set $\lambda=0.99$ (being $\lambda_{\max}(\mathbf{W}\mathbf{W}^{T})=1$) and $\tau=1$ to ensure convergence. Also, we choose $\epsilon_{1}=5\times 10^{-4}$ and $\epsilon_{2}=5\times 10^{-4}$ for the stopping rule, and $I_{\max}=1500$ as a safeguard maximum number of iterations (which is never attained in the results reported in Section 4.3), $\beta=\omega\mu^{(0)}$, where $\omega=1$ and the values of $\mu^{(0)}$ for each experiments are shown in Table 1.

All the algorithms were implemented in Matlab 8.5 (R2015a) and performed on Intel Core i5 at 2.9 GHz and CPU 8GB 1867 MHz DDR3 memory. The Haar matrix $\mathbf{W}$ is generated by using Spot–A Linear-Operator Toolbox [33]. The number of scales for the wavelet transform is set equal to 3 (see Figure 3).

4.2 A priori sparsity level

To compute the desired sparsity level, we choose $\kappa=10^{-6}$ for both the Shepp-Logan phantom and the walnut, and we apply the strategy outlined in Section 2.5. In particular, for the walnut case, since we do not have at disposal the “original” target, we compute the sparsity level from the photographs of two walnuts cut in half (see Figure 2). The a priori sparsity level $\mathcal{C}_{pr}$ for the walnut is the average of those two sparsity levels.

For the Shepp-Logan phantom, the percentage of nonzero coefficients was estimated to be $12\%$. The percentage of the nonzero coefficients for the walnut case was estimated to be $32\%$.

4.3 Reconstruction results

In this Section, we present numerical results for the CWDS method, using both simulated and real data. As a benchmark comparison, filtered back-projection (FBP) reconstructions were also computed. For both simulated and real data, we computed reconstructions for two different resolutions of the angle of view, namely 120 and 30 projection directions, respectively.

The reconstructions of the Shepp-Logan phantom are shown in Figure 4. Plots of the sparsity levels, as the iteration progresses, are reported in Figure 6. For the 120 projections case, the proposed approach converges in 885 iterations, while, in the 30 projections case, it converges in 301 iterations. As figure of merit, we use the relative error: the obtained values are summarized in Table 2, where we also report the values of the relative error obtained for the FBP reconstructions.

The reconstructions for the walnut dataset, for both 120 and 30 projections, are collected in Figure 5. The corresponding sparsity plots are shown in Figure 8. Concerning the number of iterations to convergence, the 120 projections case required 180 iterations, while in the 30 projections case convergence was reached in 206 iterations.

Lastly, the computation times for all the reconstructions are reported in Table 3.

5 Discussion

We presented results for both simulated and real X-ray data, also in the limited data case of only 30 projection views, with the fully automatic CWDS method. As it can be seen in Figures 4 and 5, the reconstructions for both the Shepp-Logan phantom and the walnut data outperform the FBP reconstructions. For the Shepp-Logan case, this is confirmed by the relative errors reported in Table 2. In detail, the reconstructions using CWDS produce sharper images, with less artefacts. Overall, the quality of the reconstruction remains good even when the number of projections is reduced to 30, while, for the FBP reconstructions, streak artefacts overwhelms the reconstructions. Finally, the presence of $\ell_{1}$-norm term combined with a sparsity transform, that produce denoising, and the non-negativity constraint (which is not enforced in the classical FBP scheme) definitively improves the reconstructions.

Concerning the behavior of the sparsity level for the walnut case, it can be seen in the first row of Figure 8 that the initial rapid oscillations decays fast. This is due to the role of the additional controller tuning $\beta$, as presented in Subsection 2.6.

For the Shepp-Logan case, it can be seen in Figure 6, and with a closer look for some iterations in Figure 7, that the ratio of nonzero wavelet coefficients produces small oscillations for many iterations. In fact, this is a behavior that can appear with the proposed method: if the controller error $e$  changes sign but the absolute difference of the error of the two consecutive iterations is small, there is very little change in $\beta$. However, in the long run, the oscillations disappear as $\beta$ is slowly decreased.

Future research could delve into alternative adaptive self-tuning controllers, such as the adaptive integral controller introduced in [34]. Such controllers might improve the system response to unexpected disturbances and help with the oscillations caused by the slow decay of $\beta$ demonstrated in Figure 7. Additionally careful analysis of the dynamics of the algorithm (3) is required to see if convergence of CWDS can always be guaranteed with the methods presented in this paper.

Anyhow, what is remarkable is that, for all numerical experiments, the sparsity level eventually converges to the desired sparsity level $\mathcal{C}_{pr}$.

6 Conclusion

In this paper, we proposed a new approach in tuning the regularization parameter, in this case the sparsity level of the reconstruction in the wavelet domain. CWDS seems to be a promising strategy, especially in real life applications where the end-users could avoid manually tuning the parameters.

In the case of sparsely collected projection data, the fully automatic CWDS outperforms the conventional FBP algorithm in terms of image quality (measured as relative RMS error).

Acknowledgments

This work was supported by the Academy of Finland through the Finnish Centre of Excellence in Inverse Problems Research 2012–2017 (Academy of Finland CoE-project 284715).