PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Magn Reson Med. Author manuscript; available in PMC 2012 January 1.
Published in final edited form as:
PMCID: PMC3005027
NIHMSID: NIHMS220627

A Fast Edge Preserving Bayesian Reconstruction Method for Parallel Imaging Applications in Cardiac MRI

Abstract

Among recent parallel MR imaging reconstruction advances, a Bayesian method called Edge-preserving Parallel Imaging with GRAph cut Minimization (EPIGRAM) has been demonstrated to significantly improve signal to noise ratio (SNR) compared to conventional regularized sensitivity encoding (SENSE) method. However, EPIGRAM requires a large number of iterations in proportion to the number of intensity labels in the image, making it computationally expensive for high dynamic range images. The objective of this study is to develop a Fast EPIGRAM reconstruction based on the efficient binary jump move algorithm that provides a logarithmic reduction in reconstruction time while maintaining image quality. Preliminary in vivo validation of the proposed algorithm is presented for 2D cardiac cine MR imaging and 3D coronary MR angiography at acceleration factors of 2-4. Fast EPIGRAM was found to provide similar image quality to EPIGRAM and maintain the previously reported SNR improvement over regularized SENSE, while reducing EPIGRAM reconstruction time by 25-50 times.

Keywords: parallel imaging, edge preserving priors, Bayesian reconstruction, graph cuts, EPIGRAM

Introduction

Over the past decade parallel imaging techniques such as SENSE (1), SMASH (2) and GRAPPA (3) have been developed for rapid magnetic resonance (MR) imaging. These techniques employ multiple receiver coils to acquire under-sampled k-space data, which causes aliasing artifacts in images, and then use coil spatial sensitivities to remove (unfold) these artifacts in either k-space or image space. This paper focuses on the image space approach where the SENSE algorithm is considered the standard. However, SENSE image reconstruction is known to achieve unfolding at the cost of noise amplification, particularly when the coil sensitivity encoding system is ill-conditioned. To reduce noise and improve unfolding performance, a new algorithm called Edge Preserving parallel Imaging reconstruction with GRAph cuts Minimization (EPIGRAM) was proposed (4), in which the image reconstruction was formulated as a Bayesian estimation problem using an edge-preserving spatial prior. This formulation imposes piecewise smoothness as opposed to global smoothness or minimum norm in standard Tikhonov type regularization schemes (5,6). The Bayesian estimation problem is converted into an energy minimization problem and solved via iterative application of the alpha-expansion move algorithm (7). The method works by discretizing the grayscale intensity levels in an image into a set of integer-valued labels – e.g. 256 labels for 8-bit precision. We follow the same terminology in this work, and use the word “label” to mean grayscale level. The number of alpha-expansion moves grows linearly with the number of these labels, making it slow for images requiring a high dynamic intensity range (e.g., more than 256 intensity labels). Algorithms that iteratively apply graph mincut (7) to solve energy minimization problems in imaging are generally referred to as graph cut algorithms, irrespective of the particular move method they use. For instance, alpha-expansion moves (7), alpha-beta swap moves (7,8) and jump moves (8) are all possible moves on labels of an image; all these moves are ultimately solved via the mincut operation on a graph, and are therefore referred to as graph cut algorithms.

This paper introduces the binary jump move-based graph cut algorithm, first proposed for energy minimization in stereo vision literature (8), that can be used to speed up EPIGRAM reconstruction while maintaining image quality.

Methods

Binary Jump Move Algorithm

A different graph cut algorithm, using a new move scheme called jump moves, is now introduced to improve reconstruction time of EPIGRAM (see detailed mathematical definition and notation in the Appendix) while maintaining similar reconstruction quality. The proposed jump-move approach is based on the original work by Vexler et al (8) and provides logarithmic speed up by reducing the number of binary graph-cut steps. We refer to this algorithm as the binary jump move algorithm and the new parallel imaging reconstruction technique as Fast EPIGRAM.

