Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Magn Reson Imaging. Author manuscript; available in PMC 2012 December 1.
Published in final edited form as:
PMCID: PMC3324824

3D Late Gadolinium Enhancement imaging of the left atrium with stack of stars and compressed sensing



To develop and test a hybrid radial (stack of stars) acquisition and compressed sensing reconstruction for efficient Late Gadolinium Enhancement (LGE) imaging of the left atrium.

Materials and Methods

Two hybrid radial acquisition schemes, kx-ky-first and kz-first, are tested using the signal equation for an inversion recovery sequence with simulated data. Undersampled data reconstructions are then performed using a compressed sensing approach with a 3D total variation constraint. The data acquisition and reconstruction framework is tested on five atrial fibrillation patients after treatment by radio-frequency ablation. The hybrid radial data is acquired with free breathing without respiratory navigation.


The kz-first radial acquisition gave improved image quality as compared to a kx-ky-first scheme. Compressed sensing reconstructions improved the overall quality of undersampled radial LGE images. An image quality metric that takes into account the signal, noise, artifact and blur for the radial images was 35%(±17%) higher than the corresponding Cartesian acquisitions. Total acquisition time for 36 slices with 1.25 × 1.25 × 2.5 mm3 resolution was under three minutes for the proposed scheme.


Hybrid radial LGE imaging of the LA with compressed sensing is a promising approach for obtaining images efficiently and offers more robust image quality than Cartesian acquisitions that were acquired without a respiratory navigator signal.

Keywords: RF ablation, left atrium, LGE imaging, stack of stars, compressed sensing


Atrial fibrillation is a growing health concern affecting over 7 million people in the US and Europe. Radio frequency (RF) ablation of the pulmonary vein ostia and the left atrium (LA) is a promising treatment option for restoring sinus rhythm. Late Gadolinium Enhancement MR imaging (LGE-MRI) of the LA has shown the potential to assist in various aspects of assessing and treating atrial fibrillation. By measuring signal enhancement in the LA post ablation, injured tissue can be identified and quantified (12). Good correlations were reported between regions of late enhancement and the sites of ablation when imaging was performed 3 months post procedure (3). The LGE information can be used to identify regions where the tissue has regenerated or was incompletely ablated. LGE imaging can also be a tool to triage patients for the RF procedure (45). It has been reported that the extent of enhancement from pre-existing fibrosis in the LA is correlated with poor outcome of the ablation procedure (4).

Currently, LGE images of the LA are obtained using a 3D Cartesian acquisition with ECG gating and a respiratory navigator (12). In every heartbeat, after the inversion pulse, segments of 3D k-space are acquired in the diastolic phase of the cardiac cycle. Acquisition of the same set of segments is repeated until it is acquired when the diaphragm position is within a specified navigator window. Parallel imaging techniques like GRAPPA can be used to reduce total scan time (2). However in a standard clinical setup with coils that are not specifically designed for imaging the atrium, the acquisition can be slow with acquisition times typically on the order of 8 minutes. Respiratory navigation using signal from the diaphragm assumes that the diaphragm and the LA positions are directly related which may not be always true. As well, the images may not be useful if there is navigator drift or poor TI selection. More rapid acquisitions would allow for repeat exams if data quality was an issue. And acquisitions that were not reliant on diaphragm position would be very useful. Here we explore a hybrid radial acquisition (radial in-plane and Cartesian in the slice direction), also termed stack of stars, that is more robust to motion (6) and to undersampling. In conjunction, we use compressed sensing reconstruction with a total variation constraint (78) to remove artifacts from radial undersampling of the kx-ky planes. Unlike parallel imaging techniques (910) that use coil sensitivities for accelerating acquisitions the compressed sensing approach exploits the sparsity of a transform of the images.


Data acquisition theory

Simulation studies were performed to determine the optimal data acquisition strategy. A 3D digital phantom was generated from a patient dataset by manually segmenting 17 different anatomical structures from 36 axial slices. Figure 1 shows the segmentations on a single slice and a 3D visualization across all of the slices. The enhancement in the posterior side of the atrial wall simulates scarred or fibrotic tissue post ablation. These structures were thresholded to generate a piecewise constant phantom. Realistic T1 and M0 values (1112) as shown in Table 1 were assigned, although T1 values can vary depending on the time delay after Gd injection and on the field strength. Then the signal evolution at each readout for each of the regions was generated using Equation [1] for an inversion recovery sequence (13).

