Abstract

We attempt to revitalize researchers' interest in algebraic reconstruction techniques (ART) by expanding their capabilities and demonstrating their potential in speeding up the process of MRI acquisition. Using a continuous-to-discrete model, we experimentally study the application of ART into MRI reconstruction which unifies previous nonuniform-fast-Fourier-transform- (NUFFT-) based and gridding-based approaches. Under the framework of ART, we advocate the use of nonlocal regularization techniques which are leveraged from our previous research on modeling photographic images. It is experimentally shown that nonlocal regularization ART (NR-ART) can often outperform their local counterparts in terms of both subjective and objective qualities of reconstructed images. On one real-world k-space data set, we find that nonlocal regularization can achieve satisfactory reconstruction from as few as one-third of samples. We also address an issue related to image reconstruction from real-world k-space data but overlooked in the open literature: the consistency of reconstructed images across different resolutions. A resolution-consistent extension of NR-ART is developed and shown to effectively suppress the artifacts arising from frequency extrapolation. Both source codes and experimental results of this work are made fully reproducible.

1. Two Cultures: Mentally versus Experimentally Reproducible Research

There are two cultures related to medical image reconstruction: theory oriented (i.e., mentally reproducible) and application oriented (i.e., experimentally reproducible). Historically, when Wilhelm Rontgen discovered X-rays, their applications into medical imaging was immediate (he took the very first picture of his wife’s hand using X-rays). By contrast, when Johann Radon studied integral geometry and measure theory in 1910s, he could not foresee the wide application of his celebrated Radon transform in tomography. The world had to wait until 1930s—when the radiologist Alessandro Vallebona first demonstrated the potential of radiography—to appreciate the usefulness of such a mathematical tool as the Radon transform. The field of modern tomography is largely founded on two schools of researchers: one deals with the exploitation of basic physical phenomenon like Röntgen (e.g., electromagnetic interaction and electron-positron annihilation) and the other studies the mathematical abstraction of imaging modality and its inverse (i.e., tomographic reconstruction) like Radon.

The invention of digital computers in 1940s has had a profound impact on the evolution of science and engineering including tomography. On the one hand, ever-increasing computing resources have dramatically facilitated the development of numerical algorithms (e.g., the invention of fast Fourier transform (FFT) in 1960s). On the other hand, computers have played an unexpected role in disconnecting the two schools: theorists (applied mathematicians) do not need to consult practicians (e.g., practical medical physicists and engineers) for real-world data anymore because they can use computers to simulate the whole forward imaging process including the nice-looking image, while practicians have become reluctant to accept “new” tools because those tools—despite their theoretical appeal—do not always lead to tangible benefits in practice (partially due to the mismatch between theoretic models and real-world data). For example, the traditional filtered back-projection (FBP) algorithm is still the choice by many practicians [1] regardless of the development of more powerful algebraic reconstruction techniques (ART) [24] by academia.

It is often argued that prohibitive computational complexity and memory requirements are the obstacles to the wide adoption of ART (e.g., in [4]). However, in view of the dramatic advances in computing technologies over the past decades, we argue that computational approaches such as ART and their extensions could strike a better tradeoff between the cost (e.g., in terms of acquisition time) and performance (e.g., as measured by subjective quality of reconstructed images). When compared with alternative physics or hardware-based approaches, computational methods are often cost effective and likely to play an important role in making the healthcare more affordable in the future. In this paper, we will take an experimentally reproducible approach toward shedding some novel insights into traditional ART and studying the impact of nonlocal regularization techniques on MRI image reconstruction. The moral of our story is, just like the application of wavelets into biomedical signal processing [5], recent advances in nonlocal image processing and nonconvex optimization could leverage into the community of biomedical imaging. In addition to their theoretical appeal, we hope that the reproducibility of our experiments could further stimulate this line of research and expedite its adoption by MRI practicians.

The rest of the paper is organized as follows. We first briefly review the current state-of-the-art of ART with an emphasis on its applicability to the so-called continuous-to-discrete models in Section 2. Then we introduce the class of nonlocal regularization techniques in Section 3 and demonstrate how nonlocal regularized ART (NR-ART) (Algorithm 1) could offer a computational approach to reduce MRI scan time in Section 4. The issue of resolution consistency with k-space data is addressed in Section 5, and we present a resolution-consistent extension of NR-ART algorithm. In Section 6, we draw some conclusions about our experimental studies and future research directions.

