Analysis of reconstruction from noisy discrete generalized Radon data

Tomographic imaging, which includes X-ray, ultrasound, seismic, and many other tomographies, is commonly used in a wide range of applications such as medical, industrial, security, and others. The task of tomographic image reconstruction can be formulated as follows. Let $\mathcal{R}:\mathcal{X}\to\mathcal{V}$ denote a linear integral transform, where $\mathcal{X}\subset\mathbb{R}^{n}$ is the image domain and $\mathcal{V}\subset\mathbb{R}^{n}$ is the data domain. Let $f$ be a function supported in $\mathcal{X}$, which represents the object being scanned. Tomographic data represent integrals of $f$ along a family of manifolds (e.g., lines, spheres, etc.). The manifolds are parametrized by points $\omega\in\mathcal{V}$, so the data are $(\mathcal{R}f)(\omega)$, $\omega\in\mathcal{V}$. In a continuous setting, the goal is to reconstruct $f(x)$, $x\in\mathcal{X}$ (or some information about $f$, e.g., its singularities) from $(\mathcal{R}f)(\omega)$, $\omega\in\mathcal{V}$.

In practice, data are always discrete and contain noise. So the task becomes to reconstruct an approximation to $f$ from discrete data $g_{j}:=(\mathcal{R}f)(\epsilon j)+\eta_{j}$, $\epsilon j\in\mathcal{V}$, $j\in\mathbb{Z}^{n}$. Here $\epsilon>0$ represents the data step size and $\eta_{j}$’s represent additive noise. For simplicity, we assume that the step size is the same along each direction. Usually, two most important factors that affect quality of reconstructed images are data discretization and strength of noise in the data.

Effects of data discretization on various aspects of image quality, e. g. spatial resolution, aliasing artifacts, are now well-understood [33, 34, 14, 40]. In particular, in a series of articles, the author developed an approach, called local reconstruction analysis (LRA), to analyze the resolution with which the singularities (i.e., jump discontinuities) of $f$ are reconstructed from discrete tomographic data in the absence of noise, see [22, 23, 20] and references therein. In those articles the singularities of $f$ are assumed to lie on a smooth curve in $\mathbb{R}^{2}$ (and surface in $\mathbb{R}^{n}$) or a rough curve in $\mathbb{R}^{2}$, denoted $\mathcal{S}$. In [24], LRA theory was advanced to include analysis of aliasing (ripple artifacts) in $\mathbb{R}^{2}$ at points away from $\mathcal{S}$.

The main idea of LRA is to obtain a simple formula to accurately approximate an image reconstructed from discrete data, $f_{\epsilon}^{\text{rec}}$, in an $\epsilon$-neighborhood of a point, $x_{0}$. For example, let $f$ be a real-valued function in $\mathbb{R}^{2}$, $\mathcal{S}\subset\mathbb{R}^{2}$ be a smooth curve, and $f$ have a jump across $\mathcal{S}$. Let $f_{\epsilon}^{\text{rec}}$ be a filtered backprojection (FBP) reconstruction of $f$ from discretely sampled Radon transform (RT) data with sampling step size, $\epsilon$. Under some mild conditions on $\mathcal{S}$, it is shown in [21, 23] that one has

where the limit is uniform with respect to $\check{x}\in D$ for any bounded set $D$. Here $\Delta f(x_{0})$ is the value of the jump of $f$ across $\mathcal{S}$ at $x_{0}$ and DTB, which stands for the Discrete Transition Behavior, is an easily computable function independent of $f$. The DTB function depends only on the curvature of $\mathcal{S}$ at $x_{0}$. When $\epsilon$ is sufficiently small, the right-hand side of (1.1) is an accurate approximation of $f_{\epsilon}^{\text{rec}}$ and the DTB function describes accurately the smoothing of the singularities of $f$ in $f_{\epsilon}^{\text{rec}}$.

As is seen from (1.1), LRA provides a uniform approximation to $f_{\epsilon}^{\text{rec}}(x)$ in domains of size $\sim\epsilon$, which is of the same order of magnitude as the sampling step size $\epsilon$. Given the DTB function, one can study local properties of reconstruction from discrete data, such as spatial resolution and strength of artifacts.

