SABER: A Systems Approach to Blur Estimation and Reduction in X-ray Imaging

TL;DR

SABER introduces a systems-based method for modeling, estimating, and reducing blur in X-ray images by combining multiple PSFs and applying deblurring algorithms, improving image clarity and contrast.

Contribution

The paper presents a novel systems approach to model and estimate the combined source and detector PSFs in X-ray imaging, enabling effective blur reduction.

Findings

Accurate estimation of source and detector PSFs from multiple radiographs.

Effective blur reduction using estimated PSFs through deblurring algorithms.

Improved image sharpness and contrast in X-ray radiographs.

Abstract

Blur in X-ray radiographs not only reduces the sharpness of image edges but also reduces the overall contrast. The effective blur in a radiograph is the combined effect of blur from multiple sources such as the detector panel, X-ray source spot, and system motion. In this paper, we use a systems approach to model the point spread function (PSF) of the effective radiographic blur as the convolution of multiple PSFs, where each PSF models one of the various sources of blur. In particular, we model the combined contribution of X-ray source and detector blurs while assuming negligible contribution from other forms of blur. Then, we present a numerical optimization algorithm for estimating the source and detector PSFs from multiple radiographs acquired at different X-ray source to object (SOD) and object to detector distances (ODD). Finally, we computationally reduce blur in radiographs using deblurring algorithms that use the estimated PSFs from the previous step. Our approach to estimate and reduce blur is called SABER, which is an acronym for systems approach to blur estimation and reduction.

Figures

Tables

Table 1. TABLE I: FWHMs W s ​ x , W s ​ y subscript 𝑊 𝑠 𝑥 subscript 𝑊 𝑠 𝑦 W_{sx},W_{sy} of the X-ray source PSF estimated using one set of horizontal and vertical edge radiographs.

horizontal	vertical	$W_{sx}$	$W_{sy}$
SOD ($mm$)	SOD ($mm$)	($\mu m$)	($\mu m$)
12.0	13.0	3.2	3.6
24.8	24.8	3.6	3.9
37.5	37.5	4.5	4.6
50.3	50.3	6.9	6.9
65.3	60.0	13.0	25.9

Table 2. TABLE II: Parameters q 𝑞 q , W d ​ 1 subscript 𝑊 𝑑 1 W_{d1} , and W d ​ 2 subscript 𝑊 𝑑 2 W_{d2} of the detector PSF estimated using one set of horizontal and vertical edge radiographs.

horizontal	vertical	$W_{d1}$	$W_{d2}$	$q$
SOD ($mm$)	SOD ($mm$)	($\mu m$)	($\mu m$)
12.0	13.0	12.6	130.6	0.89
24.8	24.8	5.7	134.9	0.91
37.5	37.5	3.4	158.3	0.92
50.3	50.3	2.3	127.0	0.92
65.3	60.0	1.9	118.3	0.93

Table 3. TABLE III: X-ray source PSF and detector PSF parameters estimated using two sets of horizontal and vertical edge radiographs.

horizontal	vertical	$W_{sx}$	$W_{sy}$	$W_{d1}$	$W_{d2}$	$q$
SOD ($mm$)	SOD ($mm$)	($\mu m$)	($\mu m$)	($\mu m$)	($\mu m$)
12.0,24.8	13.0,24.8	2.6	3.0	1.5	135.7	0.90
12.0,37.5	13.0,37.5	2.6	3.0	1.9	150.1	0.91
12.0,50.3	13.0,50.3	2.6	3.0	1.9	142.7	0.91
12.0,65.3	13.0,60.0	2.7	3.1	1.8	148.2	0.92
24.8,37.5	24.8,37.5	2.7	3.0	1.9	138.0	0.91
24.8,50.3	24.8,50.3	2.7	3.0	1.9	126.7	0.92
24.8,65.3	24.8,60.0	2.7	3.1	1.8	126.8	0.92
37.5,50.3	37.5,50.3	2.9	2.9	1.9	134.7	0.92
37.5,65.3	37.5,60.0	2.9	3.1	1.9	134.9	0.93
50.3,65.3	50.3,60.0	2.9	2.9	1.9	119.7	0.93

Table 4. TABLE IV: Mean and standard deviation of estimated parameters shown in Table III .

	$W_{sx}$	$W_{sy}$	$W_{d1}$	$W_{d2}$	$q$
Mean	2.7$\mu m$	3.0$\mu m$	1.8$\mu m$	135.7$\mu m$	0.92
Std Dev	0.12$\mu m$	0.07$\mu m$	0.11$\mu m$	9.2$\mu m$	0.01
$100\frac{\text{Std Dev}}{\text{Mean}}$	4.4%	2.3%	6.1%	6.8%	1.09%

Table 5. TABLE V: Estimated parameters before and after step (3) of algorithm 1 using radiographs at SODs of 37.5 ​ m ​ m 37.5 𝑚 𝑚 37.5mm and 50.3 ​ m ​ m 50.3 𝑚 𝑚 50.3mm .

	$W_{sx}$	$W_{sy}$	$W_{d1}$	$W_{d2}$	$q$
	($\mu m$)	($\mu m$)	($\mu m$)	($\mu m$)
Before step (3)	4.5	4.6	2.3	67.6	0.92
After step (3)	2.9	2.9	1.9	134.7	0.92

Table 6. TABLE VI: Performance comparison between various choices of blur model. Performance evaluation is on radiographs at SODs different from those used during blur estimation. Root mean squared error (RMSE) between the measured radiographs I k ​ ( i , j ) subscript 𝐼 𝑘 𝑖 𝑗 I_{k}(i,j) and its predictions I ¯ k ​ ( i , j ) subscript ¯ 𝐼 𝑘 𝑖 𝑗 \bar{I}_{k}(i,j) is used to compare blur prediction performance. RMSE between the deblurred radiographs and the transmission functions T k ​ ( i , j ) subscript 𝑇 𝑘 𝑖 𝑗 T_{k}(i,j) is used to compare the deblurring performance.

			Prediction ${\overline{I}}_k(i,j)$		Deblur, SD = $0.01$		Deblur , SD = $0.006$
PSF	Source	Detector	$24.8mm$/	$65.3mm$/	$24.8mm$/	$65.3mm$/	$24.8mm$/	$65.3mm$/
Model	Blur	Blur	$24.8mm$	$60mm$	$24.8mm$	$60mm$	$24.8mm$	$60mm$
Gaussian	Yes	No	0.0172	0.0157	0.0717	0.0192	0.0516	0.0172
Exponential	Yes	No	0.0164	0.0157	0.0387	0.0188	0.0338	0.0171
Gaussian	No	Yes	0.0228	0.0086	0.0309	0.0327	0.0296	0.0226
Exponential	No	Yes	0.0223	0.0083	0.0322	0.0161	0.0309	0.0132
Gaussian	Yes	Yes	0.0149	0.0081	0.0367	0.0168	0.0289	0.0126
Exponential	Yes	Yes	0.0147	0.0078	0.0264	0.0139	0.0243	0.0111

Table 7. TABLE VII: Variation in deblurring performance due to SOD dependent variation in estimated PSF parameters (Table III ). The last two columns show the RMSE between the transmission functions T k ​ ( i , j ) subscript 𝑇 𝑘 𝑖 𝑗 T_{k}(i,j) and the deblurred radiographs at SOD of 24.8 ​ m ​ m 24.8 𝑚 𝑚 24.8mm .

Horizontal	Vertical	Deblur	Deblur
SOD ($mm$)	SOD ($mm$)	SD = $0.01$	SD = $0.006$
12.0,37.5	13.0,37.5	0.0242	0.0222
12.0,50.3	13.0,50.3	0.0246	0.0226
12.0,65.3	13.0,60.0	0.0241	0.0220
37.5,50.3	37.5,50.3	0.0264	0.0243
37.5,65.3	37.5,60.0	0.0258	0.0236
50.3,65.3	50.3,60.0	0.0262	0.0242

