Edge-masked CT image reconstruction from limited data

TL;DR

This paper introduces an iterative CT image reconstruction algorithm that leverages edge masks derived from limited data to improve accuracy and speed, supported by theoretical and empirical evidence.

Contribution

The proposed method combines image domain and statistical techniques, using edge masks to enhance reconstruction from limited data, adaptable to various models and transforms.

Findings

High accuracy with limited data

Fast convergence in iterative reconstruction

Flexible framework for different models

Abstract

This paper presents an iterative inversion algorithm for computed tomography image reconstruction that performs well in terms of accuracy and speed using limited data. The computational method combines an image domain technique and statistical reconstruction by using an initial filtered back projection reconstruction to create a binary edge mask, which is then used in an l2-regularized reconstruction. Both theoretical and empirical results are offered to support the algorithm. While in this paper a simple forward model is used and physical edges are used as the sparse feature, the proposed method is flexible and can accommodate any forward model and sparsifying transform.

Figures

Equations

$$\arg \underset{\mathbf{u}}{\min }\left\{||\mathbf{Ru}-\mathbf{s}|{|}_2^2+\lambda ||\mathbf{Du}|{|}_1\right\}.$$

$$\arg \underset{\mathbf{u}}{\min }||\mathbf{MDu}|{|}_2^2\text{subject to}\mathbf{Ru}=\mathbf{s},$$

$$||\mathbf{MDu}|{|}_2^2=\sum _{i:{\left(\mathbf{Dx}\right)}_i=0}{\left(\mathbf{Du}\right)}_i^2\ge 0,$$