The first extension of LRA to reconstruction from noisy data has been done in [2]. The authors consider the model $g=\mathcal{R}f+\eta$, where $\mathcal{R}$ is the classical 2D RT, $\mathcal{R}f$ is the measured sinogram with entries $(\mathcal{R}f)_{j,k}$, and $\eta$ is the noise. The entries, $\eta_{j,k}$, of $\eta$ are assumed to be independent, but not necessarily identically distributed. Due to linearity of $\mathcal{R}^{-1}$, the reconstruction error is $\mathcal{R}^{-1}\eta$. Let $N_{\epsilon}^{\text{rec}}(x)$ denote the reconstruction error, i.e. the FBP reconstruction only from $\eta_{j,k}$. Similarly to previous works on LRA (cf. also (1.1)), the authors consider $O(\epsilon)$-sized neighborhoods around any $x_{0}\in\mathbb{R}^{2}$. Let $C(D)$ be the space of continuous functions on $D$, where $D\subset\mathbb{R}^{2}$ is a bounded domain. The main result of [2] is that, under suitable assumptions on the first three moments of the $\eta_{j,k}$, the boundary of $D$, and $x_{0}$, the following limit exists:

Here, $N^{\text{rec}}$ and $N_{\epsilon}^{\text{rec}}$ are viewed as $C(D)$-valued random variables (random fields), and the limit is understood in the sense of probability distributions. It is proven also that $N^{\text{rec}}(\check{x};x_{0})$ is a zero mean Gaussian random field (GRF) and its covariance is computed explicitly. The paper also contains numerical experiments, which show an excellent match between theoretical predictions and simulated reconstructions.

In this paper we extend the results of [2] to a wide class of generalized RTs $\mathcal{R}$, which act in $\mathbb{R}^{n}$ for any $n\geq 2$ and integrate over submanifolds of any codimension $N$, $1\leq N\leq n-1$. Also, we allow for a fairly general reconstruction operator $\mathcal{A}$. The main requirement is that $\mathcal{A}$ be a Fourier integral operator (FIO) with a phase function, which is linear in the phase variable. Similarly to [2] we consider the model $g=\mathcal{R}f+\eta$ and show that the reconstruction error $N_{\epsilon}^{\text{rec}}=\mathcal{A}\eta_{j,k}$ satisfies (1.2), where $D\subset\mathbb{R}^{n}$ is a bounded domain and $\eta_{j,k}$ are independent, but not necessarily identically distributed, random variables. As in [2], $N^{\text{rec}}$ and $N_{\epsilon}^{\text{rec}}$ are viewed as $C(D)$-valued random variables, and the limit is understood in the sense of probability distributions. We prove that $N^{\text{rec}}$ is a zero mean GRF and explicitly compute its covariance. We also present a numerical experiment with a cone beam transform in $\mathbb{R}^{3}$, which shows an excellent match between theoretical predictions and simulated reconstructions.

Equations (1.1) and (1.2) together provide an explicit and accurate local approximation to the image reconstructed from noisy discrete data, $f_{\epsilon}^{\text{rec}}=\mathcal{A}(\mathcal{R}f+\eta)$. Eq. (1.1) describes the effect of edge smoothing due to data discretization, and (1.2) describes the reconstruction error due to random noise. The combined approximation is uniform over domains of size $\sim\epsilon$, i.e. precisely at the scale of native resolution enabled by available data. Results of this kind are more useful for localized analysis of tomographic reconstruction than more common global analyses, which estimate some global error norm (e.g., $L^{2}$). Formulas provided by LRA, such as (1.1) and (1.2), can be used to study resolution of reconstruction, perform statistical inference about detectability of small details in a reconstructed image, and for many other tasks.

As a first example we consider the task of detection and assessment of lung tumors in CT images (see also [2]). Typically, malignant lung tumors have rougher boundaries than benign nodules [10, 11]. Therefore the roughness of the nodule boundary is a critical factor for accurate diagnosis. Reconstructions from discrete X-ray CT data produce images in which the singularities are smoothed to various degrees. A rough boundary of the tumor may appear smoother in the reconstructed image than it really is. This can lead to a cancerous tumor being misdiagnosed as a benign nodule. Likewise, due to the random noise in the data, the smooth boundary of a benign nodule may appear rougher than it actually is, which may again result in misdiagnosis. This example illustrates the need to accurately quantify the effects of both data discretization and random noise on local (i.e., near the tumor boundary) tomographic reconstruction.