equation M1

In the above equation, TD is the time between the last α pulse and the next inversion pulse, TI is the time between the inversion pulse and the first α pulse, and N is the number of k-space views acquired after each inversion pulse. Here we chose TR=6.8 msec, N=36, flip angle α = 22° and TI=280 msec to null the myocardium when sampling the center of k-space.

Figure 1
Generation of simulated data. (a) Single slice from a fully sampled 3D acquisition with manually segmented regions overlaid in color. (b) Color coded regions segmented from the entire 3D volume. (c) Single slice from a 3D phantom generated by using different ...
Table 1
T1 and M0 values of the segmented structures used for simulations (1112). 'Tissue' corresponds to tissue regions of left ventricle, right ventricle, left atrium, right atrium; 'Blood' corresponds to blood pools of descending aorta (DA), ascending ...

Figure 2 shows a schematic of the k-space data acquisition. TI corresponds to the null point of the signal from myocardium. Different segments of k-space acquired in the diastolic phase of the cardiac cycle are color-coded and are interleaved. Here we simulate two hybrid radial schemes (i) kx-ky-first and (ii) kz-first and compare them to the Cartesian approach (14). In the kx-ky first scheme, radial lines in the kx-ky plane are all acquired first before incrementing the Cartesian kz-encoding. In the kz-first scheme, radial lines are acquired first along all of the kz encodes before acquiring the next ray in kx-ky. The difference between the two hybrid schemes is that in the kx-ky first method, the rays are consistent along the kz dimension but not in kx-ky; as contrasted to the kz-first scheme rays in each kx-ky plane are consistent but there is inconsistency across the kz encodes.

Figure 2
(a) Schematic of a segmented respiratory navigated 3D Cartesian acquisition scheme. Different segments of k-space are color coded. (b) 3D k-space schematic of hybrid radial kx-ky first (for each kz partition all rays in-plane are acquired first). (c) ...

3D phantom k-space data for both schemes was generated by sampling the 3D k-space generated at each α pulse using approximate radial lines. For both the schemes we sampled 144 radial lines in each kx-ky plane, each with 512 readout samples, and 36 kz encodes to simulate the undersampling. In contrast to the standard hybrid sampling scheme in which the same set of radial lines in each kx-ky plane are acquired for different kz encodings, here each plane is rotated by an angle based on golden ratio (15) in order to create incoherent artifacts along the slice dimension as well. A similar scheme but with an angle spacing that is half of the angle spacing between two in-plane rays was used to obtain differing streaking pattern across the slices for MR angiography data (16). The golden ratio based spacing scheme also gives flexibility in terms of acquiring an arbitrary number of kz encodings without predetermining the optimal offset angles between them. For comparison, a Cartesian acquisition was simulated with a variable density sampling scheme using 144 phase encodes, each with 512 readout samples, in each kx-ky plane. Nine lines around the center of k-space were contiguously sampled and the remaining lines were discarded in a random fashion in the outer region of k-space for each kz encoding.

Data acquisition in patients

The hybrid radial kz-first acquisition was implemented on a Siemens 3T Verio scanner and five atrial fibrillation patients were imaged 3–9 months post-ablation. The 3D hybrid radial as well as a 3D Cartesian acquisition were used. Pulse sequence parameters for the radial acquisition were TR=3.8 ms, TE=2.1 ms, acquisition matrix=512 (256 after removing oversampling) × 144 with 36 slice encodings, slice thickness=2.5 mm, FoV=320×320 mm2,TI=300–400 ms, flip angle = 14°. A Gd dose of 0.1 mmol/kg was used. All the parameters for the Cartesian acquisition were identical except the acquisition matrix was 512 (256 after removing oversampling) × 256 with 36 slice encodings. For the hybrid radial scheme this corresponds to a reduction (acceleration) factor of 5.5 according to the Nyquist sampling rate (due to oversampling by a factor of two in the read dimension) and a reduction factor of 1 for the Cartesian case. The radial acquisition corresponds to an acceleration factor of 1.7 when compared with the Cartesian acquisition. No respiratory navigators were used for the Cartesian as well as the radial acquisitions. TI scout images were acquired before the first LGE acquisition to determine TI for nulling the myocardium and for the subsequent acquisition it was empirically adjusted to be slightly longer. We randomized the order of acquisition of the hybrid radial and the Cartesian scans. Patient consent was obtained in accordance with regulations of the Institutional Review Board.