Define £ as a set of all possible labels in the intensity range of an image. Define a new set of labels S = [-2m-1, -2m-22m-2, 2m-1] with m = log2 (|£|). Here the size of set S is O (log|£|) and it is trivially shown that any label in £ can be reached by a combination of labels d [set membership] S (called a d-jump). Given a jump d [set membership] S, the energy function in EPIGRAM (see Eq.[A1]) is converted into a series of jump-moves where each pixel with intensity xp can be labeled either xp or xp+d. The binary decision for all the pixels to minimize the energy is made by finding a min-cut on the graph (8). Using this new labeling the same process of binary energy minimization is repeated for each jump d [set membership] S. Each new labeling reduces (or preserves) the objective function shown in Eq.[A1], and the process stops when no d-jump can reduce it further. An illustrative example is shown in Fig.1. Note the number of proposed jump moves is logarithmic in the number of labels, compared to (4) where computational complexity is linear.

Figure 1
Example depicting the evolution of pixel labels in graph cut algorithms: a) starting labels, b) desired end labeling, c) evolution under binary jump move, d) evolution under alpha-expansion move. Bold labels are those that perform expansion moves, either ...

Imaging Experiments

Human experiments were approved by the local Institutional Review Board and written consent was obtained from all subjects prior to imaging. Five healthy volunteers (mean age of 27 years ± 3 standard deviation [SD]) were imaged using a 1.5 T GE Signa HD scanner (GE Healthcare Technologies, Waukesha, WI, USA). A standard eight-element phased array cardiac coil (4 anterior, 4 posterior elements) was used for signal reception. Subjects were imaged supine with vector electrocardiography gating for cardiac synchronization. Real-time gradient echo scout scans were performed to localize the diaphragm and the heart. Cine 2D balanced steady-state free precession (SSFP) cardiac short-axis images were acquired at the mid-ventricular level during a 16-heartbeat end-expiratory breath-hold. Imaging parameters were: TR = 4.3 ms, TE = 1.9 ms (partial echo), flip angle = 60°, readout bandwidth = ±125 kHz, FOV = 33 cm, phase FOV factor = 0.85, slice thickness = 8 mm, acquisition matrix = 160 × 218, views per segment = 16. The right coronary artery (RCA) was also imaged using an ECG-triggered segmented k-space navigator gated 3D SSFP coronary MR angiography (CMRA) pulse sequence (9) and the following imaging parameters: TR = 4.0 ms, TE = 1.2 ms (partial echo), flip angle = 60 degrees, readout bandwidth = ±62.5 kHz, FOV = 26 cm, slice thickness = 3 mm, number of slices = 16, acquisition matrix = 160 × 224, number of echoes per heartbeat = 32 (corresponding to an acquisition window of 128 ms). To monitor respiratory motion, a pencil-beam navigator was placed through the dome of the right hemi-diaphragm. The optimal delay time between cardiac trigger and diastasis was determined from a cine scout scan.

Parallel Imaging Reconstruction

Aliased images were obtained by under-sampling full-resolution image data along the phase encoding direction in Cartesian 2D k-space. The center 32 lines of k-space data of each coil were retained to estimate low-frequency coil sensitivity maps (self-calibration). These lines were first multiplied with a Hanning window to reduce ringing artifact, zero-padded to full resolution, and then Fourier transformed to image space. Relative coil sensitivities were estimated by dividing these images by the full resolution root sum of squares (RSOS) image (computed as the root sum of squares of intensity across all the coils for each pixel location). Model parameters for EPIGRAM and Fast EPIGRAM reconstructions – number of labels |£| = 256, K = |£| / 7 and λ = 0.08 |£| – were kept constant for all cases. Regularized SENSE reconstruction (4,10) was performed with several different values of the regularization parameter μ [set membership] [0,1] and the best regularized result, corresponding to the minimum root mean squared error (RMSE) with respect to the reference image, was obtained for μ = 0.1. This is the same regularization parameter selection criterion reported in (4). RMSE was calculated for both Fast EPIGRAM and regularized SENSE images for comparison. All software was implemented in C++ and run on a workstation with a 3.6 GHz Intel Xeon processor and 1 GB of RAM.

Relative Signal to Noise Ratio (SNR)

Blood SNR was calculated for the reconstructed 2D cine short-axis SSFP images of the heart and 3D CMRA images of the RCA (both data sets were under-sampled at R = 3) using the relative SNR measure for parallel imaging methods proposed by Reeder et al (11). Under this method, a single under-sampled k-space data set was replicated twice, each time adding separate records of synthetic complex Gaussian noise which was scaled and correlated according to the coupling of the coil elements. The noise covariance matrix across coil elements was estimated from a noise-only prescan (12) in one subject and was found to be well approximated by a diagonal matrix (the average off-diagonal entry was less than 10% of the average diagonal entry). By comparing these two noisy replicates it was possible to obtain estimates of relative SNR of the reconstruction. For SNR measurement, blood signal was measured from the LV cavity of a diastolic image for 2D cine imaging and from the aortic root for 3D CMRA using region of interest (ROI) analysis. Two-tailed paired-sample t-test was used to assess the difference in SNR between reconstruction techniques.