Another example is the use of CT by the energy industry for imaging of rock samples extracted from wells. The reconstructed image is segmented to identify as accurately as possible the pore space (i.e., voids) inside the sample. This is very challenging, because pore boundaries are fractal, they possess features that are below the scanner native resolution. Then one uses numerical fluid flow simulations inside the identified pore space to compute permeability of the sample [5]. The obtained results are used for formation evaluation and improving of oil recovery [36]. Thus, the key step that affects the accuracy of the entire workflow from scanning to fluid flow simulation is CT image segmentation. Errors in locating pore boundaries and, consequently, inaccurate pore space identification lead to incorrect flow simulation and errors in computed permeability [9, 37, 38]. Therefore, as in the previous example, precise quantification of the local (near the boundaries) effects of data discretization and noise on the reconstructed image is of paramount importance.

We now discuss related existing literature concerning reconstruction in a stochastic setting and compare it with our findings. A kernel-type estimator of $f$ has been derived in [26, 27] in the setting of a tomography problem with additive random errors in observations. Its minimax optimal rate of convergence to $f$, i.e. the ground truth function, has been established both at a fixed point and in a global $L^{2}$ norm. The rate of convergence for the maximal deviation of an estimator from its mean has been obtained for similar kernel-type estimators in [6]. The accuracy of pointwise asymptotically efficient kernel estimator in terms of minimax risk of a probability density $f$ from noise-free RT data sampled on a random grid has been derived in [7, 8].

An essential common feature of all the works cited above is that convergence is established by incorporating into the reconstruction/estimation algorithm of additional smoothing at a scale $\delta$ significantly larger than the data step size. This smoothing leads to a notable loss of resolution in practical applications. This is also the reason why the cited works assume that the function $f$ being estimated is sufficiently smooth. To compare with our results, we mention two points. First, none of these papers obtain the probability distribution of the reconstructed noise (i.e., the error in the reconstruction), either pointwise or in a domain. Second, our reconstruction is performed and investigated at native resolution.

Reconstruction in the framework of Bayesian inversion has been explored as well [32, 39, 43]. Global inversion in the case when continuous data are corrupted by a Gaussian white noise is investigated in [32]. Using a Gaussian prior (i.e., with Tikhonov regularization), the authors establish the asymptotic normality of the posterior distribution and of the MAP estimator for observables of the kind $\int f(x)\psi(x)\text{d}x$, where $\psi\in C^{\infty}$ is a test function. This means that reconstruction is investigated at the scales $\sim 1$. In contrast, our results quantify the pointwise reconstruction error uniformly over regions of size $O(\epsilon)$. Various other aspects of Bayesian inversion are discussed in [39, 43].

Analysis of reconstruction errors using semiclassical analysis is developed in [41]. The goal of the paper is to analyze empirical spatial mean and variance of the noise in the inversion for a single experiment in the limit as the sampling rate $\epsilon$ goes to zero. In contrast, the emphasis of our paper is on the study of the reconstruction error across multiple reconstructions. We obtain the entire probability density function (PDF) of the error in the limit as $\epsilon\to 0$ (as opposed to only the mean and variance).

Analysis of noise in reconstructed images is an active area in more applied research as well (see [44, 12] and references therein). In this research, the proposed methodologies mostly combine numerical and semi-empirical approaches. Theoretical derivation of the reconstructed noise PDF in small neighborhoods is not provided.

The paper is organized as follows. In section 2, we describe the setting of the problem and main assumptions. In section 3, we state the main results, which are broken down into three theorems. In section 4 we illustrate how to check the key assumptions for a few common CT geometries. The beginning of the proof of the first result, theorem 3.1, is in section 5. In this section we study the contribution of the leading order term of $\mathcal{A}$ to the reconstructed image, $N_{\epsilon}^{\text{rec}}(x_{0}+\epsilon\check{x})$. In section 6 we study the contribution of the lower order terms of $\mathcal{A}$ to the reconstruction. In section 7 we finish the proof of theorem 3.1 and prove theorem 3.2, our second main result. In section 8 we prove our last result, theorem 3.4. A numerical experiment is described in section 9. Finally, the proofs of some auxiliary lemmas are in Appendices A–E.

We begin with the classical RT in $\mathbb{R}^{2}$. The following conventional data parametrization is used:

where $\vec{\alpha}=(\cos\alpha,\sin\alpha)$. Then $\partial_{p}\Phi(x,\alpha,p)=1$, so Assumption 2.4(3) is satisfied. From (2.15), $Y_{1}(x_{0},\xi)=\{\alpha\in[0,2\pi):\vec{\alpha}\cdot x_{0}=0\}$, which is a set consisting of two points, i.e. it has the upper box-counting dimension of zero, as long as $x_{0}\not=0$. Also, $Y_{2}(x_{0},\check{x})=\{\alpha\in[0,2\pi):\vec{\alpha}\cdot\check{x}=0\}$ (cf. (2.17)), which is a set of measure zero. Thus, Assumptions 2.8(1,3) are satisfied.