Undersampled data reconstruction

In order to remove the artifacts from radial undersampling we used a compressed sensing approach (1718). By exploiting sparsity in the spatial gradient of the images, reconstruction was performed by iteratively minimizing a cost functional consisting of an L1 norm constraint and a data consistency constraint as shown in Equation [2]. In the data consistency term, E represents the encoding matrix that transforms the image estimate m to match the radially sampled data d. γ is a regularization parameter that controls the tradeoff between the data consistency and the L1 norm constraint term.

equation M2

The L1 norm of the gradient transform term corresponds to the total variation (TV) constraint (19), in which Dx, Dy and Dz represent finite difference/gradient operators along x, y and z dimensions respectively. While the commonly used 2D TV constraint (78) can be applied on individual slices, improved reconstructions are possible with the 3D constraint as this takes into account sparsity in all three dimensions. An iterative gradient descent scheme was used for the minimization process (20). γ was chosen using the L-curve method (2122). The initial estimate for the simulated phantom was the image obtained using the inverse Fourier Transform of the simulated undersampled data as the data were sampled in a pseudo-radial pattern on a Cartesian grid.

Reconstruction for the in-vivo acquisitions with compressed sensing was done as outlined above. γ was chosen using the L-curve method for one radial dataset and for one Cartesian dataset and the same γ was then used for all other corresponding datasets. For the Cartesian cases γ was 0.005 while that for the radial cases was 0.025. Individual coil images were reconstructed and then combined in a square root of sum of squares fashion. For the radial acquisitions k-space data were pre-interpolated before the iterative compressed sensing reconstruction onto a Cartesian grid by identifying three sampled points nearest to a specified Cartesian point using Delaunay triangulation of the data and linear interpolation (2324). If more than one point is at the same nearest distance they are averaged before linear interpolation. During the minimization of the cost functional, Fast Fourier Transforms with no density compensation functions were used instead of convolution gridding or NUFFT (2526) leading to significantly faster reconstructions.

Evaluation of images

Images obtained from the proposed framework were compared to those obtained from the Cartesian scheme both qualitatively and quantitatively. Images were visually assessed for overall quality in terms of sharpness and ability to delineate fine structures. A quantitative quality metric was computed in two steps. In the first step a ratio of the mean signal intensity of a uniform region with no structures in two adjacent slices and the standard deviation of difference of signal intensities computed from the same region was determined. This is similar to SNR but the “noise” also reflects motion effects such as streaking or ghosting that can increase the non-uniformity in the signal. In the second step we use the blur metric proposed in (27) in order to account for blurring of edges due to motion and compressed sensing reconstructions. This method is based on the assumption that for a given image its variation with respect to a low pass filtered version is proportional to its sharpness. A sharp image would result in large variations between itself and its filtered version and a blurry image would result in small variations. Blur quantification for a given image is done by applying an averaging low pass filter and normalizing the difference. This blur metric, which ranges between 0 and 1 (larger=blurrier), is shown to correlate well with subjective evaluation on natural images (27).

Here we limit the computation of the blur metric to a square region of interest around the left atrium. The final quality metric was computed as the ratio of the SNR type of metric from the first step to the blur metric in the second step. This final quality metric represents a first order quantitative image quality measure (larger=better) that can be used for comparing two different images and which is sensitive not only to image artifacts and noise but also blurring.

For all of the in vivo datasets, quantification of the enhancement was done as previously reported in (2). The atrial wall was manually segmented and then the signal intensities were thresholded to give an estimate of the percent of the left atrium that was scar.