Results

Figure 2 shows reconstruction results of the 2D cardiac cine data with acceleration factor R = 2 obtained after 5 iterations over [32, 64, 128, 256, 512] labels in the intensity range, demonstrating very similar image quality and RMSE (relative to the reference image) provided by EPIGRAM and Fast EPIGRAM reconstructions. However, the reconstruction time of EPIGRAM, which uses the alpha expansion move algorithm, was much higher than that of the Fast EPIGRAM using the binary jump move algorithm (Fig.3). In the case of 512 labels and 5 iterations, the reconstruction time was more than 20 minutes for EPIGRAM compared to only 22 seconds for Fast EPIGRAM at acceleration factor R = 4 (approximately 50-fold reduction).

Figure 2
2D cardiac cine parallel imaging reconstructions obtained with EPIGRAM (which uses alpha expansion moves) and Fast EPIGRAM (which uses binary jump moves) at acceleration factor R = 2 after 5 iterations over [32, 64, 128, 256, 512] labels in the intensity ...
Figure 3
Comparison of reconstruction times of EPIGRAM (which uses alpha expansion moves) and Fast EPIGRAM (binary jump move algorithm) for 2D cardiac cine data after 5 iterations at acceleration factor R = 4. Plots are on a log-log scale, in order to better show ...

Over five subjects, Fast EPIGRAM achieved almost 100% SNR improvement (36 ± 13 for regularized SENSE (μ = 0.1) vs. 72 ± 16 for Fast EPIGRAM, p<0.001) for 2D cine images and 50% SNR improvement (42 ± 9 for regularized SENSE vs. 61 ± 12 for Fast EPIGRAM, p<0.001) for 3D CMRA images at acceleration factor R = 3. On average Fast EPIGRAM provided 20% lower RMSE than regularized SENSE. Table 1 summarizes Fast EPIGRAM reconstruction times for various image data matrices and accelerations.

Table 1
Fast EPIGRAM reconstruction times (using 256 intensity labels) for various image data matrices and acceleration factors R.

Discussion

Our preliminary data show that the proposed Fast EPIGRAM technique reduced reconstruction time by 25-50 times compared to EPIGRAM for 2D cardiac cine and 3D CMRA imaging. At acceleration factor of 2-4 this technique demonstrates significant improvement in SNR and RMSE over regularized SENSE, while effectively removing residual aliasing artifacts associated with SENSE. The logarithmic jump moves of Fast EPIGRAM reconstruction were the key to speed improvement. The jump move scheme takes steps that are powers of two. This is equivalent to bit flipping: if a pixel jumps by 2m then it is equivalent to flipping the mth bit from either 0 to 1 or 1 to 0. In contrast, EPIGRAM performs one expansion move per label, which makes it slow for images with high dynamic intensity range. Note that an iteration of alpha expansion move and an iteration of jump move do not produce the same result, which is why we compared the end performance (Fig.2) of both methods instead of results per iteration.

While our results and those in (4) appear to support the advantages of edge preserving priors in Bayesian parallel MRI, further work is needed to rigorously substantiate these priors. The extension of our method to sampling trajectories other than Cartesian is possible but non-trivial. Optimal selection criteria are needed for heuristically chosen parameters K, λ.

In this study, an identity weighting matrix A was chosen for the regularization term similarly to that used in recent works (4,10,13). A potentially better choice for A would be a diagonal matrix whose diagonal elements are inversely proportional to the local image SNR. While such SNR-based regularization (similar to Wiener filtering) can be a powerful approach when reliable a priori image information is available and may provide improved results, its feasibility for clinical imaging remains to be investigated in our future work.

In order to reduce reconstruction time, a 4-connected neighborhood instead of 8-connected neighborhood as in (4) was used in this study. According to our prior experience with various neighborhood size (unpublished data), this would lead to a 2-fold increase in processing speed without a discernible loss of reconstruction quality. Due to the long-range effect of Markov fields, similar performance is frequently observed under somewhat different neighborhood systems (14).