Next consider the classical RT in $\mathbb{R}^{3}$. In this case we use the conventional parametrization

Then $\partial_{p}\Phi(x,\alpha,\theta,p)=1$, so Assumption 2.4(3) is satisfied. By (2.15), we compute the Hessian:

Here we have assumed without loss of generality that $\xi=1$. Clearly, $\text{det}(\partial_{(\alpha,\theta)}^{2}\Psi)\equiv 0$ if $\hat{x}=0$. Thus, all points $x=(0,0,x_{3})$, $x_{3}\in\mathbb{R}$, are exceptional, because they violate Assumption 2.8(1). Suppose now $\hat{x}\not=0$. After simple transformations we obtain

For each $\theta\in[0,2\pi)$, we can find at most finitely many solutions $|\alpha|<\pi/2$ to $\text{det}(\partial_{(\alpha,\theta)}^{2}\Psi)=0$ and, generically, they depend smoothly on $\theta$. Hence, the sets $Y_{1}(x,\xi)$, $\xi\not=0$, have the upper box-counting dimension of one, as long as $\hat{x}\not=0$. For the cone beam transform, $n=3$ and $N=1$, so Assumption 2.8(1) is satisfied.

Finally, $Y_{2}(x,\check{x})=\{\vec{\Theta}\in\mathcal{S}^{2}:\vec{\Theta}\cdot\check{x}=0\}$ (cf. (2.17)), which is a set of measure zero in the unit sphere. Thus, Assumption 2.8(3) is satisfied as well.

Our next example is the cone beam transform in $\mathbb{R}^{3}$ with a circular source trajectory. The detector coordinates are $(u,v)\in\mathbb{R}^{2}$ (the analog of $z$), the source coordinate is $s\in[0,2\pi)$ (the analog of $y$), see Figure 1, and $n=3$, $N=2$. The source trajectory is given by

The (virtual) detector is flat, passes through the origin, and rotates together with the source. Points on the detector are parametrized as follows:

For a point $x=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$, its stereographic projection from the source $P(s)$ to a point $(u,v)=(U(x,s),V(x,s))$ in the detector plane is given by

The support of $f$ and reconstruction domain $\mathcal{X}$ are contained inside the cylinder $x_{1}^{2}+x_{2}^{2}\leq c<R$. The two key functions are given by

Clearly, $\partial_{(u,v)}\Phi$ is a $2\times 2$ identity matrix, so Assumption 2.4(3) is satisfied.

From (2.15), $Y_{1}(x,\xi)=\{s\in[0,2\pi):\partial_{s}^{2}(\xi\cdot\Psi(x,s)=0\}$, $\xi\in\mathbb{R}^{2}\setminus 0$. Hence, generally, Assumption 2.8(1) is violated if the projection of $x$ onto the detector contains a straight line segment, see Figure 2, left panel. For the circular source trajectory (4.5), this is the case when $x_{3}=0$, i.e. when $x$ is in the plane of the source trajectory, see Figure 2, right panel. As is easy to see, the projection of $x$ to the detector is an ellipse as long as $x_{3}\not=0$. Indeed, from (4.7) we find

Eliminating $\sin s$ and $\cos s$ gives

This is an ellipse, because the reconstruction domain is inside the source trajectory, i.e. $x_{1}^{2}+x_{2}^{2}<R^{2}$. Therefore Assumption 2.8(1) is satisfied if $x_{3}\not=0$.

The meaning of Assumption 2.8(3) is best understood using Figures 1 and 3. We have $Y_{2}(x,\check{x})=\{s\in[0,2\pi):\partial_{t}\Psi(x+t\check{x},s)|_{t=0}=0\}$, where $\Psi(x,s)$ is the projection of $x$ to the detector. It is clear from the figures that the directional derivative is zero only when $\check{x}$ is along the line $L(x,s)$ through $x$ and the source $P(s)$. This implies that Assumption 2.8(3) is violated only if there exists a line $L$ through $x$ such that the intersection of $L$ with the source trajectory contains a line segment as illustrated in Figure 3. Clearly, for a circular trajectory this does not happen.

Change variable $z\to w=\Phi(u,z)$ and then $w\to\check{w}=w/\epsilon$, $\xi\to\hat{\xi}=\epsilon\xi$: