Computers in Biology and Medicine 43 (2013) 2256–2262
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Structured errors in reconstruction methods for Non-Cartesian MR data Fabio Gibiino a,b,n, Vincenzo Positano c,d, Florian Wiesinger b, Giulio Giovannetti d, Luigi Landini a,c,d, Maria Filomena Santarelli d,c a
Department of Information Engineering—EIT, University of Pisa, Pisa, Italy GE Global Research, Munich, Germany c Fondazione Gabriele Monasterio CNR–Regione Toscana,Pisa, Italy d Institute of Clinical Physiology, National Research Council, Pisa, Italy b
art ic l e i nf o
a b s t r a c t
Article history: Received 19 July 2013 Accepted 8 October 2013
Background: Reconstruction methods for Non-Cartesian magnetic resonance imaging have often been analyzed using the root mean square error (RMSE). However, RMSE is not able to measure the level of structured error associated with the reconstruction process. Methods: An index for geometric information loss was presented using the 2D autocorrelation function. The performances of Least Squares Non Uniform Fast Fourier Transform (LS-NUFFT) and gridding reconstruction (GR) methods were compared. The Direct Summation method (DS) was used as reference. For both methods, RMSE and the loss in geometric information were calculated using a digital phantom and a hyperpolarized 13C dataset. Results: The performance of the geometric information loss index was analyzed in the presence of noise. Comparing to GR, LS-NUFFT obtained a lower RMSE, but its error image appeared more structured. This was observed in both phantom and in vivo experiments. Discussion: The evaluation of geometric information loss together with the reconstruction error was important for an appropriate performance analysis of the reconstruction methods. The use of geometric information loss was helpful to determine that LS-NUFFT loses relevant information in the reconstruction process, despite the low RMSE. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Magnetic resonance imaging Autocorrelation Image structures Non-Cartesian MRI reconstruction Hyperpolarized 13C
1. Introduction Magnetic resonance (MR) is an imaging technique that acquires signals to fill the Fourier domain (K-space) of the desired image. In MR it is often necessary to acquire the signal very quickly before it decays after the excitation. Several methods have been proposed to perform a rapid acquisition, including echo planar imaging [1], compressed sensing [2], and parallel imaging with restriction to multi-channel receiver coils [3,4]. A rapid acquisition can also be performed using Non-Cartesian trajectories to fill the K-space, e.g. spiral trajectories [5,6]. Non-Cartesian trajectories starting the acquisition from the center of the K-space acquire most of the relevant information of the image before the signal decays. For Non-Cartesian acquired data, specialized reconstruction methods need to be used because the Fast Fourier Transform (FFT) algorithm cannot be directly applied. The Direct Summation
n Corresponding author at: Department of Information Engineering—University of Pisa, Via G. Caruso 16, 56122 Pisa, Italy. Tel.: þ39 3473048006. E-mail address: [email protected] (F. Gibiino).
0010-4825/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2013.10.013
(DS) reconstruction method [7,8] uses the definition of the Inverse Discrete Fourier Transform (IDFT) and is often taken as reference. The “Gridding Reconstruction” (GR) [9–11] is a fast reconstruction method using a Kaiser–Bessel window to interpolate frequency data onto Cartesian grid before applying the FFT. The Least Squares NonUniform Fast Fourier Transform (LS-NUFFT) reconstruction method [12–15] was developed from NUFFT [16,17] and is similar to gridding. LS-NUFFT obtains good performances thanks to a pseudo-inverse matrix that guides the choice for gridding interpolator. Usually the performance evaluation of reconstruction methods is limited to the calculation of the root mean square error (RMSE) [14]. The RMSE provides an estimate of the error amplitude but cannot recognize any geometrical structure in the error image. The presence of structures in the error image implies that the reconstruction method lost some of the interesting geometric information. Tissue shape preservation and correct depiction of image details are key issues in diagnostic image applications. The presence of structures in the error image is not dependent on the RMSE and can be detected using second-order statistical methods [18,19]. These methods analyze the relationships between pixels or groups of pixels (usually two), finding correlations in case of
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
structured images. The geometric information loss can be analyzed using the co-occurrence matrix [20], which measures the joint probability of pair of pixels. However the co-occurrence matrix limits the analysis to the relationship of pixels with a fixed geometric position relative to each other. Moreover this method is typically very sensitive to the noise level in the image. Alternatively, the 2D autocorrelation can provide information about the structures in the image as well [21]. The autocorrelation analyzes the relationships between pixels in any direction and at any distance on the image. Therefore, the autocorrelation provides an indication of how the structures are distributed globally in the image. The reconstruction methods were tested and compared using both a digital phantom and an experimental data set. As exemplary experimental dataset we chose 3D hyperpolarized 13C in vivo data [22] sampled with a spiral trajectory. The Chemical Shift Imaging (CSI) of infused hyperpolarized 13C-labeled compounds [23,24] allows real-time monitoring of the metabolic conversion of pyruvate to lactate, alanine and bicarbonate in both normal and malignant tissues [25–28]. The [1-13C]pyruvate tracer fast decay (about 60 s) requires a very fast acquisition strategy. This is why Non-Cartesian trajectories are often used to accelerate the acquisition. In this study, the performance of LS-NUFFT and GR methods were tested taking DS as reference for the comparison. The geometric information loss was measured analyzing the structures present in the error image, besides calculating the RMSE.
2. Materials and methods 2.1. Reconstruction algorithms 2.1.1. Direct summation (DS) ! Ideally, the general relationship between the space image Ið r Þ ! and the signal Sð k Þ acquired by the MR scanner is determined by the inverse Fourier Transform: Z - ! ! ! Ið r Þ ¼ Sð k Þei2π k U r d k ð1Þ ! ! The vector r represents the space location and k represents the frequency location in the K-space. In practice the acquisition has to be limited to a finite number (Ns) of samples. These samples are ! ! placed at frequency locations k ¼ k n , n ¼0…Ns 1. The inverse Fourier Transform in Eq. (1) needs to be discretized as follows: Ns 1 ! ! ! Ið r Þ ¼ ∑ Sð k n Þei2π k n U r Δ k n ;
n¼0
ð2Þ
2257
! where Δ k n represents a differential area and takes into account ! different densities of samples in the K-space. The value of Δ k n can be calculated using the Jacobian of the trajectory mathematical ! expression [7]. Alternatively Δ k n can be determined using the Voronoi diagram [29,30] as was done in this study. We used the area of the region associated with each point of the trajectory on the K-space. The spiral trajectory with one interleave is shown in Fig. 1a and the associated Voronoi diagram is shown in Fig. 1b. 2.1.2. Gridding reconstruction (GR) The use of the FFT algorithm makes the gridding reconstruction method fast. Before using FFT, sampled data need to be interpolated on a Cartesian grid. The interpolation process is called ‘gridding’. Considering the 1D case, S is the sampled data (i.e., the signal) and DCF is the density compensation function depending on the chosen trajectory [9]. DCF is proportional to the inverse of the ! Δ k n mentioned in the previous section. C is a chosen convolution kernel and III is the Cartesian grid where data are to be interpolated. The gridding process can thus be summarized as follows: S=DCF C III 1 C; ð3Þ Sg ¼ where Sg represents the gridded data and is the convolution operator. The gridding is performed with three main steps (i) convolution between density compensated data and a chosen kernel C (ii) re-sampling data onto a Cartesian grid (iii) scaling of the space image. Scaling corresponds to deconvolving ( 1) the obtained dataset with C in the frequency domain. A 4-pixel-wide Kaiser-–Bessel (KB) window was used as C kernel. The interpolation was calculated on a two-fold field-ofview (FOV) to reduce aliasing [9]. DCF was calculated using Voronoi diagram [29]. 2.1.3. Least Squares Non-Uniform Fast Fourier Transform (LS-NUFFT) The LS-NUFFT method is similar to the GR method because it performs an interpolation followed by inverse FFT. The main difference is the choice of the convolution kernel for the interpolation. LS-NUFFT is based on the approximation of a complex exponential, which placed into the 1D version of the IDFT formula (Eq. (2)), gives the following expression [13–15]: " #
I j sj 1
½ðωh mMΔx=2πÞ þ ðq=2Þ
Ns 1
k ¼ ½ðωh mMΔx=2πÞ ðq=2Þ
h¼0
∑
∑ ρkh αkh eik
j2π mM
ð4Þ
where ωh (h ¼1…Ns) is the pulsation of sampled data; sj are scaling factors to be used in the space domain; jΔx (j¼ 1…M) is the space
Fig. 1. (a) Example of spiral trajectory with one interleaf. Center of the K-space is sampled with higher density than periphery and (b) zoomed Voronoi diagram of a spiral trajectory with only one interleave.
2258
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
discrete coordinate and m corresponds to the factor of increase for the FOV to reduce aliasing. The interpolation kernel is referred as ρhk; it depends on the frequency location k and the pulsation ωh and its width is dependent on q. Ij is image value at pixel j and αkh gathers the sampled data at the pulsation ωh and the DCF term. The inner summation in Eq. (4) is an interpolation of αkh using ρhk as kernel. The outer summation is an Inverse Discrete Fourier Transformation. The IDFT can be evaluated with FFT because the complex exponentials are defined at equidistant frequency locations k. For GR, a Fourier transform binds scaling factors and interpolation kernel (C). For LS-NUFFT the calculation of the interpolating kernel is guided by Moore–Penrose pseudo-inverse matrix as follows: ρh;ls A† vh ¼ ðAH AÞ 1 AH vh ;
ð5Þ H
, A is a Fourier matrix and A is the hermitian of where νk ¼ sj e matrix A. ρh,ls is the interpolator that minimizes the worst-case approximation error as shown by Fessler and Sutton [13]. In this work a 2D version of LS-NUFFT was implemented as has been done in [14]. A Kaiser–Bessel shape was used for scaling factors, as suggested by Fessler. The parameter q was chosen equal to 4 to make the interpolation kernel width similar to the one used in GR. iωh jΔx
2.2. Reconstruction error analysis The reconstruction error can be analyzed calculating the error image between the tested method and the reference. In this study, the DS method was taken as reference. The error image can be modeled as the sum of two error components: structured and random. The amplitude of the overall error can be measured with the RMSE. The estimation of the structured component required calculating the 2D autocorrelation of the error image [18,21]. The presence of structures in an image can be detected using the width of the central peak of the 2D autocorrelation of the image. A wider peak is associated with more evident structures, hence a heavier loss of geometric information in the error image. Moreover, the geometry of the error image is transferred in the anisotropy of the autocorrelation central peak. The two main widths associated with the autocorrelation central peak were calculated with a fitting described by the following exponential model applied to the normalized error image: gðx; yÞ ¼ e ððx=τ1 Þ þ ðy=τ2 ÞÞ :
ð6Þ
Two decay constants τ1 and τ2 were calculated the autocorrelation image, taking the mean decay constant τ as representative of the width central peak. The ratio Rτ ¼τ1/τ2 was taken as index of anisotropy of the autocorrelation peak. The parameter τ is sensitive to the presence of a random component in the image, as any other second-order statistic metric [19]. In order to obtain an index of structured error relative to the random component level, the method of surrogate data was used [35]. After calculating τ for the image, the pixel location in the image was randomized and a new width of the autocorrelation peak (τrand) was calculated. The ratio G ¼ τ/τrand was taken as index of loss in geometric information. 2.3. Datasets Experiments were performed on both digital phantom data and a real experimental dataset. MR synthetic data were generated from a 64 64 Shepp-Logan digital phantom considering FOV ¼ 30 cm. The spectrum of the image was interpolated on spiral K-space locations [13]. To mimic spiral sampling used in real applications, a maximum distance of 1/FOV was imposed between adjacent turns of the
spiral to avoid aliasing, following the Nyquist criterion. A variable density spiral (VDS) was used as filling trajectory for the K-space [31]. The distance between adjacent turns was designed to linearly increase from 0.5/FOV to 1/FOV moving further from the center of the K-space (Fig. 2a). This allows acquiring the center of the Kspace more densely and thus gaining better signal-to-noise ratio (SNR). The point spread function (PSF) [32] was calculated to characterize the behavior of the reconstruction methods. Regarding the experimental dataset, metabolic information covering the heart was obtained using a 3D IDEAL spiral CSI with FOV of 30 cm and slab thickness of 10 cm. A 3D dataset was acquired using a pulse-and-acquire sequence with 2D spiral readout and Cartesian phase encoding (PE) along the slice-encode direction (z). The following acquisition parameters were used: 14 PE along z, 11 echoes per PE at different echo time (TE), constant TE increment¼0.9 ms, 1824 samples per echo, flip angle ¼71, receiving bandwidth ¼ 62.5 kHz. The animal experiment was performed on a healthy male mini-pig using a 3 T GE HDx Excite MR scanner (GE Healthcare, Waukesha, WI, USA) with a 13C quadrature birdcage coil. [1-13C] pyruvate hyperpolarization was obtained using a HyperSense DNP polarizer (Oxford Instruments, Oxford, UK) with subsequent dissolution. Before the reconstruction, the experimental raw data were processed with the IDEAL method as described in [33,34]. The reconstruction methods were applied separately for the 2D images of each slice. The acquired VDS trajectory was similar to the one used for simulations. A zero padding interpolation (Fig. 2b) was applied to obtain 64 64 64 images. A Gaussian filter was used to reduce the truncation artifact and to low-pass filter the data, obtaining an actual in-plane resolution of about 12 mm.
3. Results The performance of τ was analyzed simulating the progressive presence of structures and random error (Fig. 3a). A 128 128 Shepp-Logan phantom was used for the simulation and the results were averaged over 30 repetitions. The parameter τ is sensitive to different levels of structures present in the image. However, as the level of simulated random error increases, τ loses sensibility to the image structures. Fig. 3b shows the behavior of parameter G obtained in the same simulation. As the intensity of the random error increases, G approaches 1. Conversely, when detectable structures are present in the error image, G increases. G can be considered equal to 1 for error levels higher than about 3 times the maximum signal of the structured error. The PSF was calculated for both reconstruction methods. The PSF characterizes how the method behaves with a general reconstructed image. The error image and the autocorrelation were calculated for the PSF (Fig. 4). The autocorrelation of the PSF decreased more slowly than in the case of LS-NUFFT than GR. From the PSF, the error associated with LS-NUFFT appears more correlated than the error associated with GR. For the digital phantom experiments on the reconstruction methods, GR obtained τ¼12.9 mm, G¼ 35.0, Rτ ¼1.4 and RMSE¼ 5.1 10 4, whereas LS-NUFFT obtained τ¼52.6 mm, G¼ 9.3, Rτ ¼1.6 and RMSE¼9.6 10 5. Therefore, LS-NUFFT allows the reconstruction to have an improvement of 5.3 times in the RMSE, as previously shown [14]. However, in the case of LS-NUFFT the structured error is 35 times more relevant than the random error, whereas only 9.3 times in the case of GR. Fig. 5 shows the error images and autocorrelations for both methods. The shape of the autocorrelation follows the one of the reconstructed object (Fig. 5d and e). The parameter Rτ keeps track of the image anisotropy. In the case of an object with one dimension longer than the other one, the correlation
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
2259
Fig. 2. (a) VDS used for simulation and experiments. The distance between adjacent turns linearly increases moving toward periphery of K-space and (b) Modulus of K-space acquired data plotted onto the sampling trajectory. The width of the trajectory does not cover the inverse of spatial resolution. This implies zero padding and thus blurring on the reconstructed image.
Fig. 3. Analysis of the structure detection performance in presence of noise. The noise level is the ratio between the amplitude of the noise variance and the structure maximum signal. The bottom row depicts the decreasing levels of structure (S1, S2, S3, S4) used in the analysis. The behavior of τ and τrand are shown in (a) as the noise level increases and for different level structures. The behavior of the relative index G is shown in (b).
will be preserved for longer distance than for the shorter dimension and Rτ will be bigger than one. The values of τ, G, Rτ and RMSE were evaluated for experimental data as well (Table 1). Fig. 6 shows the exemplary result images for pyruvate on central slice.
4. Discussion
Fig. 4. Middle line of the normalized 2D autocorrelation for the PSF for both reconstruction methods. The curve for LS-NUFFT decreases more slowly than GR. Therefore, the PSF analysis shows that the error image of the LS-NUFFT is more structured than GR.
In this study two reconstruction algorithms (LS-NUFFT and GR) were examined and quantitatively compared in terms of reconstruction accuracy, taking the DS algorithm as reference. LS-NUFFT showed a lower RMSE compared to GR in both simulation and in vivo experiments, as confirmed by previous studies [14]. In this work the analysis was extended to the geometrical information loss. An index of geometric information loss (G) was presented and calculated using the 2D autocorrelation of the error images.
2260
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
Fig. 5. Simulation results. (a) Reconstructed image with DS method. The image has been normalized by its maximum value. (b and c) Error image of GR method (b) and LS-NUFFT (c) using DS as reference. FOV is circular for both images. (d and e) Image of the 2D autocorrelation for the error image in GR (d) and LS-NUFFT (e) method. Table 1 Results obtained with the experimental dataset for each metabolite (lactate, hydrate pyruvate, alanine, bicarbonate and pyruvate). The values showed are averaged over the 64 slices of the reconstructed dataset. lactate
hydrate pyruvate
alanine
bicarbonate
pyruvate
GR RMSEx10 5 τ [mm] Rτ G
3.4 7 0.4 9.3 7 0.6 1.147 0.05 6.5 7 0.5
3.4 70.4 12.4 71.7 1.13 70.05 7.9 70.8
3.9 7 0.6 16.6 7 2.5 1.20 7 0.07 8.4 7 0.8
3.6 7 0.6 10.9 7 1.1 1.08 7 0.05 7.0 7 0.6
12.0 75.5 13.3 71.6 1.16 70.06 11.8 71.8
LS-NUFFT RMSEx10 5 τ [mm] Rτ G
0.89 7 0.12 53.3 7 3.1 1.477 0.08 13.6 7 0.5
0.89 70.14 55.1 71.9 1.50 70.03 13.0 70.8
1.317 0.34 57.8 7 1.9 1.577 0.02 10.8 7 1.8
0.95 7 0.20 55.3 7 2.9 1.54 7 0.04 13.8 7 0.9
2.76 71.17 44.1 78.5 1.44 70.12 17.3 70.9
The use of the 2D autocorrelation to quantify the presence of structures in the image was shown to be effective. Fig. 3 shows the limitation associated with the use of the width of the autocorrelation peak directly. The parameter τ is an absolute measure that does not allow distinguishing structured errors from random errors. The index G was introduced in order to have relative index for geometrical loss, valid for any level of structure in the image. The index G measures how much the structured errors are relevant compared to the random errors. Fig. 3b suggests that G has a reduced dependency on structures in the image, comparing to τ. In fact, especially in the case of several structures, the surrogate method does not cause the total loss of correlation in the image even randomizing the pixel locations. A memory of the behavior of τ can still be found in τrand, thus making G less dependent on the structure level. The estimation of geometric information loss as performance index represents an important improvement in evaluation of reconstruction methods. Previous studies only focused on the RMSE of reconstruction which is only able to give information about the first-order statistics of the reconstruction error. The presence of structures in the error image can be detected only by using a second-order statistic index, like the autocorrelation. Structures may be present in the error image also in the case of
low RMSE. Therefore a performance analysis limited to the RMSE does not benefit from all the information available in the error image. Despite the lower RMSE, LS-NUFFT showed an error image more structured than GR. A comparison limited to the RMSE would not have been able to recognize that also LS-NUFFT loses relevant information in the reconstruction process. For both the simulated and experimental datasets, LS-NUFFT obtained a G value higher than GR. However, the relative difference between the G values of the two reconstruction methods is higher in the case of the simulated data. This is due to the different levels of structure and acquisition noise in the analyzed images. In fact, the acquisition noise increases the random error component, hence decreasing the value of G. The level of noise present in the reconstructed image also depends on the reconstruction method as suggested in [36]. The shape and the size of the interpolation kernel for gridding methods is a key issue in determining the accuracy of the reconstruction [11–13]. For example, the characteristics of the gridding kernel can change the amplitude of the reconstruction error or the shading on the reconstructed image. For LS-NUFFT the choice of the gridding kernel is guided by the pseudo-inverse matrix, hence minimizing the approximation error. For LS-NUFFT the size of the kernel can cause relevant reconstruction error
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
2261
Fig. 6. Experimental results shown for Pyruvate on central slice. (a) Modulus of reconstructed image using DS. Image has been normalized by its maximum value. (b and c) error image of GR method (b) and LS-NUFFT (c) using DS as reference. FOV is circular for both images. (d and e) Image of the 2D autocorrelation for the error image in GR (d) and LS-NUFFT (e) method.
related to the matrix pseudo-inverse calculation [15]. Therefore the kernel size was carefully chosen in order to make the two reconstruction methods comparable and avoid the need for a regularization method for the pseudo-inverse calculation. In this study the same weighting function was used for the data to compare the reconstructions methods. The DCF was evaluated using the Voronoi diagram for all three methods considered. DCF has a relevant influence on the reconstructed image, as proved by the extensive studies on its shape [36–39]. LS-NUFFT differs from GR only in operations coming after data weighting. The use of different DCFs would not allow distinguishing the differences in reconstructed images due to DCF, and the ones due to the interpolator. In this work the effect of the DCF on the reconstructed image was not considered. The study of the effect of DCF over the loss in geometric information could be an interesting development of the present work. The use of the 2D autocorrelation function could also be useful in parallel imaging applications [3,4]. Parallel imaging techniques commonly suffer from coherent aliasing artifacts, which could easily be detected by the 2D autocorrelation function. The presence of an image ghost would be detected with a secondary peak placed in the same spatial location of the ghost. For this application the amplitude of the autocorrelation peak, rather than the width, could be useful to quantify the intensity of the ghost.
5. Conclusions This work compared the LS-NUFFT and the GR reconstruction methods. The performance was tested on simulations with modified Shepp-Logan phantom and on hyperpolarized 13C 3D CSI in vivo data as exemplary experimental data set. Even though LSNUFFT obtained a smaller RMSE than GR, the error image of LSNUFFT appeared more structured than GR. These results strongly suggest that both reconstruction error and geometrical loss of information indices should be adopted for a complete evaluation of reconstruction methods. Even though LS-NUFFT shows better
reconstruction accuracy, it is relevant to be aware that the reconstruction process lost interesting information and not only noise.
Conflict of interest None declared. References [1] S. Hu, M. Lustig, A.P. Chen, J. Crane, A. Kerr, D.A. Kelley, R. Hurd, J. Kurhanewicz, S.J. Nelson, J.M. Pauly, D.B. Vigneron, Compressed sensing for resolution enhancement of hyperpolarized 13C flyback 3D-MRSI, J. Magn. Reson. 192 (2) (2008) 258–264. [2] M. Lustig, D. Donoho, J.M. Pauly., MRI Sparse, The application of compressed sensing for rapid MR imaging, Magn. Reson. Med. 58 (6) (2007) 1182–1195. [3] K.P. Pruessmann, M. Weiger, M.B. Scheidegger, P. Boesiger., SENSE: sensitivity encoding for fast MRI, Magn. Reson. Med. 42 (5) (1999) 952–962. [4] M.A. Griswold, P.M. Jakob, R.M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, A. Haase, Generalized autocalibrating partially parallel acquisitions (GRAPPA), Magn. Reson. Med. 47 (6) (2002) 1202–1210. [5] A.Z. Lau, A.P. Chen, N.R. Ghugre, V. Ramanan, W.W. Lam, K.A. Connelly, G. A. Wright, C.H. Cunningham, Rapid Multislice Imaging of Hyperpolarized 13C Pyruvate and Bicarbonate in the heart, Magn. Reson. Med. 64 (5) (2010) 1323–1331. [6] C.K. Turnes, J. Romberg. Spiral F.F.T.: An efficient method for 3-D FFTs on spiral mri contours. in: Proceedings of IEEE 17th International Conference on Image Processing, 2010. [7] R.D. Hoge, R.K.S. Kwan, G.B. Pike, Density compensation functions for spiral MRI, Magn. Reson. Med. 38 (1996) 117–128. [8] G.E. Sarty, Natural k-plane coordinate reconstruction method for magnetic resonance imaging: mathematical foundations, Int. J. Imag. Syst. Technol. 8 (1997) 519–528. [9] J.I. Jackson, C.H. Meyer, D.G. Nishimura, A. Macovski, Selection of a convolution function for Fourier inversion using gridding, IEEE Trans. Med. Imag. 10 (1991) 473–478. [10] J.D. O’Sullivan, A. Fast Sinc, Function gridding algorithm for fourier inversion in computer tomography, IEEE Trans. Med. Imag. 4 (4) (1985) 200–207. [11] P.J. Beatty, D.G. Nishimura, J.M. Pauly, Rapid gridding reconstruction with a minimal oversampling ratio, IEEE Trans. Med. Imag. 24 (2005) 799–808. [12] N. Nguyen, Q.H. Liu, The regular fourier matrices and nonuniform fast fourier transforms, SIAM J. Sci. Comput. 21 (1) (1999) 283–293. [13] J.A. Fessler, B.P. Sutton, Nonuniform fast fourier transforms using min-max interpolation, IEEE Trans. Sign. Proc. 51 (2) (2003) 560–574. [14] L. Sha, H. Guo, A.W. Song, An improved gridding method for spiral MRI using nonuniform fast Fourier transform, J. Magn. Res. 162 (2003) 250–258.
2262
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
[15] F. Gibiino, V. Positano, L. Landini, M.F. Santarelli, Regularization techniques on non-uniform Fast Fourier Transform, Int. J. Num. Meth. Biomed. Eng. 29 (5) (2013) 561–573. [16] A. Dutt, V. Rokhlin, Fast fourier transforms for nonequispaced data, SIAM J. Sci. Comput. 14 (6) (1993) 1368–1393. [17] M. Singh, M. Mandal, A. Basu Gaussian and Laplacian of Gaussian weighted functions for robust feature based tracking. Pattern Recognition Lett. 2005;26:1995-2005. [18] G.N. Srinivasan, Shobha G., Statistical texture analysis. in: Proceedings of world academy of science, engineering and technology, 2008, pp. 2070–3740. [19] W.K. Pratt, Digital Image Processing: PISK Inside, 16, John Wiley & Sons, Inc, New York, USA (2001) 555–577. (Chapter). [20] R.M. Haralick, L. Watson, A facet model for image data, Comput. Vision Graphics Image Process. 15 (1981) 113–129. [21] D.M. Rubin, A simple autocorrelation algorithm for determining grain size from digital images of sediment, J. Sediment. Res. 74 (1) (2004) 160–165. [22] P. Henry, G. Adriany, D. Deelchand, R. Gruetter, M. Marjanska, G. Oz, E. R. Seaquist, A. Shestov, K. Ugurbil, In vivo 13C NMR spectroscopy and metabolic modelling in the brain: a practical perspective, Magn. Res. Imaging 24 (2006) 527–539. [23] K.J. Neurohr, E.J. Barrett, R.G. Shulman, In vivo carbon-13 nuclear magnetic resonance studies of heart metabolism, Proc. Natl. Acad. Sci. USA 80 (1983) 1603–1607. [24] J.H. Ardenkjaer-Larsen, B. Fridlund, A. Gram, G. Hansson, L. Hansson, M. H. Lerche, R. Servin, M. Thaning, K. Golman, Increase in signal-to-noise ratio of 410,000 times in liquid-state NMR, Proc. Natl. Acad. Sci. USA 100 (18) (2003) 10158–10163. [25] K. Golman, R. in’t Zandt, M. Thaning, Real-time metabolic imaging, Proc. Natl. Acad. Sci. USA 103 (2006) 11270–11275. [26] L. Menichetti, F. Frijia, A. Flori, F. Wiesinger, V. Lionetti, G. Giovannetti, G. D. Aquaro, F.A. Recchia, J.H. Ardenkjaer-Larsen, M.F. Santarelli, M. Lombardi, Assessment of real-time myocardial uptake and enzymatic conversion of hyperpolarized [1-13C]pyruvate in pigs with slice selective Magnetic Resonance Spectroscopy, Contrast Media Mol. Imaging 7 (1) (2012) 85–94. [27] G.D. Aquaro, F. Frijia, V. Positano, L. Menichetti, M.F. Santarelli, J.H. ArdenkjaerLarsen, F. Wiesinger, V. Lionetti, S.L. Romano, G. Bianchi, D. Neglia, G. Giovannetti, R. F. Schulte, F.A. Recchia, L. Landini, M. Lombardi, 3D CMR mapping of metabolism by
[28]
[29]
[30]
[31]
[32] [33]
[34]
[35] [36]
[37] [38] [39]
hyperpolarized 13C-Pyruvate in ischemia-reperfusion, JACC Cardiovasc. Imaging 6 (6) (2013) 743–744. C.H. Cunningham, A.P. Chen, M. Lustig, B.A. Hargreaves, J. Lupo, D. Xu, J. Kurhanewicz, R.E. Hurd, J.M. Pauly, S.J. Nelson, BD.B. Vigneron, Pulse sequence for dynamic volumetric imaging of hyperpolarized metabolic products, J. Magn. Res. 193 (2008) 139–146. V. Rasche, R. Proksa, R. Sinkus, P. Börnert, H. Eggers, Resampling of data between arbitrary grids using convolution interpolation, IEEE Trans. Med. Imaging 18 (5) (1999) 385–392. D. Moratal, A.V. Lluch, V. Bodí, M.E. Brummer, Magnetic resonance imaging gridding reconstruction methods with and without density compensation functions, IEEE Lat. Am. Trans. 9 (1) (2011) 774–778. M.S. Sussman, J.A. Stainsby, N. Robert, N. Merchant, G.A. Wright, Variabledensity adaptive imaging for high-resolution coronary artery MRI, Magn. Reson. Med. 48 (2002) 753–764. J.G. Pipe, Reconstructing MR images from undersampled data: data-weighting considerations, Magn. Reson. Med. 43 (6) (2000) 867–875. E.K. Brodsky, J.H. Holmes, H. Yu, S.B. Reeder, Generalized k-Space decomposition with chemical shift correction for Non-Cartesian water-fat imaging, Mag. Reson. Med. 59 (2008) 1151–1164. F. Wiesinger, E. Weidl, M.I. Menzel, M.A. Janich, O. Khegai, S.J. Glaser, A. Haase, M. Schwaiger, R.F. Schulte, IDEAL Spiral CSI for dynamic metabolic mr imaging of hyperpolarized [1-13C], Pyruvate. Magn. Reson. Med. 68 (2011) 8–16. J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, Testing for nonlinearity in time series: the method of surrogate data, Physica D 58 (1992) 77–94. F. Gibiino, V. Positano, G. Giovannetti, F. Frijia, L. Menichetti, J.H. ArdenkjaerLarsen, F. Wiesinger, R. Schulte, D. De Marchi, V. Lionetti, G. Aquaro, M. Lombardi, L. Landini, M.F. Santarelli, Reconstruction methods from hyperpolarized 13C chemical shift imaging spiral 3D data: comparison between direct summation and gridding method, ISBI proceeding (2012) 614–617. J.G. Pipe, P. Menon, Sampling density compensation in mri: rationale and an iterative numerical solution, Magn. Reson. Med. 41 (1999) 179–186. D. Rosenfeld, An optimal and efficient new gridding algorithm using singular value decomposition, Magn. Reson. Med. 40 (1998) 14–23. H. Sedarat, D.G. Nishimura, On the optimality of the gridding reconstruction algorithm, IEEE Trans. Med. Imaging 19 (4) (2000) 306–317.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Structured errors in reconstruction methods for Non-Cartesian MR data Fabio Gibiino a,b,n, Vincenzo Positano c,d, Florian Wiesinger b, Giulio Giovannetti d, Luigi Landini a,c,d, Maria Filomena Santarelli d,c a
Department of Information Engineering—EIT, University of Pisa, Pisa, Italy GE Global Research, Munich, Germany c Fondazione Gabriele Monasterio CNR–Regione Toscana,Pisa, Italy d Institute of Clinical Physiology, National Research Council, Pisa, Italy b
art ic l e i nf o
a b s t r a c t
Article history: Received 19 July 2013 Accepted 8 October 2013
Background: Reconstruction methods for Non-Cartesian magnetic resonance imaging have often been analyzed using the root mean square error (RMSE). However, RMSE is not able to measure the level of structured error associated with the reconstruction process. Methods: An index for geometric information loss was presented using the 2D autocorrelation function. The performances of Least Squares Non Uniform Fast Fourier Transform (LS-NUFFT) and gridding reconstruction (GR) methods were compared. The Direct Summation method (DS) was used as reference. For both methods, RMSE and the loss in geometric information were calculated using a digital phantom and a hyperpolarized 13C dataset. Results: The performance of the geometric information loss index was analyzed in the presence of noise. Comparing to GR, LS-NUFFT obtained a lower RMSE, but its error image appeared more structured. This was observed in both phantom and in vivo experiments. Discussion: The evaluation of geometric information loss together with the reconstruction error was important for an appropriate performance analysis of the reconstruction methods. The use of geometric information loss was helpful to determine that LS-NUFFT loses relevant information in the reconstruction process, despite the low RMSE. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Magnetic resonance imaging Autocorrelation Image structures Non-Cartesian MRI reconstruction Hyperpolarized 13C
1. Introduction Magnetic resonance (MR) is an imaging technique that acquires signals to fill the Fourier domain (K-space) of the desired image. In MR it is often necessary to acquire the signal very quickly before it decays after the excitation. Several methods have been proposed to perform a rapid acquisition, including echo planar imaging [1], compressed sensing [2], and parallel imaging with restriction to multi-channel receiver coils [3,4]. A rapid acquisition can also be performed using Non-Cartesian trajectories to fill the K-space, e.g. spiral trajectories [5,6]. Non-Cartesian trajectories starting the acquisition from the center of the K-space acquire most of the relevant information of the image before the signal decays. For Non-Cartesian acquired data, specialized reconstruction methods need to be used because the Fast Fourier Transform (FFT) algorithm cannot be directly applied. The Direct Summation
n Corresponding author at: Department of Information Engineering—University of Pisa, Via G. Caruso 16, 56122 Pisa, Italy. Tel.: þ39 3473048006. E-mail address: [email protected] (F. Gibiino).
0010-4825/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2013.10.013
(DS) reconstruction method [7,8] uses the definition of the Inverse Discrete Fourier Transform (IDFT) and is often taken as reference. The “Gridding Reconstruction” (GR) [9–11] is a fast reconstruction method using a Kaiser–Bessel window to interpolate frequency data onto Cartesian grid before applying the FFT. The Least Squares NonUniform Fast Fourier Transform (LS-NUFFT) reconstruction method [12–15] was developed from NUFFT [16,17] and is similar to gridding. LS-NUFFT obtains good performances thanks to a pseudo-inverse matrix that guides the choice for gridding interpolator. Usually the performance evaluation of reconstruction methods is limited to the calculation of the root mean square error (RMSE) [14]. The RMSE provides an estimate of the error amplitude but cannot recognize any geometrical structure in the error image. The presence of structures in the error image implies that the reconstruction method lost some of the interesting geometric information. Tissue shape preservation and correct depiction of image details are key issues in diagnostic image applications. The presence of structures in the error image is not dependent on the RMSE and can be detected using second-order statistical methods [18,19]. These methods analyze the relationships between pixels or groups of pixels (usually two), finding correlations in case of
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
structured images. The geometric information loss can be analyzed using the co-occurrence matrix [20], which measures the joint probability of pair of pixels. However the co-occurrence matrix limits the analysis to the relationship of pixels with a fixed geometric position relative to each other. Moreover this method is typically very sensitive to the noise level in the image. Alternatively, the 2D autocorrelation can provide information about the structures in the image as well [21]. The autocorrelation analyzes the relationships between pixels in any direction and at any distance on the image. Therefore, the autocorrelation provides an indication of how the structures are distributed globally in the image. The reconstruction methods were tested and compared using both a digital phantom and an experimental data set. As exemplary experimental dataset we chose 3D hyperpolarized 13C in vivo data [22] sampled with a spiral trajectory. The Chemical Shift Imaging (CSI) of infused hyperpolarized 13C-labeled compounds [23,24] allows real-time monitoring of the metabolic conversion of pyruvate to lactate, alanine and bicarbonate in both normal and malignant tissues [25–28]. The [1-13C]pyruvate tracer fast decay (about 60 s) requires a very fast acquisition strategy. This is why Non-Cartesian trajectories are often used to accelerate the acquisition. In this study, the performance of LS-NUFFT and GR methods were tested taking DS as reference for the comparison. The geometric information loss was measured analyzing the structures present in the error image, besides calculating the RMSE.
2. Materials and methods 2.1. Reconstruction algorithms 2.1.1. Direct summation (DS) ! Ideally, the general relationship between the space image Ið r Þ ! and the signal Sð k Þ acquired by the MR scanner is determined by the inverse Fourier Transform: Z - ! ! ! Ið r Þ ¼ Sð k Þei2π k U r d k ð1Þ ! ! The vector r represents the space location and k represents the frequency location in the K-space. In practice the acquisition has to be limited to a finite number (Ns) of samples. These samples are ! ! placed at frequency locations k ¼ k n , n ¼0…Ns 1. The inverse Fourier Transform in Eq. (1) needs to be discretized as follows: Ns 1 ! ! ! Ið r Þ ¼ ∑ Sð k n Þei2π k n U r Δ k n ;
n¼0
ð2Þ
2257
! where Δ k n represents a differential area and takes into account ! different densities of samples in the K-space. The value of Δ k n can be calculated using the Jacobian of the trajectory mathematical ! expression [7]. Alternatively Δ k n can be determined using the Voronoi diagram [29,30] as was done in this study. We used the area of the region associated with each point of the trajectory on the K-space. The spiral trajectory with one interleave is shown in Fig. 1a and the associated Voronoi diagram is shown in Fig. 1b. 2.1.2. Gridding reconstruction (GR) The use of the FFT algorithm makes the gridding reconstruction method fast. Before using FFT, sampled data need to be interpolated on a Cartesian grid. The interpolation process is called ‘gridding’. Considering the 1D case, S is the sampled data (i.e., the signal) and DCF is the density compensation function depending on the chosen trajectory [9]. DCF is proportional to the inverse of the ! Δ k n mentioned in the previous section. C is a chosen convolution kernel and III is the Cartesian grid where data are to be interpolated. The gridding process can thus be summarized as follows: S=DCF C III 1 C; ð3Þ Sg ¼ where Sg represents the gridded data and is the convolution operator. The gridding is performed with three main steps (i) convolution between density compensated data and a chosen kernel C (ii) re-sampling data onto a Cartesian grid (iii) scaling of the space image. Scaling corresponds to deconvolving ( 1) the obtained dataset with C in the frequency domain. A 4-pixel-wide Kaiser-–Bessel (KB) window was used as C kernel. The interpolation was calculated on a two-fold field-ofview (FOV) to reduce aliasing [9]. DCF was calculated using Voronoi diagram [29]. 2.1.3. Least Squares Non-Uniform Fast Fourier Transform (LS-NUFFT) The LS-NUFFT method is similar to the GR method because it performs an interpolation followed by inverse FFT. The main difference is the choice of the convolution kernel for the interpolation. LS-NUFFT is based on the approximation of a complex exponential, which placed into the 1D version of the IDFT formula (Eq. (2)), gives the following expression [13–15]: " #
I j sj 1
½ðωh mMΔx=2πÞ þ ðq=2Þ
Ns 1
k ¼ ½ðωh mMΔx=2πÞ ðq=2Þ
h¼0
∑
∑ ρkh αkh eik
j2π mM
ð4Þ
where ωh (h ¼1…Ns) is the pulsation of sampled data; sj are scaling factors to be used in the space domain; jΔx (j¼ 1…M) is the space
Fig. 1. (a) Example of spiral trajectory with one interleaf. Center of the K-space is sampled with higher density than periphery and (b) zoomed Voronoi diagram of a spiral trajectory with only one interleave.
2258
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
discrete coordinate and m corresponds to the factor of increase for the FOV to reduce aliasing. The interpolation kernel is referred as ρhk; it depends on the frequency location k and the pulsation ωh and its width is dependent on q. Ij is image value at pixel j and αkh gathers the sampled data at the pulsation ωh and the DCF term. The inner summation in Eq. (4) is an interpolation of αkh using ρhk as kernel. The outer summation is an Inverse Discrete Fourier Transformation. The IDFT can be evaluated with FFT because the complex exponentials are defined at equidistant frequency locations k. For GR, a Fourier transform binds scaling factors and interpolation kernel (C). For LS-NUFFT the calculation of the interpolating kernel is guided by Moore–Penrose pseudo-inverse matrix as follows: ρh;ls A† vh ¼ ðAH AÞ 1 AH vh ;
ð5Þ H
, A is a Fourier matrix and A is the hermitian of where νk ¼ sj e matrix A. ρh,ls is the interpolator that minimizes the worst-case approximation error as shown by Fessler and Sutton [13]. In this work a 2D version of LS-NUFFT was implemented as has been done in [14]. A Kaiser–Bessel shape was used for scaling factors, as suggested by Fessler. The parameter q was chosen equal to 4 to make the interpolation kernel width similar to the one used in GR. iωh jΔx
2.2. Reconstruction error analysis The reconstruction error can be analyzed calculating the error image between the tested method and the reference. In this study, the DS method was taken as reference. The error image can be modeled as the sum of two error components: structured and random. The amplitude of the overall error can be measured with the RMSE. The estimation of the structured component required calculating the 2D autocorrelation of the error image [18,21]. The presence of structures in an image can be detected using the width of the central peak of the 2D autocorrelation of the image. A wider peak is associated with more evident structures, hence a heavier loss of geometric information in the error image. Moreover, the geometry of the error image is transferred in the anisotropy of the autocorrelation central peak. The two main widths associated with the autocorrelation central peak were calculated with a fitting described by the following exponential model applied to the normalized error image: gðx; yÞ ¼ e ððx=τ1 Þ þ ðy=τ2 ÞÞ :
ð6Þ
Two decay constants τ1 and τ2 were calculated the autocorrelation image, taking the mean decay constant τ as representative of the width central peak. The ratio Rτ ¼τ1/τ2 was taken as index of anisotropy of the autocorrelation peak. The parameter τ is sensitive to the presence of a random component in the image, as any other second-order statistic metric [19]. In order to obtain an index of structured error relative to the random component level, the method of surrogate data was used [35]. After calculating τ for the image, the pixel location in the image was randomized and a new width of the autocorrelation peak (τrand) was calculated. The ratio G ¼ τ/τrand was taken as index of loss in geometric information. 2.3. Datasets Experiments were performed on both digital phantom data and a real experimental dataset. MR synthetic data were generated from a 64 64 Shepp-Logan digital phantom considering FOV ¼ 30 cm. The spectrum of the image was interpolated on spiral K-space locations [13]. To mimic spiral sampling used in real applications, a maximum distance of 1/FOV was imposed between adjacent turns of the
spiral to avoid aliasing, following the Nyquist criterion. A variable density spiral (VDS) was used as filling trajectory for the K-space [31]. The distance between adjacent turns was designed to linearly increase from 0.5/FOV to 1/FOV moving further from the center of the K-space (Fig. 2a). This allows acquiring the center of the Kspace more densely and thus gaining better signal-to-noise ratio (SNR). The point spread function (PSF) [32] was calculated to characterize the behavior of the reconstruction methods. Regarding the experimental dataset, metabolic information covering the heart was obtained using a 3D IDEAL spiral CSI with FOV of 30 cm and slab thickness of 10 cm. A 3D dataset was acquired using a pulse-and-acquire sequence with 2D spiral readout and Cartesian phase encoding (PE) along the slice-encode direction (z). The following acquisition parameters were used: 14 PE along z, 11 echoes per PE at different echo time (TE), constant TE increment¼0.9 ms, 1824 samples per echo, flip angle ¼71, receiving bandwidth ¼ 62.5 kHz. The animal experiment was performed on a healthy male mini-pig using a 3 T GE HDx Excite MR scanner (GE Healthcare, Waukesha, WI, USA) with a 13C quadrature birdcage coil. [1-13C] pyruvate hyperpolarization was obtained using a HyperSense DNP polarizer (Oxford Instruments, Oxford, UK) with subsequent dissolution. Before the reconstruction, the experimental raw data were processed with the IDEAL method as described in [33,34]. The reconstruction methods were applied separately for the 2D images of each slice. The acquired VDS trajectory was similar to the one used for simulations. A zero padding interpolation (Fig. 2b) was applied to obtain 64 64 64 images. A Gaussian filter was used to reduce the truncation artifact and to low-pass filter the data, obtaining an actual in-plane resolution of about 12 mm.
3. Results The performance of τ was analyzed simulating the progressive presence of structures and random error (Fig. 3a). A 128 128 Shepp-Logan phantom was used for the simulation and the results were averaged over 30 repetitions. The parameter τ is sensitive to different levels of structures present in the image. However, as the level of simulated random error increases, τ loses sensibility to the image structures. Fig. 3b shows the behavior of parameter G obtained in the same simulation. As the intensity of the random error increases, G approaches 1. Conversely, when detectable structures are present in the error image, G increases. G can be considered equal to 1 for error levels higher than about 3 times the maximum signal of the structured error. The PSF was calculated for both reconstruction methods. The PSF characterizes how the method behaves with a general reconstructed image. The error image and the autocorrelation were calculated for the PSF (Fig. 4). The autocorrelation of the PSF decreased more slowly than in the case of LS-NUFFT than GR. From the PSF, the error associated with LS-NUFFT appears more correlated than the error associated with GR. For the digital phantom experiments on the reconstruction methods, GR obtained τ¼12.9 mm, G¼ 35.0, Rτ ¼1.4 and RMSE¼ 5.1 10 4, whereas LS-NUFFT obtained τ¼52.6 mm, G¼ 9.3, Rτ ¼1.6 and RMSE¼9.6 10 5. Therefore, LS-NUFFT allows the reconstruction to have an improvement of 5.3 times in the RMSE, as previously shown [14]. However, in the case of LS-NUFFT the structured error is 35 times more relevant than the random error, whereas only 9.3 times in the case of GR. Fig. 5 shows the error images and autocorrelations for both methods. The shape of the autocorrelation follows the one of the reconstructed object (Fig. 5d and e). The parameter Rτ keeps track of the image anisotropy. In the case of an object with one dimension longer than the other one, the correlation
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
2259
Fig. 2. (a) VDS used for simulation and experiments. The distance between adjacent turns linearly increases moving toward periphery of K-space and (b) Modulus of K-space acquired data plotted onto the sampling trajectory. The width of the trajectory does not cover the inverse of spatial resolution. This implies zero padding and thus blurring on the reconstructed image.
Fig. 3. Analysis of the structure detection performance in presence of noise. The noise level is the ratio between the amplitude of the noise variance and the structure maximum signal. The bottom row depicts the decreasing levels of structure (S1, S2, S3, S4) used in the analysis. The behavior of τ and τrand are shown in (a) as the noise level increases and for different level structures. The behavior of the relative index G is shown in (b).
will be preserved for longer distance than for the shorter dimension and Rτ will be bigger than one. The values of τ, G, Rτ and RMSE were evaluated for experimental data as well (Table 1). Fig. 6 shows the exemplary result images for pyruvate on central slice.
4. Discussion
Fig. 4. Middle line of the normalized 2D autocorrelation for the PSF for both reconstruction methods. The curve for LS-NUFFT decreases more slowly than GR. Therefore, the PSF analysis shows that the error image of the LS-NUFFT is more structured than GR.
In this study two reconstruction algorithms (LS-NUFFT and GR) were examined and quantitatively compared in terms of reconstruction accuracy, taking the DS algorithm as reference. LS-NUFFT showed a lower RMSE compared to GR in both simulation and in vivo experiments, as confirmed by previous studies [14]. In this work the analysis was extended to the geometrical information loss. An index of geometric information loss (G) was presented and calculated using the 2D autocorrelation of the error images.
2260
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
Fig. 5. Simulation results. (a) Reconstructed image with DS method. The image has been normalized by its maximum value. (b and c) Error image of GR method (b) and LS-NUFFT (c) using DS as reference. FOV is circular for both images. (d and e) Image of the 2D autocorrelation for the error image in GR (d) and LS-NUFFT (e) method. Table 1 Results obtained with the experimental dataset for each metabolite (lactate, hydrate pyruvate, alanine, bicarbonate and pyruvate). The values showed are averaged over the 64 slices of the reconstructed dataset. lactate
hydrate pyruvate
alanine
bicarbonate
pyruvate
GR RMSEx10 5 τ [mm] Rτ G
3.4 7 0.4 9.3 7 0.6 1.147 0.05 6.5 7 0.5
3.4 70.4 12.4 71.7 1.13 70.05 7.9 70.8
3.9 7 0.6 16.6 7 2.5 1.20 7 0.07 8.4 7 0.8
3.6 7 0.6 10.9 7 1.1 1.08 7 0.05 7.0 7 0.6
12.0 75.5 13.3 71.6 1.16 70.06 11.8 71.8
LS-NUFFT RMSEx10 5 τ [mm] Rτ G
0.89 7 0.12 53.3 7 3.1 1.477 0.08 13.6 7 0.5
0.89 70.14 55.1 71.9 1.50 70.03 13.0 70.8
1.317 0.34 57.8 7 1.9 1.577 0.02 10.8 7 1.8
0.95 7 0.20 55.3 7 2.9 1.54 7 0.04 13.8 7 0.9
2.76 71.17 44.1 78.5 1.44 70.12 17.3 70.9
The use of the 2D autocorrelation to quantify the presence of structures in the image was shown to be effective. Fig. 3 shows the limitation associated with the use of the width of the autocorrelation peak directly. The parameter τ is an absolute measure that does not allow distinguishing structured errors from random errors. The index G was introduced in order to have relative index for geometrical loss, valid for any level of structure in the image. The index G measures how much the structured errors are relevant compared to the random errors. Fig. 3b suggests that G has a reduced dependency on structures in the image, comparing to τ. In fact, especially in the case of several structures, the surrogate method does not cause the total loss of correlation in the image even randomizing the pixel locations. A memory of the behavior of τ can still be found in τrand, thus making G less dependent on the structure level. The estimation of geometric information loss as performance index represents an important improvement in evaluation of reconstruction methods. Previous studies only focused on the RMSE of reconstruction which is only able to give information about the first-order statistics of the reconstruction error. The presence of structures in the error image can be detected only by using a second-order statistic index, like the autocorrelation. Structures may be present in the error image also in the case of
low RMSE. Therefore a performance analysis limited to the RMSE does not benefit from all the information available in the error image. Despite the lower RMSE, LS-NUFFT showed an error image more structured than GR. A comparison limited to the RMSE would not have been able to recognize that also LS-NUFFT loses relevant information in the reconstruction process. For both the simulated and experimental datasets, LS-NUFFT obtained a G value higher than GR. However, the relative difference between the G values of the two reconstruction methods is higher in the case of the simulated data. This is due to the different levels of structure and acquisition noise in the analyzed images. In fact, the acquisition noise increases the random error component, hence decreasing the value of G. The level of noise present in the reconstructed image also depends on the reconstruction method as suggested in [36]. The shape and the size of the interpolation kernel for gridding methods is a key issue in determining the accuracy of the reconstruction [11–13]. For example, the characteristics of the gridding kernel can change the amplitude of the reconstruction error or the shading on the reconstructed image. For LS-NUFFT the choice of the gridding kernel is guided by the pseudo-inverse matrix, hence minimizing the approximation error. For LS-NUFFT the size of the kernel can cause relevant reconstruction error
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
2261
Fig. 6. Experimental results shown for Pyruvate on central slice. (a) Modulus of reconstructed image using DS. Image has been normalized by its maximum value. (b and c) error image of GR method (b) and LS-NUFFT (c) using DS as reference. FOV is circular for both images. (d and e) Image of the 2D autocorrelation for the error image in GR (d) and LS-NUFFT (e) method.
related to the matrix pseudo-inverse calculation [15]. Therefore the kernel size was carefully chosen in order to make the two reconstruction methods comparable and avoid the need for a regularization method for the pseudo-inverse calculation. In this study the same weighting function was used for the data to compare the reconstructions methods. The DCF was evaluated using the Voronoi diagram for all three methods considered. DCF has a relevant influence on the reconstructed image, as proved by the extensive studies on its shape [36–39]. LS-NUFFT differs from GR only in operations coming after data weighting. The use of different DCFs would not allow distinguishing the differences in reconstructed images due to DCF, and the ones due to the interpolator. In this work the effect of the DCF on the reconstructed image was not considered. The study of the effect of DCF over the loss in geometric information could be an interesting development of the present work. The use of the 2D autocorrelation function could also be useful in parallel imaging applications [3,4]. Parallel imaging techniques commonly suffer from coherent aliasing artifacts, which could easily be detected by the 2D autocorrelation function. The presence of an image ghost would be detected with a secondary peak placed in the same spatial location of the ghost. For this application the amplitude of the autocorrelation peak, rather than the width, could be useful to quantify the intensity of the ghost.
5. Conclusions This work compared the LS-NUFFT and the GR reconstruction methods. The performance was tested on simulations with modified Shepp-Logan phantom and on hyperpolarized 13C 3D CSI in vivo data as exemplary experimental data set. Even though LSNUFFT obtained a smaller RMSE than GR, the error image of LSNUFFT appeared more structured than GR. These results strongly suggest that both reconstruction error and geometrical loss of information indices should be adopted for a complete evaluation of reconstruction methods. Even though LS-NUFFT shows better
reconstruction accuracy, it is relevant to be aware that the reconstruction process lost interesting information and not only noise.
Conflict of interest None declared. References [1] S. Hu, M. Lustig, A.P. Chen, J. Crane, A. Kerr, D.A. Kelley, R. Hurd, J. Kurhanewicz, S.J. Nelson, J.M. Pauly, D.B. Vigneron, Compressed sensing for resolution enhancement of hyperpolarized 13C flyback 3D-MRSI, J. Magn. Reson. 192 (2) (2008) 258–264. [2] M. Lustig, D. Donoho, J.M. Pauly., MRI Sparse, The application of compressed sensing for rapid MR imaging, Magn. Reson. Med. 58 (6) (2007) 1182–1195. [3] K.P. Pruessmann, M. Weiger, M.B. Scheidegger, P. Boesiger., SENSE: sensitivity encoding for fast MRI, Magn. Reson. Med. 42 (5) (1999) 952–962. [4] M.A. Griswold, P.M. Jakob, R.M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, A. Haase, Generalized autocalibrating partially parallel acquisitions (GRAPPA), Magn. Reson. Med. 47 (6) (2002) 1202–1210. [5] A.Z. Lau, A.P. Chen, N.R. Ghugre, V. Ramanan, W.W. Lam, K.A. Connelly, G. A. Wright, C.H. Cunningham, Rapid Multislice Imaging of Hyperpolarized 13C Pyruvate and Bicarbonate in the heart, Magn. Reson. Med. 64 (5) (2010) 1323–1331. [6] C.K. Turnes, J. Romberg. Spiral F.F.T.: An efficient method for 3-D FFTs on spiral mri contours. in: Proceedings of IEEE 17th International Conference on Image Processing, 2010. [7] R.D. Hoge, R.K.S. Kwan, G.B. Pike, Density compensation functions for spiral MRI, Magn. Reson. Med. 38 (1996) 117–128. [8] G.E. Sarty, Natural k-plane coordinate reconstruction method for magnetic resonance imaging: mathematical foundations, Int. J. Imag. Syst. Technol. 8 (1997) 519–528. [9] J.I. Jackson, C.H. Meyer, D.G. Nishimura, A. Macovski, Selection of a convolution function for Fourier inversion using gridding, IEEE Trans. Med. Imag. 10 (1991) 473–478. [10] J.D. O’Sullivan, A. Fast Sinc, Function gridding algorithm for fourier inversion in computer tomography, IEEE Trans. Med. Imag. 4 (4) (1985) 200–207. [11] P.J. Beatty, D.G. Nishimura, J.M. Pauly, Rapid gridding reconstruction with a minimal oversampling ratio, IEEE Trans. Med. Imag. 24 (2005) 799–808. [12] N. Nguyen, Q.H. Liu, The regular fourier matrices and nonuniform fast fourier transforms, SIAM J. Sci. Comput. 21 (1) (1999) 283–293. [13] J.A. Fessler, B.P. Sutton, Nonuniform fast fourier transforms using min-max interpolation, IEEE Trans. Sign. Proc. 51 (2) (2003) 560–574. [14] L. Sha, H. Guo, A.W. Song, An improved gridding method for spiral MRI using nonuniform fast Fourier transform, J. Magn. Res. 162 (2003) 250–258.
2262
F. Gibiino et al. / Computers in Biology and Medicine 43 (2013) 2256–2262
[15] F. Gibiino, V. Positano, L. Landini, M.F. Santarelli, Regularization techniques on non-uniform Fast Fourier Transform, Int. J. Num. Meth. Biomed. Eng. 29 (5) (2013) 561–573. [16] A. Dutt, V. Rokhlin, Fast fourier transforms for nonequispaced data, SIAM J. Sci. Comput. 14 (6) (1993) 1368–1393. [17] M. Singh, M. Mandal, A. Basu Gaussian and Laplacian of Gaussian weighted functions for robust feature based tracking. Pattern Recognition Lett. 2005;26:1995-2005. [18] G.N. Srinivasan, Shobha G., Statistical texture analysis. in: Proceedings of world academy of science, engineering and technology, 2008, pp. 2070–3740. [19] W.K. Pratt, Digital Image Processing: PISK Inside, 16, John Wiley & Sons, Inc, New York, USA (2001) 555–577. (Chapter). [20] R.M. Haralick, L. Watson, A facet model for image data, Comput. Vision Graphics Image Process. 15 (1981) 113–129. [21] D.M. Rubin, A simple autocorrelation algorithm for determining grain size from digital images of sediment, J. Sediment. Res. 74 (1) (2004) 160–165. [22] P. Henry, G. Adriany, D. Deelchand, R. Gruetter, M. Marjanska, G. Oz, E. R. Seaquist, A. Shestov, K. Ugurbil, In vivo 13C NMR spectroscopy and metabolic modelling in the brain: a practical perspective, Magn. Res. Imaging 24 (2006) 527–539. [23] K.J. Neurohr, E.J. Barrett, R.G. Shulman, In vivo carbon-13 nuclear magnetic resonance studies of heart metabolism, Proc. Natl. Acad. Sci. USA 80 (1983) 1603–1607. [24] J.H. Ardenkjaer-Larsen, B. Fridlund, A. Gram, G. Hansson, L. Hansson, M. H. Lerche, R. Servin, M. Thaning, K. Golman, Increase in signal-to-noise ratio of 410,000 times in liquid-state NMR, Proc. Natl. Acad. Sci. USA 100 (18) (2003) 10158–10163. [25] K. Golman, R. in’t Zandt, M. Thaning, Real-time metabolic imaging, Proc. Natl. Acad. Sci. USA 103 (2006) 11270–11275. [26] L. Menichetti, F. Frijia, A. Flori, F. Wiesinger, V. Lionetti, G. Giovannetti, G. D. Aquaro, F.A. Recchia, J.H. Ardenkjaer-Larsen, M.F. Santarelli, M. Lombardi, Assessment of real-time myocardial uptake and enzymatic conversion of hyperpolarized [1-13C]pyruvate in pigs with slice selective Magnetic Resonance Spectroscopy, Contrast Media Mol. Imaging 7 (1) (2012) 85–94. [27] G.D. Aquaro, F. Frijia, V. Positano, L. Menichetti, M.F. Santarelli, J.H. ArdenkjaerLarsen, F. Wiesinger, V. Lionetti, S.L. Romano, G. Bianchi, D. Neglia, G. Giovannetti, R. F. Schulte, F.A. Recchia, L. Landini, M. Lombardi, 3D CMR mapping of metabolism by
[28]
[29]
[30]
[31]
[32] [33]
[34]
[35] [36]
[37] [38] [39]
hyperpolarized 13C-Pyruvate in ischemia-reperfusion, JACC Cardiovasc. Imaging 6 (6) (2013) 743–744. C.H. Cunningham, A.P. Chen, M. Lustig, B.A. Hargreaves, J. Lupo, D. Xu, J. Kurhanewicz, R.E. Hurd, J.M. Pauly, S.J. Nelson, BD.B. Vigneron, Pulse sequence for dynamic volumetric imaging of hyperpolarized metabolic products, J. Magn. Res. 193 (2008) 139–146. V. Rasche, R. Proksa, R. Sinkus, P. Börnert, H. Eggers, Resampling of data between arbitrary grids using convolution interpolation, IEEE Trans. Med. Imaging 18 (5) (1999) 385–392. D. Moratal, A.V. Lluch, V. Bodí, M.E. Brummer, Magnetic resonance imaging gridding reconstruction methods with and without density compensation functions, IEEE Lat. Am. Trans. 9 (1) (2011) 774–778. M.S. Sussman, J.A. Stainsby, N. Robert, N. Merchant, G.A. Wright, Variabledensity adaptive imaging for high-resolution coronary artery MRI, Magn. Reson. Med. 48 (2002) 753–764. J.G. Pipe, Reconstructing MR images from undersampled data: data-weighting considerations, Magn. Reson. Med. 43 (6) (2000) 867–875. E.K. Brodsky, J.H. Holmes, H. Yu, S.B. Reeder, Generalized k-Space decomposition with chemical shift correction for Non-Cartesian water-fat imaging, Mag. Reson. Med. 59 (2008) 1151–1164. F. Wiesinger, E. Weidl, M.I. Menzel, M.A. Janich, O. Khegai, S.J. Glaser, A. Haase, M. Schwaiger, R.F. Schulte, IDEAL Spiral CSI for dynamic metabolic mr imaging of hyperpolarized [1-13C], Pyruvate. Magn. Reson. Med. 68 (2011) 8–16. J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, Testing for nonlinearity in time series: the method of surrogate data, Physica D 58 (1992) 77–94. F. Gibiino, V. Positano, G. Giovannetti, F. Frijia, L. Menichetti, J.H. ArdenkjaerLarsen, F. Wiesinger, R. Schulte, D. De Marchi, V. Lionetti, G. Aquaro, M. Lombardi, L. Landini, M.F. Santarelli, Reconstruction methods from hyperpolarized 13C chemical shift imaging spiral 3D data: comparison between direct summation and gridding method, ISBI proceeding (2012) 614–617. J.G. Pipe, P. Menon, Sampling density compensation in mri: rationale and an iterative numerical solution, Magn. Reson. Med. 41 (1999) 179–186. D. Rosenfeld, An optimal and efficient new gridding algorithm using singular value decomposition, Magn. Reson. Med. 40 (1998) 14–23. H. Sedarat, D.G. Nishimura, On the optimality of the gridding reconstruction algorithm, IEEE Trans. Med. Imaging 19 (4) (2000) 306–317.