Although the developed Fast EPIGRAM reconstruction represents a critical improvement in speed over EPIGRAM and a major step towards clinical acceptance, it is still inadequate for our critical goal of real-time reconstruction (Table 1). To further reduce latency, we seek to achieve further time reduction through the use of parallel processing hardware architectures. In our estimation, a full volumetric reconstruction of 3D CMRA data may be possible in under 2 minutes with an 80-node computing cluster, by processing different image slices on different processor nodes in parallel. Further reduction in computational load may also be achieved by reconstructing an image of known phase (12), thereby replacing the original complex variables by real ones. Although this technique can potentially halve the reconstruction time, it is limited by the accuracy of estimated phase. Other algorithmic modifications that may reduce processing time include the use of fewer and non-uniformly distributed labels, multi-resolution approach where each iteration is initialized by a lower-resolution image from the previous iteration, and new energy minimization algorithms such as Belief Propagation (13).

In conclusion, Fast EPIGRAM provides significantly faster reconstruction than EPIGRAM as well as higher SNR than regularized SENSE for 2D cardiac cine imaging and 3D coronary MRA.

Appendix

Details of Epigram Algorithm

In EPIGRAM (4), the maximum a posteriori (MAP) Bayesian parallel reconstruction problem was recast as the following minimization:

X^MAP=argminX(YEX2+GEP(X))
[A1]

where vector X = {xp|p [set membership] P} contains the intensity values of the pixels of the target image, p indexes the set P of all pixels in the image, and N is cardinality of P. Matrix E encodes coil sensitivity responses and is a NC/R × N block-diagonal matrix (4), where R is the under-sampling or acceleration factor and C is the number of coils. Y is a vector of length NC/R, and contains the aliased images “seen” by the receiver coils. The prior cost function GEP(x) is an edge-preserving Markov Random Field prior (14,15) given as

GEP(X)=(p,q)Nsλmin(|xpxq|,K)
[A2]

which is the separation cost to assign intensities xp and xq to neighboring pixels p and q, and K, λ are user-defined constants. Ns is the spatial neighborhood system that consists of pairs of adjacent pixels, usually the 4-connected neighbors. This choice of priors is considered edge-preserving because it limits the maximum penalty for two neighboring pixels that have very different intensities, a condition that occurs at object boundaries. The resulting non-convex energy function was solved in (4) using graph cuts (7,16). We refer readers to (4) for a detailed description of the graph construction and to (17) for some interesting theoretical bounds in support of its strong performance. An expansion of Eq.[A1] into pixel-wise energy terms is described here for pixel intensities of complex value.

The computational problem in Eq.[A1] is to efficiently minimize

ɛ(x)=YEX2+(p,q)NsV(xp,xq)
[A3]

The first term in the right side of Eq.[A3] can be written as

YEX2=YHY+2real(YHEX)+(EX)HEX
[A4]

YHY is independent of X and can be ignored in minimizing Eq.[A4] over X. The output of each coil consists of only N/R pixels in the phase encoding direction where N is the total number of pixels in the phase encoding direction in the fully sampled image. Define lth coil output Yl([p with macron]) at pixel [p with macron], where [p with macron] maps to A[p with macron](ī,j) = {p(i, j): ī = mod(i,N/R), i [set membership] [1,N]} in the unfolded image that alias to each other due to under-sampling in the phase encoding direction. The phase encoding direction is indexed by i and the frequency encoding direction by j. Letting the sensitivity of coil l at pixel p be sl (p), the second term in Eq.[A4] can be written as

YHEX=p[lyl(p¯)Hsl(p)]xp
[A5]

Each pixel takes a complex value, which means that two real values (xrp,, xcp) must be resolved for each pixel, where the superscripts r and c denote the real and imaginary component, respectively. Therefore,

real(YHEX)=p[(l(ylr(p¯)slr(p)+ylc(p¯)slc(p)))xpr+(l(ylc(p¯)slr(p)ylr(p¯)slc(p)))xpc]
[A6]

The last term of Eq.[A4] can be written as

(EX)HEX=p(lsl(p)Hsl(p))xpHxp+(p,p)Na(lsl(p)Hsl(p))xpHxp
[A7]

