Blind deconvolution in model-based iterative reconstruction for CT using a normalized sparsity measure

Model-based iterative reconstruction techniques for CT that include a description of the noise statistics and a physical forward model of the image formation process have proven to increase image quality for many applications. Specifically, including models of the system blur into the physical forward model and thus implicitly performing a deconvolution of the projections during tomographic reconstruction, could demonstrate distinct improvements, especially in terms of resolution. However, the results strongly rely on an exact characterization of all components contributing to the system blur. Such characterizations can be laborious and even a slight mismatch can diminish image quality significantly. Therefore, we introduce a novel objective function, which enables us to jointly estimate system blur parameters during tomographic reconstruction. Conventional objective functions are biased in terms of blur and can yield lowest cost to blurred reconstructions with low noise levels. A key feature of our objective function is a new normalized sparsity measure for CT based on total-variation regularization, constructed to be less biased in terms of blur. We outline a solving strategy for jointly recovering low-dimensional blur parameters during tomographic reconstruction. We perform an extensive simulation study, evaluating the performance of the regularization and the dependency of the different parts of the objective function on the blur parameters. Scenarios with different regularization strengths and system blurs are investigated, demonstrating that we can recover the blur parameter used for the simulations. The proposed strategy is validated and the dependency of the objective function with the number of iterations is analyzed. Finally, our approach is experimentally validated on test-bench data of a human wrist phantom, where the estimated blur parameter coincides well with visual inspection. Our findings are not restricted to attenuation-based CT and may facilitate the recovery of more complex imaging model parameters.

. Second simulation study. The different parts of the objective function using the normalized sparsity measure are plotted over the blur parameters for the seven different measurements simulated with different blur parameters. The vertical dashed lines correspond to the blur parameters used in the respective simulations. In (a) the value of the data-fidelity is shown. The value of the normalized sparsity measure is shown in (b) and (c) depicts the overall objective function value. . First simulation study. In (a)-(c) the different parts of the objective function using conventional TV regularization are depicted over the blur parameter for different regularization strengths. In (a), the value of the data-fidelity is shown, in (b) the value of the regularization term and (c) shows the value of the objective function. In (d)-(f) the same information is shown for the proposed objective function using the normalized sparsity measure. In (g)-(k) reconstructions for different blur parameters are shown. The left sides of each figure show the reconstructions obtained with conventional TV regularization, while the right halves have been reconstructed using the normalized sparsity measure. In addition, the upper halves are shown in a more narrow grayvalue window than the lower halves.

Introduction
There is increasing interest in model-based iterative reconstruction (MBIR) techniques which include a physical model of the image formation process as well as a mathematical formulation of the noise characteristics (Nuyts et al 2013). The model accounting for the noise characteristics can, for instance, be derived from a Poisson distribution  or a multivariate normal distribution, which is also capable of modeling noise correlations (Tilley et al 2016a). The physical model is not limited to the description of the exponential decay in intensity. In particular, various effects that add to the system blur, and thus degrade the image quality, can be included, such as the finite source size or the spread of the scintillator (Nuyts et al 2013). MBIR algorithms that neglect blur have demonstrated an ability to outperform conventional approaches, such as analytical reconstruction methods (Feldkamp et al 1984), showing improved image quality or significantly reduced patient