Input: measurement data , -space location ( ), and
density compensation function ;
Output: reconstructed image
(i) Initialization: set to be the output of TV-regularized
ART;
(ii) Main loop: for , set the threshold
(a) Projection onto the regularization constraint set:
where refers to a nonlocal regularized
filter in (3) or (4);
(b) Projection onto the observation constraint set:
where refers to the ART iteration (2);

2. ART with Local Regularization: Continuous-to-Discrete Models

An issue fundamental to the mathematical modeling of tomographic imaging (forward scanning) is how continuous spatial information is encoded into discrete measured data. Figure 1 includes a diagram of the so-called continuous-to-discrete model [9, 10]. Let stand for the targeted continuous function, for example, the continuous Fourier transform representation of a cross-section of human body. Discrete measurements from MRI scanners are often called k-space data. A key observation from the practice of MRI is that k-space data are often acquired on nonuniform sampling grids (refer to Figure 1(a)). It follows that the continuous-to-discrete model is beyond the reach of Nyquist-Shannon sampling theorem [11] dealing with uniformly-sampled data only. For the simplicity of notation, we use , where , to denote the continuous data acquired on the nonuniform sampling lattice in k-space and , where the discrete data on the uniform sampling lattices in the spatial domain (target of reconstruction). The problem of tomographic image reconstruction can then be stated as the estimation of discrete pixels ’s from continuous measurements ’s.

The art of modeling lies in the heuristics of approximation and their computational implications. Linearized forward models are often adopted for their computational tractability; that is, where denotes a term of additive noise and the linear operator can admit different approximations. For example, in MRI applications, the two approaches of approximations have been called Problem 1 (related to gridding [7]) and Problems 2 (related to nonuniform FFT [12]), respectively, in the literature [13, 14]. For linear models, both noniterative (e.g., convolution or filtered back-projection) and iterative (e.g., ART originated from Kaczmarz method [15] and conjugate-gradient methods [16]) reconstruction methods have been extensively studied in the literature. The iteration of ART goes like where is a relaxation parameter controlling the convergence behavior of ART iteration. A geometric interpretation of ART iterations can be found in [17], and it can be viewed as a generalization of Kaczmarz method [15] from Hilbert space to metric space ( is not self-adjoint any more [18] for the case of undersampled MR data). Under the assumption that can be approximated by a linear operator with spectral radius (maximum eigenvalues) of , it can be shown by eigenvalue analysis [19] that the convergence condition for (2) is . The class of conjugate-gradient methods can be seen closely related to (2) by replacing with a term arising from Gram-Schmidt orthonormalization [20].

Surprisingly, experimental findings related to ART seem scarce in the open literature. Standard textbooks such as [2123] mostly focus on the theory of ART while presenting minimal experimental results (e.g., only one figure related to ART was found in [21, p. 455]). It is only briefly mentioned in [22] that ART has not been widely used for MRI reconstruction. However, experimental studies of the closely-related problem of regularized restoration of photographic images are abundant in the literature (e.g., please refer to [19] for a review of state-of-the-art by 1990s and [24] for more recent advances). With rapid advances of computer simulation for medical image reconstruction, we deem it worthwhile to have some experimental study of how the standard total-variation- (TV-) based regularized ART works on MRI image reconstruction even for the tutorial purpose (please refer to Figures 2, 3, and 4). Our experimental studies are also expected to shed some novel insights on the limitation of local regularization, thus motivating the introduction of nonlocal regularization in the next section.