where the set Na = {(p,p′): p [set membership] A[p with macron], p[set membership] A[p with macron], pp′} is defined over all pairs of aliasing pixels. This can be further factorized as follows:

(EX)HEX=A+B+CA=p(l(slr2(p)+slc2(p)))(xpr2+xpc2)B=(p,p)Na(l(slr(p)+slr(p)+slc(p)slc(p)))(xprxpr+xpcxpc)C=(p,p)Na(l(slc(p)slr(p)slr(p)slc(p)))(xprxpcxpcxpr)
[A8]

Grouping of all the single and pair-wise interaction terms in Eq.[A3] results in

ε(x)=a2+pbr(p)xpr22pcr(p)xpr+pbc(p)xpc22pcc(p)xpc+2(p,p)Nad1(p,p)(xprxpr+xpcxpc)+2(p,p)Nad2(p,p)(xprxpcxpcxpr)+(p,q)NsV(xp,xq)
[A9]

for appropriately determined functions br(p), bc(p), cr(p), cc(p), d1(p,p′) and d2(p,p′). Also, the prior function is extended for real and imaginary components:

V(xp,xq)=λ(min(|xprxqr|,K)+min(|xpcxqc|,K))
[A10]

Substituting Eq.[A10] in Eq.[A9] gives the energy function that is to be minimized for edge-preserving parallel imaging reconstruction.

References

1. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42(5):952–962. [PubMed]
2. Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38(4):591–603. [PubMed]
3. Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (GRAPPA) Magn Reson Med. 2002;47(6):1202–1210. [PubMed]
4. Raj A, Singh G, Zabih R, Kressler B, Wang Y, Schuff N, Weiner M. Bayesian parallel imaging with edge-preserving priors. Magn Reson Med. 2007;57(1):8–21. [PMC free article] [PubMed]
5. Ying L, Xu D, Liang ZP. On Tikhonov regularization for image reconstruction in parallel MRI. Proceedings of IEEE EMBS 26thAnnual International Conference; 2004. pp. 1056–1059. [PubMed]
6. Liang ZP, Bammer R, Ji J, Pelc N, Glover G. Making better SENSE: Wavelet denoising, Tikhonov regularization and total least squares. Proceedings of ISMRM 10th Scientific Meeting; 2002. p. 2388.
7. Boykov Y, Veksler O, Zabih R. Fast approximation energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2001;23(11):1222–1239.
8. Veksler O. PhD dissertation. Cornell University; 1999. Efficient graph-based energy minimization methods in computer vision.
9. Nguyen TD, Spincemaille P, Wang Y. Improved magnetization preparation for navigator steady-state free precession 3D coronary MR angiography. Magn Reson Med. 2004;51(6):1297–1300. [PubMed]
10. Lin FH, Kwong KK, Belliveau JW, Wald LL. Parallel imaging reconstruction using automatic regularization. Magn Reson Med. 2004;51(3):559–567. [PubMed]
11. Reeder SB, Wintersperger BJ, Dietrich O, Lanz T, Greiser A, Reiser MF, Glazer GM, Schoenberg SO. Practical approaches to the evaluation of signal-to-noise ratio performance with parallel imaging: application with cardiac imaging and a 32-channel cardiac coil. Magn Reson Med. 2005;54(3):748–754. [PubMed]
12. Kellman P, McVeigh ER. Image reconstruction in SNR units: a general method for SNR measurement. Magn Reson Med. 2005;54(6):1439–1447. [PMC free article] [PubMed]
13. Liu B, King K, Steckner M, Xie J, Sheng J, Ying L. Regularized sensitivity encoding (SENSE) reconstruction using Bregman iterations. Magn Reson Med. 2009;61(1):145–152. [PubMed]
14. Li SZ. Markov random field modeling in computer vision. London: Springer-Verlag; 1995.
15. Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1984;6(6):721–741. [PubMed]
16. Kolmogorov V, Zabih R. What energy functions can be minimized via graph cuts? IEEE Transactions on Pattern Analysis and Machine Intelligence. 2004;26(2):147–159. [PubMed]
17. Raj A, Singh G, Zabih R. MRF's for MRI's: Bayesian reconstruction of MR images via graph cuts. Proceedings of IEEE Conference on Computer Vision and Pattern Recognition. 2006;v1:1061–1068.