Equations

$$\frac{I_s(i,j)}{I_b(i,j)}={\int }_{E}S(E)\exp \left(-{\int }_{L(i,j)}\mu \left(E,r\right)dr\right)dE,$$

$$\tilde{T}(i,j)={\int }_{E}S(E)\exp \left(-{\int }_{L(i,j)}\mu \left(E,r\right)dr\right)dE.$$

$$I(i,j)=\frac{I_s(i,j)-I_d(i,j)}{I_b(i,j)-I_d(i,j)}=\tilde{T}(i,j)+w(i,j),$$

$$I_k(i,j)=T_k(i,j)\ast p_k^{(s)}(i,j)\ast p^{(d)}(i,j)\ast p_k^{(m)}(i,j),$$

$$T_k(i,j)=l_k+(h_k-l_k){\tilde{T}}_k(i,j),$$

$$I_k(i,j)=l_k+(h_k-l_k){\tilde{T}}_k(i,j),\text{ where }w(i,j)=1.$$

$$p^{(s)}\left(i,j\right)=\frac{1}{Z_s}\exp \left(-{\Delta }^r{\left(i^2{s}_{sx}^2+j^2{s}_{sy}^2\right)}^{\frac{r}{2}}\right),$$

$$Z_s=\sum _i\sum _j\exp \left(-{\Delta }^r{\left(i^2{s}_{sx}^2+j^2{s}_{sy}^2\right)}^{\frac{r}{2}}\right).$$

$$p_k^{(s)}(i,j)=\frac{1}{Z_{s,k}}\exp \left(-{\Delta }^r\frac{SO{D}_k^r}{OD{D}_k^r}{\left(i^2{s}_{sx}^2+j^2{s}_{sy}^2\right)}^{\frac{r}{2}}\right),$$

$$Z_{s,k}=\sum _i\sum _j\exp \left(-{\Delta }^r\frac{SO{D}_k^r}{OD{D}_k^r}{\left(i^2{s}_{sx}^2+j^2{s}_{sy}^2\right)}^{\frac{r}{2}}\right).$$

$$\begin{gathered} p^{(d)}(i,j)=q\frac{1}{Z_{d1}}\exp \left(-s_{d1}^r{\Delta }^r{\left(i^2+j^2\right)}^{\frac{r}{2}}\right) \\ +(1-q)\frac{1}{Z_{d2}}\exp \left(-s_{d2}^r{\Delta }^r{\left(i^2+j^2\right)}^{\frac{r}{2}}\right) \end{gathered}$$

$$Z_{d1}$$

$$Z_{d2}$$