Figure 3 shows results of kx-ky first and kz-first hybrid schemes with compressed sensing reconstructions on simulated phantom data. The image corresponding to the kz-first scheme outperforms the kx-ky first scheme in terms of overall quality and preserving image contrast. The poor quality of the kx-ky first scheme is because k-space data are acquired while signal intensities in the myocardium and different regions are changing from negative to positive making it significantly inconsistent, even with interleaving. For both cases TI is defined as the time between the inversion pulse and the α pulse for reading the center of k-space (null point of the myocardium). In the kx-ky-first case the null point is defined here to be when the 19th radial line (of total 36 rays) is acquired in a heartbeat. Whereas with a centric kz-first scheme the null point is defined as the kz=0 encode, which is acquired first in each heartbeat. Thus TI is essentially longer for the kz-first method so there are fewer negative values than with the kx-ky acquisition. Also the inconsistency in this case is across the kz spatial frequencies and not within the kx-ky plane.

Figure 3
Simulation results. (a) Compressed sensing reconstruction of kx-ky first scheme with 144 rays. (b) Corresponding reconstruction of a kz-first scheme with labeled left atrium (LA). Arrow points to enhanced signal region in the atrium after ablation. (c) ...

While the Cartesian image is better than the kx-ky first scheme it is inferior to the kz-first scheme. Fig. 3c shows the residual ghosting and blurring obtained with the Cartesian acquisition due to undersampling. This is also illustrated in the difference images shown in Figs 3d and 3e. RMS error for the Cartesian acquisition was 72% higher than that for the radial acquisition. The image quality metric for the radially-sampled images was 21% higher than the Cartesian-sampled images. Coronal cuts shown in Figures 3f–3i show that while the kz-first image has blurring it is closer to the original image as compared to kx-ky-first and to the Cartesian images with compressed sensing.

Patient acquisitions

Figure 4 shows the results obtained with the Cartesian and radial acquisitions with total variation compressed sensing reconstruction. Two axial slices from two patients are shown. The phase encoding direction for the Cartesian images in the top row was chosen in the left-right direction and that for the bottom row was in the anterior-posterior direction. Using compressed sensing reconstruction on fully sampled Cartesian data can be equivalent to performing a total variation denoising. Cartesian images overall have poor quality as compared to radial for both phase encoding directions even though they took 1.7 times as long as the radial acquisitions. The images are blurred and/or have ghosting artifacts while the radial images are sharper with better delineation of different structures. The image quality metric for radial slices was 28% and 58% higher than the Cartesian images respectively for the top and the bottom row images. A mean overall improvement of 35% (±17) was observed in all five patients. For the Cartesian datasets an average of 15% (±5) of the left atrium was classified as scar while for the radial datasets it was 18% (±7).

Figure 4
In-vivo results from two patients (top and bottom rows). (a, c) Two axial slices from Cartesian acquisitions. (b, d) Closely matching slices from hybrid radial acquisitions. Phase encoding direction for the top row Cartesian images is left-right and that ...


Here we implemented a kz-first hybrid radial acquisition, which results in a smooth modulation of the non steady state inversion recovery signal in 3D k-space. This results in improved overall image quality as compared to a kx-ky first scheme. A weighted view sharing method of the kx-ky data as proposed in (28) within each plane could possibly improve kx-ky first results. However, it was reported that such a method for 2D LGE imaging of the left ventricle had poor quality as compared to a Cartesian acquisition (28). 3D imaging with kz-first can overcome some of these limitations. While we present results here for imaging the LA, the 3D hybrid radial and reconstruction method can be applied to LGE imaging of the left ventricle.

Respiratory navigation

In this study, we acquired radial data with ECG gating in the diastolic phase of cardiac cycle without respiratory navigator gating. This resulted in improved image quality and faster scan times as compared to a Cartesian acquisition scheme without respiratory navigation. Using a respiratory navigator can help reduce motion artifacts however at the cost of increased scan time. Figure 5 compares images from a Cartesian acquisition with respiratory navigation (acquired with a GRAPPA acceleration factor of 1.7 and a five channel cardiac coil) and the hybrid radial non-navigated scheme on a patient. While the image quality is slightly superior with the navigated scheme in terms of overall sharpness and contrast (due to the selection of appropriate TI), left atrial scar in both the images is visible (arrows). Quantification of the percent of scar in the left atrium for the radial acquisition was 28% and that for the Cartesian acquisition was 22%. Total acquisition time for the hybrid radial scheme was faster by about four times than the navigated approach. As well, the adaptive nature of the acceptance window can sometimes result in poor navigated image quality as the heart/LA might not have moved even though the diaphragm has drifted.