Abstract
Model-based iterative reconstruction techniques for CT that include a description of the noise statistics and a physical forward model of the image formation process have proven to increase image quality for many applications. Specifically, including models of the system blur into the physical forward model and thus implicitly performing a deconvolution of the projections during tomographic reconstruction, could demonstrate distinct improvements, especially in terms of resolution. However, the results strongly rely on an exact characterization of all components contributing to the system blur. Such characterizations can be laborious and even a slight mismatch can diminish image quality significantly. Therefore, we introduce a novel objective function, which enables us to jointly estimate system blur parameters during tomographic reconstruction. Conventional objective functions are biased in terms of blur and can yield lowest cost to blurred reconstructions with low noise levels. A key feature of our objective function is a new normalized sparsity measure for CT based on total-variation regularization, constructed to be less biased in terms of blur. We outline a solving strategy for jointly recovering low-dimensional blur parameters during tomographic reconstruction. We perform an extensive simulation study, evaluating the performance of the regularization and the dependency of the different parts of the objective function on the blur parameters. Scenarios with different regularization strengths and system blurs are investigated, demonstrating that we can recover the blur parameter used for the simulations. The proposed strategy is validated and the dependency of the objective function with the number of iterations is analyzed. Finally, our approach is experimentally validated on test-bench data of a human wrist phantom, where the estimated blur parameter coincides well with visual inspection. Our findings are not restricted to attenuation-based CT and may facilitate the recovery of more complex imaging model parameters. dose (Thibault et al 2007, Noël et al 2013. However, including models for the system blur directly into the forward model has the potential to further increase image quality, especially in terms of resolution (Tilley et al 2018).
Prior efforts to model and correct for the system blur have concentrated largely on preprocessing steps prior to tomographic reconstruction (La Rivière et al 2006, Riviere and Vargas 2008, Zhang et al 2014. However, such preprocessing steps, at best, preserve information. In addition, noise properties have to be propagated through these preprocessing steps in order to leverage them in subsequent (iterative) tomographic reconstructions, which often requires additional simplifications (Tilley et al 2016a).
By directly integrating the system blur into the physical model, the deblurring step is simultaneously performed during reconstruction with regularization in the image domain. For instance, total-variation (TV) regularization (Rudin et al 1992) is preferably applied in the image domain, where the assumption of having piece-wise constant structures is more valid than in the projection domain. In addition, such modelling directly accounts for the noise characteristics of the measurements in the projection domain (Tilley et al 2018).
MBIR methods, that directly incorporate the blur into the forward model have already been used in SPECT (Yu et al 2000, Feng et al 2006, digital breast tomosynthesis (Zheng et al 2018), and flat-panel cone-beam CT (Tilley et al 2016a, Hashemi et al 2017, Tilley et al 2018. However, these techniques require a precise characterization of all components that contribute to the system blur. Overestimation of blur can yield ringing artifacts and noise amplifications, while residual blur in reconstructions is obtained when the system blur is underestimated. Such characterizations can be difficult and time consuming, especially when the system blur changes over time. For instance, the influence of the x-ray focal spot can change with system geometry and tube current, or for different samples, due to depth dependency of the blur (Riviere and Vargas 2008, Behling 2015, Tilley et al 2016b. High-resolution CCD detectors with optics have to be aligned, resulting in a varying detector response. Other limiting factors are more complicated to characterize. For example, blur can be caused by the elements of the reconstruction algorithm, e.g. interpolation in projection operations or the regularization (Tward and Siewerdsen 2008, Gang et al 2014, Hashemi et al 2017. Because there are penalties for improper modeling, it may be advantageous to parameterize and jointly estimate the blur model. This is similar to blind deconvolution in optical imaging, where algorithms exist that simultaneously estimate the system blur and recover the underlying image without the need of prior characterization (Chan and Wong 1998). These algorithms employ similar regularization techniques as those found in MBIR, making them a promising approach for the transfer to tomographic reconstruction.
In the following, we establish a MBIR algorithm, which includes a parameterized blur model to account for the system blur. We introduce a novel sparsity measure for CT, which targets unbiased reconstructions in terms of blur. This enables us to formulate an objective function, which can be used to jointly estimate the blur parameters during tomographic reconstruction. We outline an efficient optimization strategy for jointly estimating lowdimensional blur parameters with the proposed objective function. Finally, we validate and discuss our findings in a simulation study and apply our approach to physical bench data for a low-dimensional blur model.

MBIR using a parameterized blur model
We use an approach similar to Tilley et al (2018) and replace the static blur operator by a parameterized version. The measurement model can then be written as where ȳ is a vector of mean measurements and µ describes the distribution of the linear attenuation coefficients. In addition, A models the x-ray transform, the exponential function operates element-wise, I 0 is the (in this case scalar) source flux and B is a matrix accounting for the system blur, which depends on some blur parameters σ. The objective function L consists of a data fidelity term that is derived from a multivariate normal distribution and a regularization term, and has the following form We denote y as the projection measurements and W approximates the inverse of the covariance matrix, which is estimated from the measurements according to W = D[y −1 ] with D denoting a diagonal matrix with the vector in brackets as its diagonal elements. This approximation holds if the measurements are Poisson distributed and uncorrelated. Noise correlations or additive noise can also, in principle, be included in the inverse of the covariance matrix (Tilley et al 2016a). As the tomographic reconstruction is an ill-conditioned problem, a regularization term is included, which incorporates prior knowledge about the sample. This term consists of a scalar coupling parameter β, specifying the relative strength of the regularization, and a regularizer R ψ . Although various regularization strategies exist (Zhang et al 2018), pair-wise Gibbs priors that penalize the differences in neighboring voxel values, are commonly used for CT applications. This is motivated by penalizing high noise levels as they are predominantly present in the high-frequency components. In general, such a regularizer can be written as where the index i iterates over all entries of µ. The index n loops over all spatially neighboring voxels of i (usually 26) denoted by the set N i . The distance to the corresponding neighbors is weighted by ∆ i,n and ψ is a (scalar) penalty function. One commonly used choice is the total variation (TV) regularization R 1 , which employs the 1 norm as penalty function. TV arises from the assumption of the object being piece-wise constant. This reduces the noise content, while additionally imposing sparse solutions. Tomographic reconstruction algorithms that employ TV regularization have been very successful for CT applications (Sidky and Pan 2008), especially in the field of compressed sensing.
However, regularization strategies, which penalize the difference in neighboring voxel values are biased in terms of blur, as the high-frequency components are attenuated and thus decrease the response of any derivativetype filter. Consequently, the regularization values of blurred reconstructions (obtained with less deblurring) are lower than for sharp reconstructions (with more deblurring), as the noise levels of blurred reconstructions are in general lower.

Blind deconvolution of photographs
For digital processing of images, TV regularization is also widely used (Rudin et al 1992, Levin et al 2011, Ahmed et al 2014. TV-based methods have become popular in the field of noise removal (Rudin et al 1992) and have then been extended to blind deconvolution (Chan and Wong 1998).
Many new algorithms that adapt the framework of Chan and Wong (1998) have been developed since. Though these approaches have been observed to have a degree of success, it has also been shown that any algorithm that minimizes the respective objective functions provides suboptimal results, as the desired solution has a higher cost than the 'no-blur' solution (Levin et al 2011). This behavior results from bias of the regularization that favors noiseless and blurry images over noisy and sharp ones. An analysis of the original approach of Chan and Wong (1998) has been performed in Perrone and Favaro (2014), where they show how subtle choices in the optimization procedure avoid the 'no-blur' solution, despite having lowest cost. Various optimization methods have been developed to yield the desired solution, even though those solutions do not have the lowest cost, such as marginalization of the cost function, adaptive cost functions, alpha-matte extraction or edge location (Krishnan et al 2011).
Ideally, we would like to construct a metric that captures prior information about the sample such as sparsity without biasing the blur model. With respect to figure 1, we sketch how such a metric can be constructed. The 2 norm as depicted in figure 1(a) for a 2D problem only penalizes the magnitude (centrical contour lines) favoring blurred solutions and thus only biasing the solution without enforcing sparsity. The 1 norm as depicted in figure 1(b) enforces sparsity. For the same magnitude (distance from the origin), sparse solutions are favored, where one of the components is small (near the axes). However, globally, small magnitudes are favored, which again biases the blur model towards blurred and noiseless solutions. In Krishnan et al (2011) a new metric was proposed called the normalized sparsity measure. They propose to scale the 1 norm by the 2 norm to be scaleindependent-enforcing sparsity at all magnitudes. As depicted in figure 1(c) minimum penalties are located along the axis with no dependence on the magnitude (note radial contour lines). Practically, using such a modified regularizer permits enforcement of edges (via sparsity at pixel differences) but those edges may be created at all contrast levels-unlike, for example, a 1 penalty that increasingly penalizes large pixel differences.
This regularization approach has been successfully applied for blind image deconvolution of photographs (Krishnan et al 2011). In that work, two regularization terms are employed. The normalized sparsity measure is applied to the gradient images (the concatenation of the forward difference in both dimensions) and a TV regularization is applied on the blur kernel. In a first step the blur kernel is estimated and then a second non-blind deconvolution is used to recover the desired images.

Tomographic reconstruction using a normalized sparsity measure
Although regularization techniques and the reconstruction approaches in blind-deconvolution of photographs differ significantly from those techniques applied in CT, we use the underlying concepts of Krishnan et al (2011) to construct a TV regularizer for CT, which does not encourage blurry solutions.
For normalization, we employ quadratic regularization which uses the 2 2 -norm in the penalty function. To correctly account for the reduction in magnitude imposed by TV regularization, the square-root of the quadratic normalization is used according to We refer to S as the normalized sparsity measure for CT. This measure is designed to capture the prior information about the sample being sparse (piece-wise constant), without reducing the magnitude and thus imposing less bias on the solution. However, substituting the conventional regularizers R with the sparsity measure S in the objective function poses a challenge as the regularization term is non-convex. As proposed in Krishnan et al (2011), one can fix the denominator of the regularizer from the previous iterate. We adapt that approach here. We can interpret the objective function with conventional TV regularization R 1 and varying regularization strength. The normalization of the sparsity measure is attributed to the regularization strength according to defining the denominator of the normalized sparsity measure as Algorithm (1) outlines the generic solving strategy to solvê The optimization algorithm requires the projection measurements y, a non-zero start image µ (0) , the blur parameters σ, the regularization strength β, and the total number of iterations N. In every iteration n ∈ {0, ..., N − 1}, the normalization γ is computed from the current estimate µ (n) . Afterwards, the update ∆µ is calculated according the the objective L γ , using a gradient-based optimization algorithm to obtain the new estimate µ (n+1) . This process is repeated until the total number of iterations N is exceeded.

Blind deconvolution reconstruction
To simultaneously reconstruct the volume and perform a blind deconvolution, the objective has to be jointly estimated. We denote this as Figure 1. Comparison of different norms for a two-dimensional (2D) problem. The potential is color-coded and the contour lines are drawn in white. In (a) the 2 norm is depicted, which only penalizes the magnitude and is therefore scale-dependent. In (b) the corresponding 1 norm is shown. The magnitude is again penalized, also making this norm scale-dependent. However, at the same distance from the origin, sparse solutions (located at the edges) are preferred. In (c) we show the 1 norm normalized by the corresponding 2 norm. This normalized measure has its minima along the axis, therefore enforcing sparsity again. However, there is no penalty on the magnitude.
and we refer to this as blind deconvolution CT reconstruction. We note that the normalized sparsity measure does not directly depend on the blur parameter σ, but instead indirectly via the volume µ. As a result, there is no gradient information of the normalized sparsity measure available. Thus with a pure gradient-based method, the correct blur parameters cannot be recovered from the prior knowledge imposed by the regularization. However, if we restrict ourselves to a low dimensional problem, e.g. σ = σ is a scalar and define which is a one-dimensional objective function, we may then solvê to estimate the blur parameter. This approach is not efficient in practice, as a whole (iterative) tomographic reconstruction has to be performed to yield L γ (σ) at any point. Therefore, we also propose the following as a first step towards a computationally more efficient approach. The key idea is that we do not compute L γ (σ) globally, but only evaluate the function locally and update σ after a certain amount of iterations. The details are outlined in algorithm 2. In contrast to the previous optimization algorithm, we now require only a rough estimation of the blur parameter as initialization σ (0) . In addition, the number of iterations, M, and the step size of the blur parameter, ∆σ must be provided. In every iteration m ∈ {0, ..., M − 1}, three reconstructions are performed in parallel for N intermediate iterations (see algorithm 1), one with the blur parameter of the initial guess, one with a blur parameter reduced by ∆σ and one with a higher blur parameter increased by the same amount. Subsequently, the new estimate of the volume and the blur parameter are chosen according to the lowest objective value of the three reconstructions. In the next iteration, three new reconstructions distributed around the new estimate of the blur parameter are performed initialized with the new estimate of the volume. This scheme is repeated M times.
Algorithm 2. Optimization algorithm for blind deconvolution reconstruction.
In practice, we chose the parameters N and M such that their product N · M matches the total number of iterations used for a reconstruction with known blur parameter, because the reconstructions in the inner loop do not have to be performed until convergence. Thus for blind deconvolution CT reconstruction, three times as many operations have to be performed. In most cases this will result in a tripling of the reconstruction time. However, in principle, the three reconstructions in the inner loop can be performed in parallel. In this scheme, the precision of σ is predefined, although one could devise an adaptive method, where the stepsize decreases with iterations. In addition, the number of intermediate iterations, N, must not be too small, such that it is biased by the initial guess of the previous iteration. However, more intermediate iterations also increase the computational cost.

Evaluation of the regularization strategies in a simple 3D blur example
To gain intuition for 3D regularization in CT and to explore the tradeoffs of the normalized sparsity measure versus the conventional TV regularization, we introduce a simple 3D deblurring problem. Specifically, we use the measurement model of equation (1) and set A = D[1] as the identity and B(σ) = B 3D (σ) as three-dimensional symmetric Gaussian blur.
Simulated measurements were produced using the following forward model with σ true = 1 px and N denoting the multivariate normal distribution defined by the given mean vector and covariance matrix. To avoid partial volume effects, µ was upsampled by a factor of 8. This simulation permits to directly recover µ using one step Fourier inversion of the blur.
For the simulation study, we evaluated different numbers of photons I 0 ∈ {10 5 , 10 6 , ∞}. We defined µ according to the FORBILD head phantom (FORBILD phantoms 1999), which is commonly used in CT simulations, with 256 × 256 voxels and 32 slices at an (monochromatic) energy of 100 keV. We investigated the values of the TV regularization and the normalized sparsity measure for the restorations obtained using σ ∈ [0.8, 1.1] px.

Evaluation of blind deconvolution reconstruction for CT
For tomography, we used again the FORBILD head phantom at an (monochromatic) energy of 100 keV with 256 × 256 voxels and 32 slices. This corresponds to a voxelsize of 0.1 cm. The measurements were created according to where µ was up-sampled for data generation by a factor of 8 to avoid partial volume effects. In equation (12), A models a parallel beam geometry with 403 projections equidistantly distributed between 0 • to 180 • , fulfilling the Nyquist criterion. The flux was set to I 0 = 10 4 and the matrix B(σ) accounted for a 2D symmetric Gaussian detector blur parameterized by σ. The values were drawn from an uncorrelated multivariate normal distribution N with equal mean and variance to approximate photon statistics.
Three studies were conducted. The first study evaluates the differences between CT reconstructions using conventional TV regularization and the normalized sparsity measure by solving the respective objectives functions for µ at different values of the blur parameter σ and regularization strength β. For conventional TV regularization, the objective is given by equation (2) with σ = σ and ψ = 1 . The objective with the normalized sparsity measure is given by equation (5). For the simulation of the projections, we chose σ true = 1 px. Exhaustive reconstructions were performed for σ ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4} px. We denote β 0 to be the regularization strength such that the variance in a uniform region of the reconstruction with the true blur parameter matches the variance of the FBP in the same region. Consequently, the numerical values for the normalized sparsity measure are larger than for the conventional TV regularization to yield the same variance. Reconstructions were performed for β ∈ {1/2, 1, 2, 4}β 0 . A L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) routine initialized by filtered backprojection was used. To choose the number of iterations, the difference maps between the current reconstruction and the reconstruction obtained 100 iterations earlier using β = β 0 and σ = 1.0 px (both for the conventional TV regularization and the normalized sparsity measure) were analyzed. The RMSE values between the reconstructions obtained after 2000 iterations and the ones obtained after 1900 iterations were around 0.1 HU, which is orders of magnitude below the noise level, and the reconstructions were visually identical. Thus, the number of iterations was set to 2000 in the subsequent experiments.
The second study investigates blind deconvolution reconstruction for different values of σ true . Seven independent measurements were simulated according to equation (12) with σ true ∈ {0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6} px. For every simulation, seven reconstructions with σ ∈ σ true ± {0.0, 0.1, 0.2, 0.3} px were performed. As the noise characteristics changed with the different simulation studies, we scaled β = β 0 , which we used for reconstruction for σ true = 1.0 px, with the relative difference in variance of other measurements, to approximately achieve similar variances in the reconstructed volumes. The L-BFGS algorithm was used (initialized with a filtered backprojection) to solve the objective given by equation (5). The number of iterations was chosen by the same procedure as defined above. With a RMSE of below 1 HU (both for the conventional TV regularization and the normalized sparsity measure) and visually identical reconstructions at iteration 1000 and iteration 900, we concluded that 1000 iterations are sufficient for this study.

Physical experiments
The sample used for the experimental validation was a human wrist phantom. The data was acquired on a cone-beam CT test bench. The geometric parameters were given by a source-detector distance of 1184 mm and a source-axis distance of 595 mm. The flat-panel detector (4030CB, Varian, Palo Alto CA) has a pixelsize of 0.388 mm. The x-ray tube (Rad-94, Varian, Salt Lake City UT) was operated at 100 kVp with 100 mA and an exposure time of 6.3 ms. The average number of detected photons in a sample free region was 2.4 · 10 5 ph px −1 . The acquired projections were gain and offset corrected. The resulting 360 projections were spatially binned by a factor of two, to 340 × 100 px. The dimensions of the reconstructed volume were 340 2 × 120 voxels.
For blind deconvolution reconstruction, we again solved the respective objective function (5), using algorithm 1 with the L-BFGS routine. The RMSE between the reconstruction obtained after iteration 390 and the reconstruction obtained after iteration 400 was below 0.0001 cm −1 and both reconstructions were visually identical. We thus concluded that 400 iterations are sufficient for the experimental study. The regularization strength was set to β = 0.01 according visual inspection. To demonstrate the dependency of the objective function and the reconstructed slices on the blur parameter, reconstructions for all σ ∈ {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0} px were performed.

Evaluation of the regularization strategies in a simple 3D blur example
The conventional TV regularization and the normalized sparsity measure were evaluated for different noise parameters and numbers of photons in the simple 3D blur example.
In figures 2(a)-(c) the values of the conventional TV regularization and the normalized sparsity measure are shown for restorations obtained with different blur parameters. The values attributed to the conventional TV regularization are shown in blue on the left axis and the corresponding values for the normalized sparsity measure in orange on the right. The vertical lines at σ = 1.0 px mark the true blur parameter used for simulation. The results of the low flux scenario with I 0 = 10 5 are shown in figure 2(a). The dependency of both regularization approaches is similar, increasing with higher blur parameters without any distinct optima at the true blur parameter. Due to the high noise, a low value of σ would be estimated. In figure 2(b) the flux was increased by one order of magnitude. Although the behavior of the conventional TV regularization remains the same, the normalized sparsity measure yields a global optimum at the true blur parameter used for the simulation. Thus, the normalized sparsity measure correctly captures the prior assumptions imposed onto the sample. The results without noise are shown in figure 2(c). Overall, the conventional TV regularization still increases with increasing blur parameters. The optimum at the true blur parameter of the normalized sparsity measure is prominent. The asymmetric nature of the penalty can be seen by the negative curvature for underestimating the blur parameter and the positive curvature for an overestimation.
In figure 2(d) a zoomed image of the ground truth phantom used for simulation is depicted. Three reconstructions obtained at I = 10 6 for three different blur parameters (indicated by the vertical dashed lines in figure 2(b)) are shown in figures 2(e)-(g). For the reconstruction in figure 2(e), the blur parameter was underestimated. The noise level is reduced and the structures in the reconstruction are blurred. Matching the blur parameter for reconstruction with the one in the simulations yields a sharp image with moderate noise levels as shown in figure 2(f). An overestimation of the blur parameter as depicted in figure 2(g) gives overshoots at the edges of the structures and amplified noise levels. Finally, line profiles along the horizontally averaged regions marked by the colored rectangles in figures 2(d)-(g) emphasize the effect of the edge sharpness on the blur parameter. For the correct blur parameter the ground truth is accurately recovered, while an underestimation of the blur parameters blurs the structures and overshoots are visible for a corresponding overestimation.
We conclude that the conventional TV regularization results in a bias of the blur estimate towards smaller values. With increasing blur parameter the noise level in the restorations increases and thus the value of the conventional TV regularization. If the flux is sufficiently high, the normalized sparsity measure correctly captures the prior information about the sample.

Evaluation of blind deconvolution reconstruction for CT
In the first CT study, reconstructions with the conventional TV regularization and the normalized sparsity measure were performed for different values of the blur parameter and the regularization strength.
In figures 3(a)-(c), different parts of the objective function for the conventional TV regularization are shown, while in figures 3(d)-(f) the same parts of the objective with the normalized sparsity measure are shown. From left to right, the values of the data-fidelity, the regularization term, and the whole objective function are plotted over the blur parameter for the different regularization strengths, respectively. The location of the true blur parameter is indicated by a vertical solid gray line. To make the numerical values between different regularization strengths comparable, they are normalized to the first data point at σ = 0.5 px.
For conventional TV regularization, the data-fidelity, the regularization term, and thus also the objective function value increase monotonically with an increasing blur parameter. This results from the fact that with increasing blur parameter, the reconstructed noise level also increases. The change in the objective function values increases more rapidly if the blur parameter is overestimated than if it is underestimated. The monotonicity of the objective function implies that a joint optimization of both the reconstructed volume and the underlining blur parameter, is not possible as the 'no-blur' solution would be recovered.
For the reconstruction using the normalized sparsity measure, the data-fidelity still increases overall with the blur parameter. However, for an underestimation of the true blur parameter, the values are relatively constant. Around the true blur parameter a local minimum can be observed. For an increasing blur parameter, the datafidelity increases as the noise level rises. In contrast to the data-fidelity, the normalized sparsity measure yields a global optimum around the true value of the blur parameter, which matches the observation of the simple 3D blur model. As in the simple 3D example, the curvature for an underestimation of the blur parameter is negative and positive elsewhere. The sparsity of the object is captured for all the different noise levels. However, we observe a slight overestimation of the blur parameter. We believe that the overestimation can be explained by the additional blur components, which we do not explicitly model, such as algorithmic components like interpolation in the projection operations. Due to the comparably small dependency of the data-fidelity on the blur parameter compared to the regularization, the objective function value gives a distinct optimum around the true blur parameter. However, with increasing regularization strength, the blur parameter is slightly overestimated. This is  . First simulation study. In (a)-(c) the different parts of the objective function using conventional TV regularization are depicted over the blur parameter for different regularization strengths. In (a), the value of the data-fidelity is shown, in (b) the value of the regularization term and (c) shows the value of the objective function. In (d)-(f) the same information is shown for the proposed objective function using the normalized sparsity measure. In (g)-(k) reconstructions for different blur parameters are shown. The left sides of each figure show the reconstructions obtained with conventional TV regularization, while the right halves have been reconstructed using the normalized sparsity measure. In addition, the upper halves are shown in a more narrow grayvalue window than the lower halves.
in accordance with the assumption that the normalized sparsity measure also accounts for the blurring induced by the regularization, which becomes more prominent with increasing regularization strength. In general, the optimum becomes more distinct with increasing regularization strength as the prior knowledge about the sample becomes more important. For the weakest regularization strength β = 13, the two reconstructions with the largest blur parameter σ ∈ {1.3, 1.4} diverge due to noise amplifications and can therefore not be depicted.
In figures 3(g)-(k) reconstructions for the two approaches with a regularization strength of β = β 0 are depicted for different blur parameters σ ∈ {0.6, 0.8, 1.0, 1.2, 1.4} px. The left halves of the images were obtained using conventional TV regularization, while the right halves show reconstructions obtained with the normalized sparsity measure. In addition, the upper halves are shown with a more narrow window, while the lower halves utilize the full range of the phantom. Due to the choice of the regularization strength, the noise level for σ = 1.0 px is matched for both approaches. For low blur parameters, the noise level of the approach using the normalized sparsity measure is lower compared to the conventional TV regularization, which can be seen from the narrow windows. Respectively, higher noise levels are present for an overestimation of the blur parameter as the effective regularization strength is reduced by the normalization of the sparsity measure. The edge sharpness can best be analyzed from the lower halves. An underestimation blurs the edges of the structures, while an overestimation yields overshoots at the edges.
We find that the conventional TV regularization is not suited for recovering the blur parameter. Using the normalized sparsity measure we have a global optimum around the true blur parameter. This study validates the feasibility of blind deconvolution reconstruction, which estimates the optimal blur parameter jointly during reconstruction.
The second study examines blind deconvolution reconstructions of seven independent measurements simulated with varying σ true . In figure 4, the different parts of the objective functions are shown for the different measurements over the values of the blur parameter. To make the values between the different simulations comparable, we again normalized the values to the first data-point of the respective curves. Globally, the data-fidelity increases with increasing blur parameters. However, for larger blur parameters a minimum can be observed for the true blur parameter. As for the previous cases, the sparsity measure used for regularization has an optimum which slightly overestimates the true blur parameter. Again, this can be potentially explained by additional blur components, such as the interpolation that is inherent to projection operations or the fact that the non-binary edges of the phantom are being compensated. The overestimation becomes slightly less prominent with an increasing blur parameter, as the effect of the blur components that are not modelled becomes less significant. Finally, the objective function yields an optimum for the true blur parameter for all simulated scenarios, which further validates that the objective can be used for blind deconvolution reconstruction to estimate the underlying model parameter of the blur.
Finally, the last simulation study evaluates the proposed efficient search routine for blind deconvolution reconstruction. The results for two different initial guesses of the blur parameter are seen in figure 5. In gray, the objective function L γ (σ) is plotted after 100, 200, 300, 400, 500 (dashed) and 2000 (solid) iterations over the blur  Colored dots denote reconstructions that have been calculated by the optimization algorithm. Small dots refer to omitted reconstructions, while the larger dots denote the best reconstructions in each iteration. Dashed colored lines illustrate, which estimates were used for the subsequent reconstructions. In (a) σ (0) = 0.7 px was used as an initial guess for the blur parameter. In (b), the respective initial guess was set to σ (0) = 1.3 px.
parameters. In practise, these points would not be calculated. Points that are explicitly calculated are marked in blue and orange, respectively.
In figure 5(a) we used σ (0) = 0.7 px as initialization. Consequently, after N = 100 (m = 0) iterations, intermediate reconstructions for σ ∈ {0.6, 0.7, 0.8} px are obtained. The values coincide with those from the reference reconstructions depicted in gray. As the objective function value for σ (1) = 0.8 px is the lowest, the corresponding volume is used as initialization for further reconstruction at σ ∈ {0.7, 0.8, 0.9} px. The values for the reconstructions do not match the reference anymore, as a different initialization is used and for the case of σ = 0.8 px a slight mismatch is observed due to the fact that the optimization algorithm has been restarted. After an additional iteration m = 2, the true blur parameter σ (3) = σ true is already obtained. Additional iterations do not change the blur parameter, but refine the reconstructed volume further.
The course of the search routine for an overestimated initial guess of σ (0) = 1.3 px is depicted in figure 5(b). After four iterations m = 3, the true blur parameter σ (4) = σ true is obtained and the reconstruction is further refined in additional iterations.
We find that the true blur parameter is recovered in both cases after only a few iterations. The increase in computational time is moderate, because all reconstructions at the same iteration m can be performed independently and in parallel. Given the steeper slope of the objective function value, it might be advantageous to overestimate the initial guess for the value of the blur parameter.

Physical experiments
The results for the experimental validation of blind deconvolution reconstruction obtained by applying algorithm 1 for fixed values of σ are shown in figure 6.
In figures 6(a)-(c), the different parts of the objective function are plotted over the blur parameter again. The data fidelity increases monotonously with increasing blur parameter. This trend has already been observed in our second simulation study as depicted in figure 4(a), if the underlying blur parameter is small. The values of the intermediate iterations approach those of the converged solution relatively fast. The values of the nor malized sparsity measure have a global optimum around σ = 0.8 px. With increasing iterations the optimum of the regularization drifts slightly to higher values until convergence. In later iterations the improvements of the regularization are more prominent than for the data fidelity. The objective function has an optimum at σ = 0.6 px. Due to the regularization, the estimated blur parameter of subsequent iterations is shifted to larger values. We have already examined the evolution of the objective function in our third simulation study in figure 5. There, we observed the same shift of the optimum with increasing iterations.
For reference, in figure 6(d) a slice of the reconstructed human wrist phantom is depicted without the blur model using σ = 0.0 px. In figures 6(e)-(i) the corresponding reconstructed slices are shown for σ ∈ {0.4, 0.5, 0.6, 0.7, 0.8} px. The slices are cropped to fit the entire wrist. Additional zooms are provided. In blue, zooms on the edges of the wrist with a narrow colorbar are shown to better study the quality of the edges. For small blur parameters σ ∈ {0.0, 0.4, 0.5} px these edges are blurred out. In contrast, larger blur parameters σ ∈ {0.7, 0.8} px lead to overshooting artifacts. For a physically accurate depiction of the data, we would prefer the sharpest edges without any additional overshoots. This can be clearly seen for σ = 0.6 px. The zooms marked in orange highlight lower contrast features in the soft-tissues with a narrow gray-scale window. As the regularization strength was not optimized for the soft-tissue window, the improvements are less prominent in this case. However, the best separation between the low-contrast soft-tissue structures and the bone without having overshoots can also be observed for σ = 0.6 px. These results coincide nicely with the optimum of the objective function in convergence. Thus, we may jointly estimate this parameter during reconstruction with the proposed search routine for blind deconvolution reconstruction. In contrast to the simulation study, we showed that our approach works even with a mismatched blur model, since the Gaussian blur assumption does not perfectly match the physical system.

Conclusion
We have developed an objective function for CT with a distinct optimum at the desired blur parameters, enabling the simultaneous recovery of these parameters together with the deblurred reconstructed volume. Specifically, we used a normalized sparsity measure, which correctly captures the prior information imposed onto the sample. We demonstrated an efficient search routine to reconstruct the volume while jointly estimating underlying lowdimensional blur parameters. We compared our approach to conventional TV regularization, where we find that the objective function always favors a 'no-blur' solution. An extensive simulation study was performed to evaluate the performance of our approach for different regularization strengths and system blurs. We further validated our search routine for the joint estimation of the blur parameter. Finally, we applied our findings to experimental test bench data demonstrating the feasibility of this approach with an imperfect blur model.
As we have not established a convergence proof, which can be difficult for non-convex objectives, we instead have investigated empirically that the solutions are relatively robust. Therefore, we used the simulation study which was summarized in figure 3 with σ = 1.0 px and β = β 0 (using the normalized sparsity measure) and found that different initializations result in visually identical solutions with mutual RMSEs of around 1-10 HU (which is below the noise level). In particular, starting with high-frequency noise as initialization (which has a high normalization factor) or a blurred structure (which has a low normalization factor) yield solutions that are visually identical. Future work will include investigation of the convergence properties, the influence of the adaptive regularization strength on convergence speed and specialized optimization algorithms tailored to the exact forward model and objective function. For example, an optimization algorithm for the physical mean model with known blur model has been developed in Tilley et al (2018).
Better blur models with more parameters are possible. However, an extension of this approach to highdimensional blur models is not trivial. Although the objective function still gives a distinct optimum for the desired blur model parameters, the gradient-free optimization scheme is challenging, due to the exponential increase of the computational cost with the number of blur model parameters. For example, in ongoing work, low-dimensional blur models, e.g. a Lorentzian + Gaussian model, are investigated. The blur parameters within the search routine may be updated in an alternating fashion, again requiring only three reconstructions in parallel. In this work, we assumed shift-invariant blur models. This is reasonable for detector blur, for cone-beam geometries with large objects, source blur becomes depth-dependent. It has been demonstrated that employing shift-variant blur models in MBIR can improve image quality further (Tilley et al 2016b), and this work has possible extensions to such models.
Inherent to deconvolution problems, high noise levels or large system blur will eventually render this approach unstable. As seen in the evaluation of the regularization techniques, high noise levels corrupt the influence of the sample edges in such a way that a distinct optimum cannot be found. The limitations of jointly estimating the system blur during tomographic reconstruction will be investigated in future work. In particular, the impact of noise in physical systems with realistic detectors and appropriate fluence levels should be studied. Possible attempts to amend the corruption of the sample edges would be modified sparsity measures, which account differently for small differences in neighboring voxel values based on for instance the Huber penalty (Huber 2011).
In blind deconvolution of photographs as introduced in Krishnan et al (2011), blur kernel estimation is often performed on the high-frequencies of the image, whereas our data fidelity term accounts for all frequencies. For sufficiently weak regularization strengths, the data-fidelity could possibly dominate the behavior of the objective function, which would not guarantee an optimum at the correct parameter for all cases, e.g. for small blur parameters. Given the proposed solution, the estimation of the blur parameters can solely be determined from the value of the regularization alone. However, a gradient-based joint optimization approach would then not be feasible. One possible solution to this problem is to modify the data-fidelity term to emphasize the high-frequency components. For instance, including covariance matrices as introduced in Tilley et al (2016a) act as a high-pass filter suppressing the energy contained in the low-and mid-frequencies, which we will examine in future work.
In summary, our approach of jointly estimating the blur parameters during tomographic reconstruction has some distinct benefits compared to the conventional approaches. Firstly, we omit the need for time-consuming prior characterization of the system blur parameters, if a parameterized blur model is available. This can prove useful, especially when the blur parameters change frequently, for instance if the x-ray source parameters or the geometry is changed. Secondly, in contrast to any prior characterization of the system blur, which acts on the projection data, our integrated approach estimates the blur directly from the reconstructed volume. Thereby, more subtle sources of system blur can be directly compensated. These include interpolation in the projection operations or smoothing due to the regularization. Investigating to what extent these additional components can be accounted for is subject to future work. Our approach is neither limited to the examined forward model nor to the established data-fidelity model. We believe that the introduction of the normalized sparsity measure for CT might also prove beneficial in different areas, in particular, where model parameters that change the noise characteristics have to be chosen carefully. One example are the material specific constants in propagation-based phase-contrast CT, which might be simultaneously estimated with this approach.