$${\mathbf{M}}_{ii}=\{\begin{array}{cc}1& |{\mathbf{y}}_i|<\tau \\ 0& |{\mathbf{y}}_i|\ge \tau \end{array},$$

$$\arg \underset{\mathbf{u}}{\min }\left\{||\mathbf{Ru}-\mathbf{s}|{|}_2^2+\lambda ||\mathbf{MDu}|{|}_2^2\right\}.$$

Full text

***Victor Churchill, \linkableVictor.A.Churchill.GR@dartmouth.edu **

1 Introduction

One goal in computed tomography (CT) imaging is to reduce dose while maintaining image quality and reconstruction speed. When collecting less data to reduce dose, direct inversion techniques such as filtered back projection (FBP) produce artifacts. Iterative inversion techniques, like those arising from compressed sensing, have shown that a priori information about physical edges in an image can be extremely useful for reconstruction from limited data. For example, assuming edges are sparse in a CT image, one can use total variation (TV) regularization [1], which solves

[TABLE]

Here $\mathbf{R}$ is the CT forward operator111The Radon transform is used for data collection and the inverse Radon transform with a Ram-Lak filter for FBP., $\mathbf{s}$ is the collected sinogram data, $\mathbf{D}$ is the sparsifying transform222The anisotropic 2D TV transform defined by $\mathbf{Du}=\sum_{i,j}|\mathbf{u}_{i+1,j}-\mathbf{u}_{i,j}|+|\mathbf{u}_{i,j+1}-\mathbf{u}_{i,j}|$ is used., and $\lambda$ is the regularization parameter that balances fidelity, edge sparsity, and noise reduction. In many cases, Eq. (1) is able to reconstruct images of similar quality to FBP using less data. The drawback is the reconstruction speed. In general there is no direct solution to Eq. (1), so an iterative convex optimization method must be used. In this paper, a more direct computational method is proposed, similar to the edge-adaptive $\ell_{2}$-regularized least squares scheme [2] that was formulated for recovery from non-uniform Fourier data. The method combines an image domain technique with statistical reconstruction. First, the resulting image from an initial FBP reconstruction is used to create an edge mask. Then, an iterative $\ell_{2}$-regularized reconstruction is performed using the edge mask to regularize away from edges. In what follows, theoretical and empirical evidence is offered to show that with an accurate edge detection, this method has the potential to produce near-perfect reconstructions with extremely limited data at a lower computational cost. In addition, recognizing that the ideal sparse feature for CT images may not be physical edges in the TV domain, this algorithm is flexible and can in fact accommodate a “wavelet edge” or any another non-physical “edge”, where the definition is expanded to mean a nonzero location in a sparse domain.

2 Algorithm

First, the importance of edge locations in image reconstruction is demonstrated by proving that given exact edge locations in the sparsity domain, Eq. (2) returns a unique edge-sparse solution if there is one.

Theorem 1. Assume there is a unique solution $\mathbf{x}$ to $\mathbf{R}\mathbf{u}=\mathbf{s}$ with $||\mathbf{D}\mathbf{u}||_{0}\leq K$ where $||\cdot||_{0}$ counts the number of nonzero elements and $K>0$ is a constant. Let $\mathbf{u}^{*}$ be the solution to the problem

[TABLE]

where $\mathbf{M}$ is a diagonal mask matrix with diagonal element $1$ if $(\mathbf{D}\mathbf{x})_{i}=0$ and [math] otherwise. Then $\mathbf{u}^{*}=\mathbf{x}$.

Proof. We have that

[TABLE]

which is minimized (at zero) for $\mathbf{u}=\mathbf{x}$, since we know $\mathbf{R}\mathbf{x}=\mathbf{s}$. Since $\mathbf{x}$ is the unique solution with $||\mathbf{D}\mathbf{u}||_{0}\leq K$, it must be that $\mathbf{u}^{*}=\mathbf{x}$.$\quad\square$

It is critical to note that this theorem would also hold using the $\ell_{1}$ norm instead of the squared $\ell_{2}$ norm in Eq. (2). That being said, $\ell_{2}$ is appropriate here because $\mathbf{MDu}$ is not sparse but zero, therefore doesn’t require a sparsity-encouraging norm like $\ell_{1}$ to regularize. Therefore, using $\ell_{2}$ gives the same result while being more computationally efficient.

Figure 1 shows an example of this theorem in practice compared with FBP and Eq. (1). In this experiment, the exact edge locations are known so the mask $\mathbf{M}$ is prescribed as in Theorem 1 and reconstructed via Eq. (2). However, only a single-view sinogram (one projection line) is used. Still, the reconstruction is near-perfect. This example is not intended to serve as a rigorous lower bound for the amount of data required to reconstruct given the exact edge locations, but rather it simply demonstrates the potential of Eq. (2) – that with accurate edge detection, masking provides highly accurate results even with a single-view sinogram.

Theorem 1 does not tell us how to collect edge information from CT sinogram data (especially sparse-view data), however. Indeed this is the limiting factor in the method at this point. In this proof-of-concept paper, the edge map is obtained by transforming the FBP image, which is efficient to compute, into the edge domain by the sparsity transform $\mathbf{D}$, and thresholding. There are likely more accurate options. For example, it is possible to design a filter in the FBP process that will suppress most information besides edges. Once the edge map is obtained, the image is constructed by solving the noise-present variant of Eq. (2). The full method is given by Algorithm 1.

The newly introduced threshold parameter $\tau$ is of particular importance. As in Theorem 1, the success of the algorithm hinges on having an accurate edge map. Empirically, choosing a higher value $\tau$ may remove some smaller edges but overall yields better results than a less-restrictive edge map, with lower value $\tau$, that might include some false edges created by noise or artifacts from the limited number of views and the FBP process (see Figures 2 and 3 later). The regularization parameter, $\lambda$, on the other hand was empirically shown to be not as critical when edge-masking was used [2].

3 Results

As an example, an image is reconstructed from a 45-view sinogram. Figure 2 shows the edge mask creation process. In this example, the edge map to be used in Algorithm 1 only includes the outer edges of the phantom. Despite this edge mismatch, in Figure 3 the image is compared with FBP and Eq. (1) and still achieves the most accurate result. This demonstrates the robustness of the algorithm to perform well even when not all edges are present in the mask.

In these results the conjugate gradient algorithm is used to solve Eq. (5). To solve Eq. (1) the Split Bregman method [3], which requires multiple conjugate gradient iterations, is used. It was set to use $10$ iterations, such that the solve takes about $10$ times longer than Algorithm 1, but more could be necessary. Hence Algorithm 1 has a distinct advantage in terms of reconstruction speed.

4 Conclusion

This paper offers theoretical evidence to explain the potential importance of edge locations in CT image reconstruction. A practical scheme is then developed to utilize this theory. Finally, empirical evidence is presented for using edge-masked regularization to improve accuracy and speed in CT image reconstruction using limited data in order to reduce dose.

There is still much to explore. First and foremost, realizing that the optimal sparse feature for many CT images may not be physical edges, this algorithm is flexible and can accommodate any non-physical “edge”, where the definition is expanded to mean a nonzero location in a sparse domain. Hence we hope to explicitly explore other sparsifying transforms used in CT. In addition, we hope to test our method on more realistic CT forward models. Finally, we believe that because of its reconstruction speed our method has excellent potential for 3D reconstruction.