Simulation of forwarding MRI encoding is based on the image reconstruction toolbox (available at http://www.eecs.umich.edu/~fessler/code/) and SparseMRI toolbox (available at http://www.stanford.edu/~mlustig/SparseMRI.html). In our experiments, additive white Gaussian noise with zero mean and variance of is used to simulate . Two test images—a synthetic and a more realistic —are used in our experiments (refer to Figures 3(a) and 4(a)). For the k-space trajectory shown in Figure 1(a), reconstruction is based on 13392 samples which implies that matrix is sized by (therefore computing pseudoinverse directly is impractical); both implementations of operator and its adjoint are available from the image reconstruction toolbox. SparseMRI toolbox offers an efficient implementation of TV regularization, which will be used as the benchmark in the next section.

Figure 2 includes the comparison of peak sign-to-noise-ratio (PSNR) profiles of ART with and without TV regularization for the two test images, respectively. It can be observed that (1) iterative reconstruction such as ART does outperform zero-padding FFT (note that the starting point corresponds to zero-padding FFT). This is not surprising because it is well known that terminating ART after a finite number of iterations can be viewed as a strategy of regularization as well [25]; (2) the gap between the solid (w/o TV regularization) and dashed (with TV regularization) curves is larger on than , which can be explained away by noting that TV model better fits the class of images that are piecewise constant. However, as we inspect the visual quality of reconstructed images in Figures 3 and 4, we note that spiraling artifacts are the primary nuisance in computer simulation of continuous-to-discrete models. The severity of those artifacts depends on the parameter setting of k-space data simulation (e.g., the variable density spiral and its associated Gaussian density decay). Since those artifacts are not associated with the object of imaging, it is critical to suppress them via enforcing the prior constraint in ART. Apparently TV-based regularization is not sufficient to suppress the spiraling artifacts. What other tool can we use?

3. ART-Based Reconstruction with Nonlocal Regularization

In this paper, we consider a new class of nonlocal regularization techniques. Unlike previous variational formulation such as [26] which introduced the nonlocal term using level-set evolution of boundary curves, we propose to work directly with the self-similarity of important structures (e.g., translation invariance of edges, bilateral symmetry of objects) in medical images. The rationale is that the target functional could be implicitly defined by nonlocal similarity-based projection operators; for example, in our previous work [24], we have rigorously shown that such nonlocal projection operators are nonexpansive maps and therefore admit fixed points [27]. We argue that the main strength of such a fixed-point based approach is that it better fits an engineer’s intuition; for example, it is not always possible to analytically articulate the target functional especially when nonlinearity is involved but often more feasible to come up with a filtering solution (a projection operator) to suppress undesirable noise. Next, we will elaborate on the two specific examples of nonlinear projection operators for nonlocal regularization (note that, as the consequence of nonlocal regularization, the implicitly defined cost functional is inevitably nonconvex).

3.1. Nonlocal Perona-Malik Diffusion (PMD)

In our recent work [28], we proposed a nonlocal extension of PMD [29] as follows: where is a relaxation parameter, denotes the difference operator, is the nonlinear edge-stopping function, and the last term involves nonlocal diffusion. In the context of biomedical imaging, where the imaging object is often approximately bilaterally symmetric, we have adopted a choice of setting the nonlocal neighbor to be the pixel at the location of the mirrored position. Two types of edge-stopping function are suggested in the original PMD paper [29]: , and where is a constant controlling the tradeoff between backward and forward diffusions. We have empirically found that usually achieves better PSNR performance than and therefore adopted in our experiments.

3.2. Nonlocal Translation-Invariant Thresholding

In recent years, there has been a flurry of works on nonlocal image processing (e.g., nonlocal mean denoising [30], K-SVD denoising [31], block-matching 3D denoising [32], nonlocal total-variation restoration [33], learned simultaneous sparse coding [34], and clustering-based sparse representation [35]). It has become clear that nonlocal regularization techniques are capable of exploiting the global translational invariance property of important image structures such as edges and textures, which are complementary to the class of local regularization techniques. Several nonlocal thresholding strategies can be described as where denotes the nonlocal thresholding operator ( is a standard hard thresholding operator and represent the forward/inverse patch transforms as described in [31, 35]). We have empirically found that block-matching 3D (BM3D) thresholding [36] (the first step of BM3D denoising algorithm [32]) is computationally efficient due to its MEX-based implementation under MATLAB.

It is straightforward to incorporate nonlocal regularization into ART under the framework of alternating projections [37]. Observation data and prior knowledge, respectively, determine two constraint sets whose intersection defines the boundary of searching for the unknown target. In theory, projection onto convex sets would guarantee the convergence, while, in practice, it has been shown alternating projections could still produce useful results despite lack of convexity in some constraint sets (e.g., in the application of texture synthesis [38]). Here, we propose to experimentally study the performance of an iterative image reconstruction algorithm called nonlocal regularized ART (NR-ART). The flow chart of our NR-ART algorithm is described as follows. It should be noted that the NR-ART algorithm has some connection with the idea of augmented Lagrangian method as described in [39]. More specifically, our projection-based algorithm also admits a variational interpretation though the objective functional is likely to be nonconvex [40], which motivates our adoption of deterministic annealing strategy [41] in NR-ART. The role played by the decreasing thresholds is similar to that of an augmented Lagrangian multiplier [42]. Following the experimental setup in the previous section, we can validate the gain of nonlocal regularization on the two test images and compare it against that achieved by TV regularization [6]. Figure 5 compares the PSNR profiles of NR-ART algorithm with two choices of regularization operators: nonlocal PMD versus BM3D thresholding. We note that, despite the competing performance achieved by nonlocal PMD for the image, its PSNR performance is strikingly inferior to that of BM3D thresholding for the image. Again such discrepancy can be interpreted as a warning sign for designing regularization techniques; working with image might be misleading because computer-generated and real-world images often have highly different characteristics!

Figures 7 and 8 include comparisons of reconstructed images by various regularization techniques. It is encouraging to observe consistent gain achieved by nonlocal regularization over its local counterpart. For the image, visual inspection of subjective quality also strongly favors those obtained by NR-ART; for the image, subjective quality difference is less obvious due to visual masking. To facilitate visual inspection, we have compared the magnitude of reconstruction errors in Figure 6. The comparison of error magnitude maps more clearly shows that NR-ART is capable of more faithfully reconstructing the image than TV-regularized ART around regular edges and textures.

4. Nonlocal Regularized ART: From Simulated to Real-World Data

What is wrong with using simulated k-space data? To say the least, there is little known about how well software such as SaprseMRI could approximate the actual MRI process in the real world. The spiraling artifacts we have observed in previous sections only appear to be correlated with the type of noisy spike artifacts known in the literature of MRI [22]; taking other artifacts (e.g., ghost or data clipping) into account could have called for a more sophisticated investigation of computer simulation. Therefore, it seems inevitable to validate the strength of any MRI reconstruction algorithms on real-world k-space data even though the original (or ground truth) is not available.

For real-world k-space data, the inverse of (1) is often implemented by a so-called gridding algorithm [7]. Gridding algorithms have been extensively studied by both the medical imaging and magnetic resonance communities in recent decades. On the one hand, engineers have been working on various aspects of the gridding algorithm such as the selection of convolution function [43], computationally efficient implementations [44, 45], and its connection with nonuniform FFT [12]. On the other hand, MRI practicians have independently discovered the connection of gridding to the singular value decomposition [46], invented reverse gridding [47] (the adjoint operator of gridding interpolation), and explored its use with parallel MRI [4850] in which data from multiple, spatially distinctive coils are combined to provide additional spatial information to the reconstruction. Without exception, TV-based regularization has been adopted as the standard (e.g., [6, 51, 52]).

In this section, we want to demonstrate the potential benefit of nonlocal regularization in reducing the scanning time of MRI. Fast MRI scanning has been a hot topic in the MRI community, and various parallel imaging techniques (e.g., SMASH [53], SENSE [54]) have been developed. As an alternative to hardware-based approaches, we advocate a software-based approach here for its low cost and flexibility. Figure 9 shows the reconstructed image from one-third of the original k-space data (i.e., 4464 instead of 13392 locations). We use experimental results on one set of k-space data (downloaded from http://www.stanford.edu/class/ee369c/data/) to demonstrate the potential benefit of nonlocal regularization. Our implementation of adjoint gridding operator (i.e., reverse gridding in [47]) is an ad hoc extension of original gridding implementation (available at http://www.stanford.edu/class/ee369c/mfiles/). A set of manually chosen parameters are used for NR-ART algorithm: .

Figure 9 includes the comparison of reconstructed images from various competing techniques. Due to downsampling, gridding algorithm starts to fall apart and renders noticeable artifacts; the implementation available from SparseMRI toolbox does not produce satisfactory result due to its adoption of NUFFT-based implementation (In other words, we do not think this is the problem caused by TV-regularization but mismatch between simulated and real-world MRI models). By contrast, NR-ART is still capable of generating an image with acceptable quality from one-third of k-space samples. We note that the scan time saving offered by nonlocal regularization is complimentary to hardware-based parallel imaging strategies. Therefore, it would be highly desirable to validate the benefit of NR-ART for multicoil MRI data such as those acquired by SENSE [54].

5. Resolution-Consistent NR-ART: A Cross-Validation Approach

For a collection of real-world k-space data, a user can arbitrarily specify the spatial resolution of the reconstructed images because the continuous-to-discrete model can be viewed as a generalized analog-to-digital conversion process. It has been well known from the theory of Fourier imaging that spatial interpolation is equivalent to frequency extrapolation [10]. Specifically, the time-frequency duality implies that resolution enhancement in the spatial domain can be implemented by compressing the k-space data in the Fourier domain. However, there is no free lunch; as more k-space data are compressed into the range of lower-frequencies, the absence of higher frequencies often causes various artifacts as shown in Figures 11(a) and 11(b). It is natural to ask: how can we obtain resolution-consistent reconstruction results? A related issue of resolution and noise properties has been addressed for PET image reconstruction in [55], but the issue of resolution consistency does not appear to have been addressed in the open literature.

Figure 10(a) shows one way of achieving resolution consistency by a cross-validation approach. The basic idea is to alternate the projections of NR-ART between the low resolution (LR) and the high resolution (HR). Note that nonlocal regularization is enforced at both resolutions; moreover, the observation data at two different resolutions are connected by a pair of projection and back-projection operators. Such a projection point of view allows us to further generalize the continuous-to-discrete model in (1)—that is, two discrete models at varying resolutions can be connected by the bridge—the continuous-function . In the language of operators, we can construct a new operator and its adjoint to achieve a resolution-consistent NR-ART.

Figures 11(c) and 11(d) contains the LR and HR images reconstructed by the resolution-consistent NR-ART algorithm. It can be observed that artifacts at the HR associated with the original gridding reconstruction have been effectively suppressed by our cross-validation approach. One potential advantage of resolution-consistent NR-ART is that we can choose the scaling factor in such a way that the compressed k-space data in the Fourier domain could possess desirable properties (e.g., lying closer on average to the locations on an integral lattice, therefore producing smaller interpolation errors; please refer to Figure 10(b)). Another interesting and promising issue is related to the sharpness of reconstructed images. When most high-frequency components are absent, image blurring is also inevitable. Such a drawback is beyond the hope of linear reconstruction strategies; however, recent advances in blind image deconvolution have shown the promise of exploiting a priori knowledge about the blurring kernel and image source. Figure 12 includes two examples of blindly-deblurred images for Figures 9(c) and 11(d). Visual inspection appears to support their plausibility, but more objective evaluation is necessary to confirm that only the imaging object (not interfering artifacts) gets enhanced. We reserve this issue for future study.

6. Conclusions and Discussions

In this paper, we have experimentally studied the extensions of ART into MRI reconstruction: NR-ART for simulated k-space data and resolution-consistent NR-ART for real-world k-space data. The key take-home messages include (1) both NUFFT-based and gridding-based reconstructions can be unified into the framework of ART with the introduction of a continuous-to-discrete model; their adjoint operators admit a conjugate-gradient implementation of reversing the continuous-to-discrete model [47]; (2) nonlocal regularization techniques have the potential to outperform their local TV-based regularization because important image structures in medical images can be approximately characterized by translational invariance of local patches [32]; (3) enforcing resolution consistency offers a plausible approach for suppressing artifacts arising from frequency extrapolation; when artifacts are properly under control, the sharpness of reconstructed images can be further enhanced by newly-developed blind deconvolution techniques such as [8].

Several significant issues remain open. First, mathematical modeling of the MRI forward scanning process is the source of systematic errors causing varying interpretations of experimental results with simulated and real-world k-space data. More investigation about the continuous-to-discrete model (e.g., the approximation quality of operator and its adjoint) is necessary to sharpen our understanding about strengths and limitations with digital representations of analog imaging objects. We do need a new nonlinear sampling theorem for the continuous-to-discrete model (in contrast to the one established for discrete-to-discrete model in CS theory [56, 57]). Second, a priori knowledge in the context of human imaging is likely to significantly differ from that used to image the natural world. For example, it has been recently shown that only 8 radial lines in the Fourier domain (sampling density: ) are sufficient for a perfect reconstruction of Shepp-Logan phantom image [28]. Such experimental finding calls for deeper theoretical understanding of nonlocal sparsity (e.g., from translation invariance to transformation invariance). Third, despite visually compelling experimental results achieved using blind deconvolution techniques, their validity remains to be confirmed for more real-world k-space data.

To make sure this research is fully experimentally reproducible, we have established a homepage for this project at http://www.csee.wvu.edu/~xinl/demo/NRART.html where both MATLAB source codes and saved experimental results can be accessed. It is also our sincere hope that this work can further stimulate the medical imaging community’s interest in reproducible research [58]. A section devoted to sharing source codes and test data related to medical imaging has been created at http://www.csee.wvu.edu/~xinl/source.html and http://www.csee.wvu.edu/~xinl/database.html. It is this author’s opinion that the lack of test data is a major obstacle to the advance of medical imaging research. Unlike photographic images, the collection of real-world k-space data is seldom publicly available. Even though the privacy of patients is often a concern, one could argue that we can surely get around this issue by various modern technological means (e.g., deidentification and access control). If the two-culture phenomenon persists within the community of medical imaging researchers and practicians, the dearth of publicly available test data and reproducible research is likely to be a contributing factor.

Acknowledgment

This work was partially supported by Grant NSF-CCF-0914353.