Figure 5
(a,c) Images acquired with respiratory navigated Cartesian acquisition with GRAPPA reconstruction. Arrows points to bright signal in the ablated atrium. (b,d) Closely matching slices from the hybrid radial acquisition with compressed sensing reconstruction. ...

Image quality assessment

Here we proposed a quantitative image quality metric that can gauge image artifacts, noise and blurring and can be used when there is no ground truth available. This metric however does not substitute for the human visual system and should be used in conjunction with expert visual inspection for diagnosis. The blur metric used here does not use the edge information in the images as this can vary depending on the threshold chosen for the gradient image (29). More sophisticated methods can also be used to quantify image quality as proposed in (30).

In conclusion, the 3D hybrid radial acquisition with compressed sensing offers a promising framework for efficient LGE imaging of the LA. This method is more robust to motion than Cartesian sampling without a respiratory navigator. The approach may be useful for identifying and tracking the effect of RF ablation in patients with atrial fibrillation. Future work of including self-navigating techniques may further improve the approach.


This work was supported in part by the Ben B. and Iris M. Margolis Foundation


1. Peters DC, Wylie JV, Hauser TH, Kissinger KV, Botnar RM, Essebag V, Josephson ME, Manning WJ. Detection of pulmonary vein and left atrial scar after catheter ablation with three-dimensional navigator-gated delayed enhancement MR imaging: initial experience. Radiology. 2007;243(3):690–695. [PubMed]
2. McGann CJ, Kholmovski EG, Oakes RS, Blauer JJE, Daccarett M, Segerson N, Airey KJ, Akoum N, Fish E, Badger TJ, DiBella EVR, Parker D, MacLeod RS, Marrouche NF. New Magnetic Resonance Imaging-Based Method for Defining the Extent of Left Atrial Wall Injury After the Ablation of Atrial Fibrillation. J Am Coll Cardiol. 2008;52(15):1263–1271. [PubMed]
3. Ishihara Y, Nazafat R, Wylie JV, Linguraru MG, Josephson ME, Howe RD, Manning WJ, Peters DC. MRI Evaluation of RF Ablation Scarring for Atrial Fibrillation Treatment. Proceedings of SPIE Medical Imaging: Visualization and Image-Guided Procedures; San Diego: 2007. p. 65090Q.
4. Oakes RS, Badger TJ, Kholmovski EG, Akoum N, Burgon NS, Fish EN, Blauer JJ, Rao SN, DiBella EV, Segerson NM, Daccarett M, Windfelder J, McGann CJ, Parker D, MacLeod RS, Marrouche NF. Detection and quantification of left atrial structural remodeling with delayed-enhancement magnetic resonance imaging in patients with atrial fibrillation. Circulation. 2009;119(13):1758–1767. [PMC free article] [PubMed]
5. Akoum N, Daccarett M, McGann C, Segerson N, Vergara G, Kuppahally S, Badger T, Burgon N, Haslam T, Kholmovski E, Macleod ROB, Marrouche N. Atrial Fibrosis Helps Select the Appropriate Patient and Strategy in Catheter Ablation of Atrial Fibrillation: A DE-MRI Guided Approach. Journal of Cardiovascular Electrophysiology. 2010 no-no. [PMC free article] [PubMed]
6. Glover GH, Pauly JM. Projection reconstruction techniques for reduction of motion effects in MRI. Magn Reson Med. 1992;28:275–289. [PubMed]
7. Block KT, Uecker M, Frahm J. Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. Magn Reson Med. 2007;57:1086–1098. [PubMed]
8. Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med. 2007;58:1182–1195. [PubMed]
9. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42(5):952–962. [PubMed]
10. Griswold MA, Jakob PM, Heidemann RM, et al. Generalized auto-calibrating partially parallel acquisitions (GRAPPA) Magn Reson Med. 2002;47:1202–1210. [PubMed]
11. Messroghli DR, Niendorf T, Schulz-Menger J, Dietz R, Friedrich MG. T1 mapping in patients with acute myocardial infarction. J Cardiovasc Magn Reson. 2003;5(2):353–359. [PubMed]
12. Messroghli DR, Greiser A, Fröhlich M, Dietz R, Schulz-Menger J. Optimization and validation of a fully-integrated pulse sequence for modified look-locker inversion-recovery (MOLLI) T1 mapping of the heart. Journal of Magnetic Resonance Imaging. 2007;26(4):1081–1086. [PubMed]
13. Larsson HB, Rosenbaum S, Fritz-Hansen T. Quantification of the effect of water exchange in dynamic contrast MRI perfusion measurements in the brain and heart. Magn Reson Med. 2001;46(2):272–281. [PubMed]
14. Adluru G, Vijayakumar S, Burgon N, Kholmovski EG, Marrouche NF, DiBella EVR. 3D Hybrid radial acquisition with compressed sensing for LGE imaging of left atrium: A simulation study. Proceedings of the 18th Annual Meeting of ISMRM; Stockholm: 2010. (abstract 1285)
15. Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. An Optimal Radial Profile Order Based on the Golden Ratio for Time-Resolved MRI. IEEE Transactions on Medical Imaging. 2007;26(1):68–76. [PubMed]
16. Knoll F, Unger M, Diwoky C, Clason C, Pock T, Stollberger R. Fast reduction of undersampling artifacts in radial MR angiography with 3D total variation on graphics hardware. Magnetic Resonance Materials in Physics, Biology and Medicine. 2010;23(2):103–114. [PubMed]
17. Candes E, Romberg J, Tao T. Robust uncertainity principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Information Theory. 2006;52:489–509.
18. Donoho DL. Compressed sensing. IEEE Trans Info Theory. 2006;52:1289–1306.
19. Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D. 1992;60:259–268.
20. Adluru G, McGann C, Speier P, Kholmovski EG, Shaaban A, Dibella EV. Acquisition and reconstruction of undersampled radial data for myocardial perfusion magnetic resonance imaging. J Magn Reson Imaging. 2009;29(2):466–473. [PMC free article] [PubMed]
21. Hansen PC. In: Computational inverse problems in electrocardiology. Johnston P, editor. Southampton: WIT Press; 2001.
22. Adluru G, Awate SP, Tasdizen T, Whitaker RT, DiBella EVR. Temporally constrained reconstruction of dynamic cardiac perfusion MRI. Magn Reson Med. 2007;57(6):1027–1036. [PubMed]
23. McLeish K, Kozerke S, Crum WR, Hill DLG. Free-breathing radial acquisitions of the heart. Magnetic Resonance in Medicine. 2004;52(5):1127–1135. [PubMed]
24. Atkinson D, Hill DLG. Reconstruction after rotational motion. Magnetic Resonance in Medicine. 2003;49(1):183–187. [PubMed]
25. Fessler JA, Sutton BP. Nonuniform fast Fourier transforms using min-max interpolation. IEEE Transactions on Signal Processing. 2003;51(2):560–574.
26. Fessler JA. [Accessed Aug 15 2010];Image reconstruction toolbox. 2010
27. Crete F, Dolmiere T, Ladret P, Nicolas M. In: The blur effect: perception and estimation with a new no-reference perceptual blur metric. Bernice ER, Thrasyvoulos NP, Scott JD, editors. SPIE; 2007. p. 64920I.
28. Peters DC, Botnar RM, Kissinger KV, Yeon SB, Appelbaum EA, Manning WJ. Inversion recovery radial MRI with interleaved projection sets. Magnetic Resonance in Medicine. 2006;55(5):1150–1156. [PubMed]
29. Marziliano P, Dufaux F, Winkler S, Ebrahimi T. A no-reference perceptual blur metric. Proceedings of the IEEE International Conference on Image Processing; Rochester: 2002. pp. III-57–III-60.
30. Xiang Z, Milanfar P. Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content. IEEE Transactions on Image Processing. 2010;19(12):3116–3132. [PubMed]