$$p_k^{(m)}(i,j)=\{\begin{array}{ll}1,& \text{ if }i=0,j=0\\ 0,& \text{ otherwise}\end{array}.$$

$$\left(\begin{array}{c}{\hat{s}}_{sx},{\hat{s}}_{sy},{\hat{s}}_{d1},{\hat{s}}_{d2},\hat{q},\\ {\hat{l}}_1,\cdots ,{\hat{l}}_K,{\hat{h}}_1,\cdots ,{\hat{h}}_K\end{array}\right)=\underset{\begin{array}{c}{s}_{sx},s_{sy},s_{d1},s_{d2},q\\ l_1,\cdots ,l_K,h_1,\cdots ,h_K\end{array}}{\mathrm{argmin}}\sum _{k=1}^K{E}_k$$

$$\begin{gathered} \text{where }{E}_k=0.5\sum _{i,j}{w}_k(i,j)(I_k(i,j)- \\ {T_k(i,j)\ast p_k^{(s)}(i,j)\ast p^{(d)}(i,j)\ast p_k^{(m)}(i,j))}^2, \\ \text{subject to }{s}_{sx}\ge 0,s_{sy}\ge 0,s_{d1}\ge 0,s_{d2}\ge 0, \\ {q}_{low}\le q\le 1,-1\le l_k\le 0.5,0.5\le h_k\le 2, \end{gathered}$$

$$\left({\hat{s}}_{sx},{\hat{s}}_{sy}\right)=\underset{s_{sx}\ge 0,s_{sy}\ge 0}{\mathrm{argmin}}\sum _{k\in {\Omega }_s}{E}_k.$$

$$\left({\hat{s}}_{d1},{\hat{s}}_{d2},\hat{q}\right)=\underset{s_{d1}\ge 0,s_{d2}\ge 0,q_{low}\le q\le 1}{\mathrm{argmin}}\sum _{k\in {\Omega }_d}{E}_k$$

$$\begin{gathered} \left(\begin{array}{c}{\hat{s}}_{sx},{\hat{s}}_{sy},{\hat{s}}_{d1},{\hat{s}}_{d2},\hat{q},\\ {\hat{l}}_1,\cdots ,{\hat{l}}_K,{\hat{h}}_1,\cdots ,{\hat{h}}_K\end{array}\right)=\underset{\begin{array}{c}{s}_{sx},s_{sy},s_{d1},s_{d2},q\\ l_1,\cdots ,l_K,h_1,\cdots ,h_K\end{array}}{\mathrm{argmin}}\sum _{k=1}^K{E}_k \\ \text{subject to the constraints in equation eqrefeq:optconstr}. \end{gathered}$$

$${\hat{p}}_k\left(i,j\right)={\hat{p}}_k^{(s)}(i,j)\ast {\hat{p}}^{(d)}(i,j)\ast {\hat{p}}_k^{(m)}(i,j),$$

$$\begin{gathered} {\hat{T}}_k(i,j)= \\ \underset{T_k(i,j),\forall i,j}{\mathrm{argmin}}\{\sum _{i,j}{w}_k(i,j){\left(I_k(i,j)-T_k(i,j)\ast {\hat{p}}_k\left(i,j\right)\right)}^2\text{vspace}-0.2in \\ +\beta \sum _{((i,j),(m,n))\in \mathcal{N}}\tilde{w}(i,j,m,n)|T_k(i,j)-T_k(m,n){|}^{1.2}\}, \end{gathered}$$

$$p^{(d)}(i,j)=\{\begin{array}{ll}1,& \text{ if }i=0,j=0\\ 0,& \text{ otherwise}\end{array}.$$

$$p_k^{(s)}(i,j)=\{\begin{array}{ll}1,& \text{ if }i=0,j=0\\ 0,& \text{ otherwise}\end{array}.$$

$$\frac{\partial f}{\partial v}=\sum _k\sum _{i,j}{w}_k(i,j)\left(I_k(i,j)-{\bar{I}}_k(i,j)\right)\frac{\partial {\bar{I}}_k(i,j)}{\partial v}$$

$$\begin{gathered} \frac{\partial {\bar{I}}_k(i,j)}{\partial v}=T_k(i,j)\ast \frac{1}{Z_{s,k}^2}(Z_{s,k}\frac{\partial g_{s,k}(i,j)}{\partial v}- \\ {g}_{s,k}(i,j)\frac{\partial Z_{s,k}}{\partial v})\ast p^{(d)}(i,j)\ast p_k^{(m)}(i,j), \end{gathered}$$

$$\frac{\partial g_{s,k}(i,j)}{\partial v}=\{\begin{array}{l}-r{g}_{s,k}(i,j){\Delta }^r\frac{SO{D}_k^r}{OD{D}_k^r}{s}_{sx}{\left(s_{sx}^2{i}^2+s_{sy}^2{j}^2\right)}^{\frac{r}{2}-1}{i}^2,\\ \text{ if v is ssx},\\ -r{g}_{s,k}(i,j){\Delta }^r\frac{SO{D}_k^r}{OD{D}_k^r}{s}_{sy}{\left(s_{sx}^2{i}^2+s_{sy}^2{j}^2\right)}^{\frac{r}{2}-1}{j}^2,\\ \text{ if v is ssy}.\end{array}$$

$$\frac{\partial {\bar{I}}_k(i,j)}{\partial v}=T_k(i,j)\ast p_k^{(s)}(i,j)\ast \frac{\partial p^{(d)}(i,j)}{\partial v}\ast p_k^{(m)}(i,j),$$

$$\frac{\partial p^{(d)}(i,j)}{\partial v}=\{\begin{array}{ll}\frac{q}{Z_{d1}^2}\left(Z_{d1}\frac{\partial g_{d1}(i,j)}{\partial v}-g_{d1}(i,j)\frac{\partial Z_{d1}}{\partial v}\right),& \text{if v is sd1},\\ \frac{1-q}{Z_{d2}^2}\left(Z_{d2}\frac{\partial g_{d2}(i,j)}{\partial v}-g_{d2}(i,j)\frac{\partial Z_{d2}}{\partial v}\right),& \text{if v is sd2},\\ \frac{1}{Z_{d1}}{g}_{d1}(i,j)-\frac{1}{Z_{d2}}{g}_{d2}(i,j),& \text{if v is q},\end{array}$$

$$\frac{\partial {\bar{I}}_k(i,j)}{\partial v}=\frac{\partial T_k(i,j)}{\partial v}\ast p_k^{(s)}(i,j)\ast p^{(d)}(i,j)\ast p_k^{(m)}(i,j),$$

$$\text{where }\frac{\partial T_k(i,j)}{\partial v}=\{\begin{array}{ll}1-{\tilde{T}}_k(i,j),& \text{if v is lk},\\ {\tilde{T}}_k(i,j),& \text{if v is hk}.\end{array}$$

Full text

I Introduction

X-ray imaging systems are widely used for 2D and 3D non-destructive characterization and visualization of a wide range of objects. The ability of X-rays to penetrate deep inside a material makes it a useful tool to visualize the interior morphology of objects. X-ray imaging is very popular in applications such as industrial imaging [1, 2, 3, 4], medical diagnosis [5, 6], and border security [4, 7, 8, 9]. A schematic representation of a typical cone-beam X-ray imaging system is shown in Fig. 1. Here, an object is exposed to a diverging beam of X-ray radiation that is generated by a X-ray source. The X-rays get attenuated as they propagate through the object and an image of the attenuated X-ray intensity is recorded by a detector consisting of a flat-panel 2D array of sensors.

The resolution of an X-ray image, also called a radiograph, is limited by several factors. The detector used to record the radiograph imposes a fundamental limit on the resolution by fixing the smallest pixel size. Also, detector cross-talk where energy from one sensor pixel leaks into its neighboring sensor pixels [10] causes blur and contrast reduction. Furthermore, the finite non-zero area of the X-ray source aperture manifests as additional blur in the radiograph [11, 12]. Other causes of blur include motion or vibrations in the sample stage or the imaging equipment. Blur is a major detriment in dimensional metrology applications where it is critical to accurately estimate the relative physical distances between image features. In medical applications, blur hinders the ability to resolve small features such as tumors that may only be a few pixels in size [13]. Hence, mitigating blur is vital to many applications especially when quantitative accuracy is important.

In order to quantify radiographic blur, we use a parametric mathematical model to describe the blurring process. Several papers address the problem of mathematically modeling numerous types of blur such as X-ray source spot blur, detector cross-talk, motion blur, and scatter [14, 15, 11]. A popular strategy is to model blur as a convolution operation with a certain point spread function (PSF). Here, the PSF used to model blur is derived using either simulation [16] or data-driven approaches [17, 18, 19, 20, 21, 22, 23]. Simulation of PSF relies on precise engineering knowledge of the relevant imaging equipment that may not always be available. Alternatively, data-driven approaches estimate PSF directly from radiographs of known well characterized objects.

Data-driven approaches to blur estimation relies on calculation of the PSF from radiographs of an object with known composition and shape such as a rollbar, slit, pinhole, or other test objects [17, 18, 19, 20, 21]. However, these methods only estimate one PSF that either models the X-ray source blur, the detector blur, or the total effective blur. They do not address the problem of disentangling and estimating the PSFs of all the different types of blur that simultaneously affect a radiograph. The magnitude of each blur varies depending on the relative placement of the X-ray source, object, and detector. The ability to estimate each individual blur PSF will allow us to predict the radiograph blur for any spatial configuration of X-ray source, object, and detector by appropriately recombining the estimated blur PSFs.

Blur in X-ray radiographs can be reduced either by hardware upgrades or using computational algorithms to reduce blur after the experiment. The former approach may not always be feasible due to cost or physical constraints. Alternatively, deblurring algorithms are a cheaper solution to computationally reduce blur in radiographs. The estimated blur PSFs can be used in a wide variety of deblurring algorithms [24, 25, 26, 27, 28, 29] to reduce blur in radiographs.

In this paper, we present an extensible framework called systems approach to blur estimation and reduction (SABER) for modeling, characterizing, and reducing the various types of blur in a X-ray imaging system. We apply a systems approach to model blur by expressing the PSF of the effective blur in a radiograph as the convolution of multiple blur PSFs with varying origins. Then, we present a numerical optimization algorithm to disentangle and estimate the individual blur PSFs from radiographs of a Tungsten plate rollbar. In particular, we focus on the simultaneous estimation of the X-ray source and detector PSFs while assuming negligible motion blur and scatter. Preliminary results using this approach was previously published as an extended abstract [12]. Software that implements SABER is freely available in the form of a open-source python package called pysaber that is documented at the link https://pysaber.readthedocs.io/.

In section II, we present the underlying principles and mathematics of X-ray imaging. We formulate our blur model in section III and estimate its parameters using an optimization algorithm in section IV. In section V, we present two approaches to deblurring radiographs. Finally, results using real experimental data from a Zeiss Xradia Versa X-ray imaging scanner are presented in section VI.

II Background

Beer’s law [30] is used to express the magnitude of X-ray attenuation within an object in terms of its thickness and material properties. X-ray attenuation is dependent on a material property called linear attenuation coefficient (LAC) which is a function of the object’s chemical composition, density, and X-ray energy. At detector pixel $(i,j)$ shown in Fig. 2, according to Beer’s law, the ratio of the X-ray measurement with the object, $I_{s}(i,j)$, and the X-ray measurement without the object, $I_{b}(i,j)$, is given by,

[TABLE]

where $\mu(E,r)$ is the LAC of the object at X-ray energy $E$ and spatial location $r$, $L(i,j)$ is the line along which $\mu(E,r)$ is integrated, and $S(E)$ is the X-ray spectral density such that $\int_{E}S(E)dE=1$ [4]. We will call the expression on the right side of the equality in equation (1) as the ideal transmission function $\tilde{T}(i,j)$, i.e.,

[TABLE]

The relation in equation (1) is valid in the absence of imaging non-idealities such as noise, blur, and scatter. Dark current is one such non-ideal characteristic where the detector measurements do not drop to zero when there is no X-ray radiation. To compensate for this effect, detector measurements are made without the X-ray beam and subtracted from measurements with the X-ray beam. Detector measurements are also corrupted by electronic noise and photon counting noise that are often modeled as additive Gaussian noise. Let $I(i,j)$ denote the normalized measurement defined as the ratio of the dark noise corrected measurements with and without the object. If $I_{d}(i,j)$ denotes the measurement at detector pixel $(i,j)$ without X-rays, then,

[TABLE]

where $w(i,j)$ is Gaussian noise. In the next section, we will modify equation (3) to account for the effects of more non-idealities such as X-ray source blur and detector blur.

The term “radiograph” will henceforth be used to refer to the normalized radiograph i.e., $I(i,j)$ over all pixel locations $(i,j)$. Also, “bright field” is used to refer to $I_{b}(i,j)$ and “dark field” is used to refer to $I_{d}(i,j)$.

III Formulation of Blur Model

The effect of various forms of blur on X-ray radiographs is modeled as a linear space-invariant phenomenon. The blur model expresses the radiograph $I(i,j)$ as the convolution of a transmission function $T(i,j)$ with multiple two dimensional PSFs each of which models a different form of blur. In the absence of non-idealities such as scatter or temporal drift in values of $I_{b}(i,j)$ and $I_{d}(i,j)$, the transmission function $T(i,j)$ is equal to the ideal transmission function $\tilde{T}(i,j)$ given in equation (2).

First, we will consider the effect of blur due to X-ray source, detector, and system motion [31]. The size of the source and motion PSFs are a function of the X-ray source to object distance (SOD) and object to detector distance (ODD). Since detector blur is due to cross-talk between detector pixels, it does not vary with SOD and ODD. Let $p^{(s)}_{k}(i,j)$ and $p^{(m)}_{k}(i,j)$ denote the PSFs of the X-ray source blur and motion blur respectively on the detector plane at a SOD of $SOD_{k}$ and ODD of $ODD_{k}$. If $I_{k}(i,j)$ is the radiograph and $T_{k}(i,j)$ is the transmission function at a SOD of $SOD_{k}$ and ODD of $ODD_{k}$, then,

[TABLE]

where $\ast$ denotes 2D convolution (Fig. 3). Based on the analysis in [32, 33, 34], we will ignore the effect of motion PSF since its full width half maximum (FWHM) was determined to be much smaller than our detector pixel width. Thus, our goal is to estimate the PSFs of X-ray source and detector blur given equation (4) and radiographs $I_{k}(i,j)$ at various SOD and ODD.

III-A Transmission Function

Simultaneous estimation of the transmission function, $T_{k}(i,j)$, and the blur PSFs, $p^{(s)}_{k}(i,j)$ and $p^{(d)}(i,j)$, is an ill-posed problem. Hence, we scan a object with known dimensions and chemical composition. Since our object is known, we can compute the ideal transmission function $\tilde{T}_{k}(i,j)$ for each scan (indexed by $k$) using equation (2). The transmission function $T_{k}(i,j)$ is then modeled as an affine transformation of the ideal transmission function $\tilde{T}_{k}(i,j)$ such that

[TABLE]

where $l_{k}$ and $h_{k}$ are scalars. The parameters $l_{k}$ and $h_{k}$ are used to compensate for uniform shifts in the measured X-ray intensity due to non-idealities such as drift in the values of $I_{b}(i,j)$ and $I_{d}(i,j)$ over time, inaccuracies in the calculation of $\tilde{T}_{k}(i,j)$, and low-frequency effects such as scatter. Temporal drifts in $I_{b}(i,j)$ and $I_{d}(i,j)$ can occur if there is any change in dark current or X-ray source intensity.

III-A1 Calibration Data

For our experiments, we fabricated a Tungsten plate with a sharp edge and uniform thickness that is sufficient to block all incoming X-rays (Fig. 4). The plate is then exposed to X-rays such that the sharp edge passes close to the center of the radiograph as shown in Fig. 5 (a). The plane of the Tungsten plate is aligned such that it is perpendicular to the direction of X-ray propagation. Since perfect alignment is very challenging, the sharp edges are rolled [19] to one meter radius of curvature, which permits alignment errors of up to one degree. Fig. 5 (b) shows the ideal transmission function for the radiograph in Fig. 5 (a). The procedure used to compute the ideal transmission function is presented in Appendix A.

III-A2 Initialization of $l_{k}$ and $h_{k}$

Due to scatter and variations in $I_{b}(i,j)$ and $I_{d}(i,j)$, $T_{k}(i,j)$ is not the same as $\tilde{T}_{k}(i,j)$. The values of $l_{k}$ and $h_{k}$ that determine the relation between $T_{k}(i,j)$ and $\tilde{T}_{k}(i,j)$ in equation (5) are jointly estimated during the blur estimation procedure that will be presented in the next section. However, to aid this procedure, we have to supply initial estimates for $l_{k},h_{k}$ that are determined by finding the best values for $l_{k},h_{k}$ that solve the following over-determined set of linear equations,

[TABLE]

Here, $w(i,j)$ is the weight term that is [math] for pixels that are close to the boundaries of the image $I_{k}(i,j)$ and $1$ elsewhere (Fig. (5) (c)). The padded region of $\tilde{T}(i,j)$ is also excluded from consideration during parameter estimation (Fig. (5) (d)). Equation (6) is solved using RANSAC regression [35, 36]. We used the scikit-learn [36] implementation of RANSAC111RANSAC Parameters: Minimum number of randomly chosen samples was $10$. For a data sample to be classified as an inlier, the maximum residual threshold was $0.1$. in this paper.

III-B Source Blur

The PSF of source blur is mathematically modeled using a 2D density function. This function is parameterized by two scale parameters $s_{sx}$ and $s_{sy}$ that are a measure of the spatial width of the PSF along the $x-$axis and $y-$axis in the plane of the X-ray source. If $\Delta$ denotes the width of each pixel, the PSF of source blur in the plane of the X-ray source is modeled as,

[TABLE]

where $r\geq 1$ and $Z_{s}$ is the normalizing constant given by,

[TABLE]

The constant $Z_{s}$ ensures that $p^{(s)}\left(i,j\right)$ sums to one when summed over all $(i,j)$. In equation (7), a Gaussian density function is obtained by setting $r=2$ and an exponential density function is obtained when $r=1$. The scale parameters are related to the full width half maximums (FWHM) along the $x-$axis as $W_{sx}=2\log(2)^{\frac{1}{r}}/s_{sx}$ and along the $y-$axis as $W_{sy}=2\log(2)^{\frac{1}{r}}/s_{sy}$. By definition, FWHM $W_{sx}$ is the distance between two points along the $x-$axis where the PSF drops to half of its maximum value and FWHM $W_{sy}$ is the corresponding distance along the $y-$axis.

At the detector plane, the FWHM of source blur is scaled by a factor of ODD/SOD. The change in FWHM of source blur with varying ODD and SOD is depicted in Fig. 6. For the $k^{th}$ radiograph, the PSF of source blur on the detector plane is given by,

[TABLE]

where $Z_{s,k}$ is the normalizing constant given by,

[TABLE]

III-C Detector Blur

The detector blur is modeled using a mixture of two density functions with scale parameters $s_{d1}$ and $s_{d2}$. If $p^{(d)}(i,j)$ denotes the PSF of detector blur, then,

[TABLE]

where $q$ is the mixture parameter and $Z_{d1},Z_{d2}$ are normalizing constants such that

[TABLE]

Also, $q$ and $(1-q)$ are parameters that function as weights for the two density functions. The scale parameters are related to the FWHMs of the two exponential functions in equation (11) as $W_{d1}=2\log(2)^{\frac{1}{r}}/s_{d1}$ and $W_{d2}=2\log(2)^{\frac{1}{r}}/s_{d2}$. A Gaussian mixture model for detector blur is obtained by setting $r=2$ and an exponential mixture model is obtained when $r=1$.

In our experiments, we expect the detector PSF to be dominated by the first exponential with a large weight of $q\approx 0.9$ and a very small FWHM $W_{d1}$ that spans only a few pixels. The second exponential is expected to have a smaller weight of $1-q\approx 0.1$ and a very large FWHM $W_{d2}$ that spans several hundreds of pixels. Depending on the X-ray imaging system, the values for $W_{d1},W_{d2}$ and $q$ may change but equation (11) is still expected to be a good model for detector PSF.

III-D Motion Blur

An additional source of blur is the relative motion of all components in the X-ray system. If these components move during the data acquisition period, then they will result in motion of the object’s image on the detector surface. Such motion will distribute the apparent location of an edge or feature over a range of locations, appearing to blur it in a manner that can be captured with a motion PSF that is convolved with the image. The system motion is calculated via a geometric analysis of the relative locations of components, propagated to the detector screen. This is done using an uncertainty propagation model presented in [32, 33, 34]. Motion blur is denoted by $p^{(m)}_{k}(i,j)$ and is typically modeled as a Gaussian density function. This form of blur also depends on the relative positions of X-ray source, object, and detector. However, we will effectively ignore the motion PSF by modeling it as a delta function since its FWHM was determined to be much smaller than the detector pixel width in our experiments. Thus, the form of $p^{(m)}_{k}(i,j)$ is given by,

[TABLE]

IV Estimation of Blur Model

Blur model estimation is the process of estimating all the parameters $s_{sx},s_{sy},s_{d1},s_{d2},$ and $q$ of the PSFs given known values for $I_{k}(i,j)$ and $\tilde{T}_{k}(i,j)$ in equations (4) and (5). We will only estimate $p^{(s)}_{k}(i,j)$ and $p^{(d)}(i,j)$ in equation (4) since $p^{(m)}_{k}(i,j)$ is assumed to be a constant delta function. By constraining $p^{(s)}_{k}(i,j)$ to take the shape of the density function in equation (7) and $p^{(d)}(i,j)$ to take the shape of the mixture density function in equation (11), blur estimation reduces to the problem of estimating the parameters $s_{sx}$, $s_{sy}$, $s_{d1}$, $s_{d2}$, and $q$. For every $k^{th}$ radiograph, the parameters $l_{k}$ and $h_{k}$ in equation (5) are also treated as unknowns and jointly estimated along with the PSF parameters.

Note that the form of source PSF as evaluated on the detector plane given by equation (9) depends on the ratio of the object to detector distance $ODD_{k}$ and the source to object distance $SOD_{k}$. Thus, while the amount of source blur in the radiograph $I_{k}(i,j)$ is a function of $ODD_{k}$ and $SOD_{k}$, the detector blur does not change with $ODD_{k}$ and $SOD_{k}$. To estimate both source and detector PSF parameters, it is necessary to acquire radiographs at a minimum of two different values of $ODD_{k}/SOD_{k}$. Also, since the source PSF has two scale parameters modeling the width along $x-$axis and $y-$axis, we need radiographs for cases when the Tungsten edge is horizontal and vertical.

We use numerical optimization to estimate all blur and transmission parameters given radiographs at different values of ODD and SOD. The parameters are estimated by solving,

[TABLE]

[TABLE]

where $q_{low}=0.8$ and $K$ is the total number of radiographs. The parameters $\hat{s}_{sx}$, $\hat{s}_{sy}$, $\hat{s}_{d1}$, $\hat{s}_{d2}$, $\hat{q}$, $\hat{l}_{k}$, and $\hat{h}_{k},\forall k$ are the estimated values of $s_{sx}$, $s_{sy}$, $s_{d1}$, $s_{d2}$, $q$, $l_{k}$, and $h_{k},\forall k$ respectively. The lower bound of $0.8$ for $q$ ensures that the first exponential in equation (11) is chosen to have the largest weight.

We solve the minimization problem in equation (15) using the L-BFGS-B algorithm [37]. L-BFGS-B is a quasi-Newton method for solving optimization problems and requires information of the gradient of the objective function, $\sum_{k=1}^{K}E_{k}$, with respect to all the variables being optimized. It also supports bound constraints on each of the variables that are optimized. However, since the optimization problem in equation (15) is non-convex, the optimized variables may be stuck in a local minima or a saddle point. Thus, good initialization of each variable is necessary to ensure reliable convergence to a solution that best fits the measured radiographs. Our approach to solving equation (15) is outlined in algorithm 1. We use the L-BFGS-B implementation in the python package scipy [38, 39]. The gradients that must be supplied to the L-BFGS-B algorithm are derived in Appendix B.

In algorithm 1, steps $1$ and $2$ are used to produce good initial estimates of the source and detector PSF parameters for initializing the optimization in step $3$. In step $1$, only source PSF parameters are estimated from half of all radiographs with the highest $ODD_{k}/SOD_{k}$ where source blur is significant if not dominant. In step $2$, only detector PSF parameters are estimated from half of all radiographs with the lowest $ODD_{k}/SOD_{k}$ where detector blur is significant if not dominant. In step $3$, all parameters from the source PSF, detector PSF, and transmission function are jointly optimized to arrive at the final estimates.

V Deblurring Algorithms

Deblurring is the process of reducing blur in radiographs using computational algorithms. We focus on using Wiener filter [25, 28] and regularized least squares deconvolution [27, 26, 28] for deblurring radiographs. These algorithms take as input the convolution of all PSFs given by,

[TABLE]

which is a function of the source to object distance $SOD_{k}$ and object to detector distance $ODD_{k}$ of the input radiograph $I_{k}(i,j)$. In equation (20), $\hat{p}^{(s)}_{k}(i,j)$ and $\hat{p}^{(d)}(i,j)$ are obtained by substituting the estimated values $\hat{s}_{sx},\hat{s}_{sy},\hat{s}_{d1},\hat{s}_{d2},$ and $\hat{q}$ from step 3 of algorithm 1 in equations (9) and (11). The PSF $\hat{p}^{(m)}_{k}(i,j)$ is assumed to be a discrete delta function (equation (14)) since FWHM of motion blur was determined to be much smaller than the pixel width using the analysis in [32, 33, 34].

V-A Wiener Filter

Wiener filter [25, 28] reduces blur by deconvolving the convolution of all PSFs, $\hat{p}_{k}\left(i,j\right)$ in equation (20), from the blurred radiograph $I_{k}(i,j)$. Deconvolution is implemented in Fourier space by dividing the Fourier transform of the radiograph with the Fourier transform of the PSFs. To reduce noise, regularization is used to enforce smoothness. We use the function skimage.restoration.wiener in the python package scikit-image [40] that implements the method in [25].

V-B Regularized Least Squares Deconvolution (RLSD)

In RLSD [27, 41, 28, 26], we solve the following optimization problem to deblur a radiograph $I_{k}(i,j)$,

[TABLE]

where $\hat{T}_{k}(i,j)$ is the deblurred radiograph, $\beta$ is the regularization parameter, and $\mathcal{N}$ is the set of all pairs of neighboring pixel indices i.e., $((i,j),(m,n))\in\mathcal{N}$ if pixel $(m,n)$ lies within a $3\times 3$ neighborhood of pixel $(i,j)$. The regularization weight parameter $\tilde{w}(i,j,m,n)\propto\frac{1}{\sqrt{(i-m)^{2}+(j-n)^{2}}}$ is inversely proportional to the distance between neighboring pixels and normalized such that $\sum_{(m,n)\in\mathcal{N}_{(i,j)}}\tilde{w}(i,j,m,n)=1$, where $\mathcal{N}_{(i,j)}$ is the set of all voxel indices that are neighbors to voxel index $(i,j)$. For simplicity, the weight parameter $w_{k}(i,j)$ is chosen to be $1$ for all $(i,j)$ in our experiments. Ideally, $w_{k}(i,j)$ should be set such that it is inversely proportional to the variance of noise in $I_{k}(i,j)$.

The regularization function in equation (21) enforces smoothness in $T_{k}(i,j)$ by penalizing the difference in values between neighboring pixels [41, 42]. The optimization problem in equation (21) is solved using the L-BFGS-B algorithm [37].

VI Results

In this section, we will estimate the blur model for a Zeiss Xradia 510 Versa X-ray imaging system. We will compare the efficacy of our blur model against conventional approaches using quantitative metrics. After estimating the blur model, we use Wiener filtering and RLSD (described in section V) to deblur radiographs.

The Versa is a commercial micro-CT system that consists of a transmissive X-ray tube with a Tungsten target anode that exhibits bremsstrahlung X-ray spectral characteristics. The accelerating voltage was selected to be $160kV$ and the tube current was $62.5\mu A$. This resulted in an average flux of $1.37\times 10^{11}counts\,mA^{-1}nstr^{-1}sec^{-1}$. The detector consists of an optically coupled Thallium-doped Cesium Iodide scintillator and a $5MP$ CCD with a $13.5\mu m$ pixel width, $16$-bit depth, and a maximum dark current noise of about $12\,counts\,sec^{-1}$. Bright field images (radiographs without object) $I_{b}(i,j)$ and dark field images (radiographs with X-rays off) $I_{d}(i,j)$ were acquired to appropriately normalize each radiograph image of the object using equation (3). All radiographs were acquired using a $20\times$ magnification lens resulting in an effective pixel width of $\Delta=0.675\mu m$.

VI-A Blur PSF Parameter Estimation

Radiographs of a Tungsten plate rollbar are used to estimate the blur model parameters $s_{sx},s_{sy},s_{d1},s_{d2}$, and $q$ by solving the optimization problem in equation (15). The source to detector distance (SDD) for all radiographs was fixed at $71mm$. First, the Tungsten plate edge is oriented in the vertical direction (Fig. 5 (a)) and radiographs are acquired at source to object distances (SOD) of $13mm$, $24.8mm$, $37.5mm$, $50.3mm$, and $60mm$. Next, the edge is oriented in the horizontal direction and radiographs are acquired at SOD of $12mm$, $24.8mm$, $37.5mm$, $50.3mm$, and $65.3mm$. Since SDD is $71mm$, the object to detector distance (ODD) for each radiograph is given by ($71-$SOD)$mm$. For each radiograph, the Tungsten plate is oriented such that its edge is slightly tilted away from the horizontal or vertical directions as shown in Fig. 5 (a). This is done to ensure that the Tungsten edge is never parallel to a row or column of detector pixels.

The traditional approach to estimating the PSF of X-ray source blur involves assuming that there is no detector blur. To demonstrate this approach, we estimate only the source and transmission parameters, $s_{sx},s_{sy},l_{k},$ and $h_{k}$ in equation (15), using a single set of horizontal and vertical edge radiographs while assuming there is no detector blur, which is enforced by setting

[TABLE]

The estimated values of X-ray source parameters for various source to object distances (SOD) are shown in Table I. The SOD of radiographs used during estimation are shown in the first two columns of Table I. Since FWHMs $(W_{sx},W_{sy})$ are more easily interpretable than scale parameters $(s_{sx},s_{sy})$, we show the FWHMs instead of the scale parameters in the third and fourth columns of Table I. We see that the estimated values for the FWHMs consistently increase with increasing values of SOD in Table I since the source FWHMs were estimated while assuming there is no detector blur. Hence, the detector blur in the input radiographs get interpreted as X-ray source blur, which causes the FWHM of X-ray source blur to increase when SOD is increased.

The traditional approach to detector blur estimation involves assuming that there is no source blur. To evaluate this approach, we estimate only the detector and transmission parameters, $s_{d1},s_{d2},q,l_{k},$ and $h_{k}$ in equation (15), using a single set of horizontal and vertical edge radiographs while assuming there is no source blur, which is enforced by setting

[TABLE]

The estimated parameters of detector blur are shown in Table II. The SOD of radiographs used during estimation are shown in the first two columns of Table II. The estimated detector blur parameters are shown in the last three columns of Table II. In this case, we see that the estimated value for FWHM $W_{d1}$ decreases with increasing SOD since we assumed that there is no X-ray source blur. Note that the first exponential with FWHM of $W_{d1}$ in equation (11) is dominant since its weight given by $q$ is approximately $0.9$ in Table II. As SOD increases and ODD decreases, the FWHM of the source blur on the detector plane reduces. However, since source blur is interpreted as detector blur, the estimated FWHM $W_{d1}$ also decreases as SOD is increased. In contrast, $W_{d2}$ and $q$ do not seem to have a significant dependence on SOD.

To evaluate our proposed approach, we use horizontal edge radiographs at two SODs and vertical edge radiographs at two SODs to simultaneously estimate both source and detector blur parameters using algorithm 1. The estimated values of blur parameters after step 3 of algorithm 1 are shown in Table III. The SOD of radiographs used during estimation are shown in the first two columns of Table III. The estimated PSF parameters are shown in the last five columns of Table III. In this case, the estimation of X-ray source and detector parameters are stable without any noticeable dependence on SOD. Table IV shows the mean (row label Mean), standard deviation (row label Std Dev), and standard deviation as a percentage of the mean (row label $100\frac{\textit{Std Dev}}{Mean}$) of the estimated parameters in Table III computed across various SODs. From Tables III and IV, we can see that by simultaneously accounting for both the X-ray source and detector blur, we are able to perform stable estimation of all PSF parameters.

The source PSF as evaluated in the plane of the source and detector are shown in Fig. 7. The PSF is obtained by substituting the mean values of $W_{sx}$ and $W_{sy}$ shown in Table IV in equations (7) and (9). Since $W_{sx}$ and $W_{sy}$ are approximately the same, we can conclude that the source PSF is approximately circular in shape. Using the method in [21], the manufacturer of the X-ray system estimated the source PSF to have a FWHM of $3.63\mu m$ using calibration data at a SOD of $10mm$ and ODD of $30mm$. For this case, the ratio ODD/SOD lies in between the data acquisition parameters of the second and third rows in Table I. Since the manufacturer’s estimate did not account for detector blur, their values are a good match with our estimates that also ignore detector blur.

The detector PSF is shown in Fig. 8 and is evaluated by substituting the mean values of $W_{d1}$, $W_{d2}$, and $q$ shown in Table IV in equation (11). Since the FWHM $W_{d1}$ of the first exponential density in equation (11) with the high weight of $q$ is only a few pixels wide (pixel width being $0.675\mu m$), we can deduce that most of the detector blur is limited to only a few pixels. Alternatively, the large value of FWHM $W_{d2}$ of the second exponential density with the smaller weight of $(1-q)$ suggests that there is significant leakage of electric charge all the way to the corners of the detector array.

Fig. 9 and Tables I, II, III, and V validate the necessity to simultaneously optimize the source blur, detector blur, and transmission function. While step 1 and step 2 of algorithm 1 independently estimate source and detector blur, step 3 estimates both forms of blur simultaneously. The estimated parameters before and after step (3) of algorithm 1 are shown in Table V. We can see that $W_{d2}$ is significantly larger after step (3) than before step (3). This increase in $W_{d2}$ is due to the simultaneous estimation of $(l_{k},h_{k})$ of the transmission function along with $s_{d2}$ of detector blur in step 3 of algorithm 1. To obtain a better fit, the optimizer decreases $l_{k}$, increases $h_{k}$, which in turn allows for $W_{d2}$ to increase. Since $W_{d2}$ models the long tails of the detector PSF in equation (11), it is most sensitive to changes in $l_{k}$ and $h_{k}$. Fig. 9 (b-d,f-h) shows line profile plots that compare the input radiograph $I_{k}(i,j)$ with the prediction estimate $\bar{I}_{k}(i,j)=T_{k}(i,j)\ast p^{(s)}_{k}(i,j)\ast p^{(d)}(i,j)\ast p^{(m)}_{k}(i,j)$ before and after step 3 of algorithm 1. The improvement in fit after step (3) as demonstrated in Fig. 9 (c,g) is due to the change in parameters $W_{sx}$, $W_{sy}$, and $W_{d1}$. Similarly, the improvement in fit after step (3) in Fig. 9 (b,d,f,h) is due to the change in parameters $W_{d2}$, $l_{k}$, and $h_{k}$.

Fig. 9 (j-l) validates the necessity to model detector PSF as a mixture of two exponential density functions with FWHMs $W_{d1}$ and $W_{d2}$ (equation (11)). The first exponential function with FWHM $W_{d1}$ and a high weight of $q$ ($q\approx 0.9$ in our experiments) models the blur close to the sharp edge that only span a few pixels as shown in Fig. 9 (k). The second exponential function with FWHM $W_{d2}$ and a low weight of $(1-q)$ ($1-q\approx 0.1$) models the slow variation in intensity values further away from the sharp edge as shown in Fig. 9 (j,l). From Fig. 9 (j,l), we can see that the detector PSF is not able to fit the slowly varying intensity that span several hundreds of pixels when $q$ is not estimated and fixed at $q=1$. Fixing $q=1$ is equivalent to modeling detector PSF as consisting of only one exponential density function with a FWHM of $W_{d1}$.

VI-B Validating the Blur Model

The optimal choice of density function used to model the source blur in equation (7) and the detector blur in equation (11) is system dependent. In Table VI, we compare the performance of various choices of blur model using quantitative evaluation metrics. Horizontal and vertical edge radiographs at SODs of $37.5mm$ and $50.3mm$ were used when simultaneously estimating both source and detector blurs. When estimating only source blur, only radiographs at SOD of $37.5mm$ were used since source blur is dominant at the lower SOD. When estimating only detector blur, only radiographs at SOD of $50.3mm$ were used since detector blur is dominant at the higher SOD. The $1^{st}$ column of Table VI indicates whether the PSF model was chosen to be either Gaussian ($r=2$) or exponential ($r=1$). The $2^{nd}$ column indicates if source blur was modeled (marked as “Yes”) or ignored (marked as “No”). The $3^{rd}$ column indicates if detector blur was modeled (marked as “Yes”) or ignored (marked as “No”). To determine the most appropriate blur model, we evaluate the blur prediction and deblurring performance on radiographs that were acquired at SOD values different than those used for blur estimation. Computation of the performance evaluation metric, root mean squared error (RMSE), in the $4^{th},6^{th}$, and $8^{th}$ columns is over both the horizontal and vertical edge radiographs at SOD of $24.8mm$. Similarly, the RMSEs in the $5^{th},7^{th}$, and $9^{th}$ columns are computed over the horizontal edge radiograph at SOD of $65.3mm$ and vertical edge radiograph at SOD of $60mm$.

The accuracy of blur prediction is quantified in the $4^{th}$ and $5^{th}$ columns under the header “Prediction $\bar{I}_{k}(i,j)$” by computing the RMSE between the measured radiographs $I_{k}(i,j)$ and its predictions $\bar{I}_{k}(i,j)=T_{k}(i,j)\ast p^{(s)}_{k}(i,j)\ast p^{(d)}(i,j)\ast p^{(m)}_{k}(i,j)$ using the blur model. In order to compute $\bar{I}_{k}(i,j)$, we must compute the transmission function $T_{k}(i,j)$ by setting appropriate values for $l_{k},h_{k}$ in equation (5) for radiographs that were not used during blur estimation. As an approximate solution to this problem, the $l_{k}$ and $h_{k}$ for horizontal/vertical edge prediction $\bar{I}_{k}(i,j)$ are chosen to be the average of the corresponding values for the horizontal/vertical edge transmission functions estimated during blur optimization. The lowest RMSE is obtained along the last row of Table VI, which corresponds to simultaneous estimation of both source and detector blurs while assuming an exponential density function for both PSFs.

Next, we compare the deblurring performance using Wiener filtering for the different choices of blur model. Since Wiener filter uses a regularization parameter that trades off noise with resolution, it is not meaningful to compare the quality of the edge when the deblurred radiographs under comparison have varying noise levels. Hence, the regularization parameter for Wiener filter is adjusted until the noise level in the radiograph is close to a pre-specified constant $\gamma$. Fig. 10 and $6^{th}$ to $9^{th}$ columns under the headers “Deblur, SD = $\gamma$” in Table VI compare the deblurring performance for various model choices. In Fig. 10, deblurred results are shown for horizontal edge radiograph at a SOD of $24.8mm$ while blur estimation was performed using radiographs at SODs of $37.5mm$ and $50.3mm$. The noise level $\gamma$ is computed as the average of the standard deviations within rectangular boxes along each of the four corners of a radiograph (similar to the blue squares in Fig. 10 (a)). The line profiles in Fig. 10 (b,d) are along the red line in Fig. 10 (a,c) respectively. The RMSE for deblurring comparison in Table VI is computed between the deblurred radiographs and its corresponding transmission functions given by equation $\eqref{eq:transvsideal}$. The lowest RMSE under the header “Deblur, SD = $\gamma$” is achieved along the last row of Table VI. In Fig. 10, we can see that using exponential over Gaussian density function results in less intense ringing artifacts in (b,d). Thus, the best choice of blur model is to use exponential density function ($r=1$) over Gaussian density function ($r=2$) while simultaneously modeling the effects of both X-ray source and detector blur. Fig. 11 shows that reducing noise and ringing artifacts in deblurred radiograph using regularization also reduces the sharpness of the Tungsten edge.

Table VIII compares the variation in deblurring performance as a consequence of the SOD dependent variation in estimated PSF parameters in Table III. Both source and detector blur PSFs were modeled using exponential density function ($r=1$) and estimated using algorithm 1. From Table VIII, we can see that irrespective of the SOD used to estimate blur, the RMSE is close to the corresponding value in the last row of Table VI (under column header “Deblur, SD = $\gamma$” and sub-header “24.8mm/24.8mm”). Remarkably, this RMSE for a joint source and detector blur model with exponential density is always lower than the RMSE for all other model choices considered in Table VI. Table VIII computes the mean, standard deviation, and standard deviation as a percentage of the mean for each column in Table VIII. The standard deviation as a percentage of the mean shown in the last row of Table VIII is in the same general range as the corresponding value computed for the PSF parameters in the last row of Table IV.

VI-C Deblurring Algorithms

A variety of deblurring algorithms that use the estimated X-ray source and detector PSFs can be used to deblur radiographs. We acquired radiograph (Fig. 12 (c,d)) of a star shaped object at a SOD of $10mm$ and ODD of $61mm$. This star shaped test object was composed of a $1\mu m$ thick Tungsten layer on a SiN membrane. The blurry radiograph in Fig. 12 (c) is then deblurred using Wiener filtering and RLSD algorithm. The regularization of both Wiener filter and RLSD methods are adjusted until the noise variance in the red square region in Fig. 12 (e) and (g) are the same and equal to $0.005$. The deblurred radiographs shown in Fig. 12 (e-h) are much sharper than the input radiograph in Fig. 12 (c,d). A scanning electron microscopy (SEM) image of the star artifact in Fig. 12 (a,b) show the slits inside each spoke of the star pattern. By comparing Fig. 12 (f) with Fig. 12 (h), we see that the sharpness of RLSD image is better than the Wiener image since the slits in the spokes of the star pattern are more clear in the RLSD image.

VII Conclusion

In this paper, we presented a method to estimate both X-ray source and detector blur from radiographs of a Tungsten plate rollbar. Importantly, our method is able to distentangle and estimate the parameters of both X-ray source and detector blur from radiographs that are simultaneously blurred by both forms of blur. We show that blur estimation can be performed using horizontal edge and vertical edge radiographs, each of which is measured at two different values of the ratio of object to detector distance (ODD) and source to object distance (SOD). Using the estimated blur model, we demonstrated the ability to deblur radiographs using various deblurring algorithms.

Appendix A Computing the Ideal Transmission Function

For every radiograph $I_{k}(i,j)$, the ideal transmission function $\tilde{T}_{k}(i,j)$ is estimated from $I_{k}(i,j)$ using traditional image processing algorithms. The first step to computing $\tilde{T}_{k}(i,j)$ is to determine the location of the Tungsten plate’s boundary. Even with the rolled edge, the plate is designed such that the transition region from Tungsten to air is no more than one pixel thick. But, from Fig. 5 (a), we can see that it is very difficult to determine the exact pixel locations of the plate’s boundary due to blur from the X-ray source and detector. Hence, we scale the radiograph $I_{k}(i,j)$ such that it is in the range of [math] to $1$ and assume that the edge lies along a iso-valued contour with a level value of $0.5$. Such a iso-valued contour is estimated using the marching squares algorithm [43] implemented in the python package scikit-image [40].

The ideal transmission function $\tilde{T}_{k}(i,j)$ is assigned a value of [math] for the pixels belonging to the Tungsten plate and $1$ for the pixels with no plate (or air pixels) since the plate is designed to completely attenuate all X-rays. For the edge boundary pixels, their values are linearly interpolated given a value of $0.5$ along the iso-valued contour and neighboring pixel values of 1 and 0. The iso-valued contour produced by the marching squares algorithm is in the form of a list of real valued $(k,l)$ coordinates. However, the pixel coordinates $(i,j)$ of $\tilde{T}_{k}(i,j)$ are integer valued. Hence, the values at edge pixel coordinates $(\lfloor k\rfloor,\lfloor l\rfloor)$, $(\lfloor k\rfloor,\lceil l\rceil)$, $(\lceil k\rceil,\lfloor l\rfloor)$, and $(\lceil k\rceil,\lceil l\rceil)$ are determined by interpolation.

To prevent aliasing during convolution, $\tilde{T}_{k}(i,j)$ is padded to three times the size of $I_{k}(i,j)$ by a special padding procedure that takes into account the orientation of the Tungsten plate (Fig. 5 (b)). The same amount of padding is applied to the top, bottom, left, and right edge of the images. First, a straight line is fit to the sharp edge of the Tungsten plate (Fig. 5 (a)). We extend the straight lines outside the image to account for the plate extending outside the image’s field of view. The straight line will split the padded image into two regions, one with the Tungsten plate and another without the plate. Every padded pixel will have a value of [math] or $1$ depending on whether it lies in the region with the plate or without the plate.

Appendix B Gradient Computation

To solve the optimization problem in equations (17), (18), and (19), we use the L-BFGS-B algorithm [37]. Similar to many optimization algorithms such as gradient descent and conjugate gradient, L-BFGS-B needs a routine to calculate the gradient of the objective function in equations (17), (18), and (19) with respect to the variables that are being optimized. For example, to solve the optimization problem in (17), L-BFGS-B needs to know the gradient of the objective function $\sum_{k\in\Omega_{s}}E_{k}$ with respect to $s_{sx}$ and $s_{sy}$.

We will derive the gradient with respect to all variables for the objective function in (19) (same as (15)). Gradients for solving (17) and (18) will be a straightforward extension of the gradient derived for (19). Let the objective function be denoted by $f$ i.e, $f=\sum_{k=1}^{K}E_{k}$. To solve (19), we need the gradient of $f$, which is a vector consisting of the partial derivatives $\frac{\partial f}{\partial s_{sx}}$, $\frac{\partial f}{\partial s_{sy}}$, $\frac{\partial f}{\partial s_{d1}}$, $\frac{\partial f}{\partial s_{d2}}$, $\frac{\partial f}{\partial q}$, $\frac{\partial f}{\partial l_{k}}$, and $\frac{\partial f}{\partial h_{k}}\,\,\forall k$.

To compute the partial derivatives, we use the chain rule and quotient rule of calculus. First, we will partially compute only that part of the derivative which is common to all parameters irrespective of whether it belongs to the source PSF, the detector PSF, or the transmission function. Let $v$ represent any one parameter among $s_{sx}$, $s_{sy}$, $s_{d1}$, $s_{d2}$, $q$, $l_{k}$, and $h_{k}$. Then, the derivative of $f$ with respect to $v$ is,

[TABLE]

where $\bar{I}_{k}(i,j)=T_{k}(i,j)\ast p^{(s)}_{k}(i,j)\ast p^{(d)}(i,j)\ast p^{(m)}_{k}(i,j)$. In equation (22), every term and operator except $\frac{\partial\bar{I}_{k}(i,j)}{\partial v}$ is independent of whether $v$ is $s_{sx}$, $s_{sy}$, $s_{d1}$, $s_{d2}$, $q$, $l_{k}$, or $h_{k}$. Next, we shall expand the derivative $\frac{\partial\bar{I}_{k}(i,j)}{\partial v}$. This requires us to consider source, detector, and transmission function parameters separately.

B-A Variable $v$ is a source PSF parameter

Let $v$ be one among the source PSF parameters $s_{sx}$ or $s_{sy}$. For compactness of representation, we shall assume $g_{s,k}(i,j)=\exp\left(-\Delta^{r}\frac{SOD^{r}_{k}}{ODD^{r}_{k}}\left(i^{2}s^{2}_{sx}+j^{2}s^{2}_{sy}\right)^{\frac{r}{2}}\right)$ and $Z_{s,k}=\sum_{i,j}g_{s,k}(i,j)$ from equation (10). Then, using the quotient rule,

[TABLE]

where $\ast$ denotes discrete 2D convolution, $\partial Z_{s,k}/\partial v=\sum_{i,j}\partial g_{s,k}(i,j)/\partial v$, and

[TABLE]

B-B Variable $v$ is a detector PSF parameter

Let $v$ be one among the detector PSF parameters $s_{d1}$, $s_{d2}$, or $q$. For compactness of representation, we shall assume $g_{d1}(i,j)=\exp\left(-s^{r}_{d1}\Delta^{r}\left(i^{2}+j^{2}\right)^{\frac{r}{2}}\right)$, $g_{d2}(i,j)=\exp\left(-s^{r}_{d2}\Delta^{r}\left(i^{2}+j^{2}\right)^{\frac{r}{2}}\right)$, $Z_{d1}=\sum_{i,j}g_{d1}(i,j)$, and $Z_{d2}=\sum_{i,j}g_{d2}(i,j)$ from equations (12) and (13). Using the quotient rule, we get,

[TABLE]

where $\ast$ denotes discrete 2D convolution,

[TABLE]

$\partial g_{d1}(i,j)/\partial s_{d1}=-g_{d1}(i,j)\Delta^{r}\left(i^{2}+j^{2}\right)^{\frac{r}{2}}rs_{d1}^{r-1}$, $\partial g_{d2}(i,j)/\partial s_{d2}=-g_{d2}(i,j)\Delta^{r}\left(i^{2}+j^{2}\right)^{\frac{r}{2}}rs_{d2}^{r-1}$, $\partial Z_{d1}/\partial v=\sum_{i,j}\partial g_{d1}(i,j)/\partial v$, and $\partial Z_{d2}/\partial v=\sum_{i,j}\partial g_{d2}(i,j)/\partial v$.

B-C Variable $v$ is a transmission function parameter

Let $v$ be one among the transmission function parameters $l_{k}$ or $h_{k}$. Then,

[TABLE]

[TABLE]

Acknowledgment

LLNL-JRNL-760284. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LDRD funding with tracking numbers 16-ERD-006 and 19-ERD-022 were used for this project. We thank Kyle Champley from LLNL for useful discussions. We also thank Markus Baier and Professor Simone Carmignato from University of Padova for fabricating and providing the star pattern test object utilized in this manuscript.

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.