Matrix pencil decomposition of time-resolved proton MRI for robust and improved assessment of lung ventilation and perfusion
Grzegorz Bauman1,2 and Oliver Bieri1,2

1Radiological Physics, University Hospital of Basel, Basel, Switzerland, 2Department of Biomedical Engineering, University of Basel, Basel, Switzerland


In the contemporary Fourier decomposition lung MRI, time-resolved registered 2D image series are Fourier transformed to identify in a power spectrum the underlying respiratory and cardiac frequencies. Subsequently, the amplitudes corresponding to the respiratory and cardiac motion are extracted voxel-wise to eventually produce ventilation and perfusion images. However, the analysis of truncated oscillatory signals and the peak search in the Fourier spectrum is usually very unstable and inaccurate. Here, we propose to use a robust and fully-automated method of signal analysis using a matrix pencil decomposition in combination with a linear least squares analysis for improved quantitative pulmonary function assessment.


Fourier decomposition (FD)1 MRI provides regional ventilation and perfusion information from a single free-brething acquisition series. In the current FD MRI implementation, the ventilation and perfusion images are produced by integration of frequencies corresponding to respiratory and cardiac cycles in a power spectrum. For proper quantification of FD MRI, a precise estimation of the underlying frequencies is mandatory; however, the peak search in the Fourier spectrum is usually very unstable, inaccurate and slow. Furthermore, the analysis of oscillatory truncated signals using simple fast Fourier transform often results in spectral leakage and Gibbs ringing. In this work we propose to use matrix pencil (MP) method2 to get a proper estimation of the respiratory and cardiac frequency components in combination with a linear least squares approach to find the corresponding amplitudes.


As for the contemporary FD, the MP analysis is performed voxel-wise on registered image series i.e. with lung structures aligned along the time axis. MP method has been already suggested as an efficient tool to estimate the spectral components for a sum of $$$M$$$ damped exponentials of form:


where: $$$s_n$$$- NMR signal in a voxel, $$$n=0...N-1$$$- equidistant data points, $$$\eta_n$$$- noise-term, $$$a_p$$$- complex amplitude, $$$\omega_p$$$- frequency, and $$$R_p$$$-damping factor of $$$M$$$ distinct exponentials with signal poles $$$z_p{\equiv}e^{(i\omega_p-R_p)\Delta{t}}$$$ (see Ref.2).

For MP the pulmonary signal is first summed over the whole image and subsequently low and high-pass filtered to separate cardiac and respiratory cycles. The poles corresponding to the respiratory and cardiac cycles are then retrieved from MP decomposition of low-pass filtered $$$S_r(n)$$$ and high-pass filtered $$$S_c(n)$$$ signals:


where $$$u$$$ and $$$v$$$ indicate number of detected respiratory and cardiac frequency components. Note that generally, one pole for the low-pass filtered signal will have $$$\omega\sim{0}$$$ and will thus retrieve the local signal offset. From the corresponding respiratory and cardiac poles, two matrices are defined:


Subsequently, the complex respiratory and cardiac amplitudes can now be obtained voxel-wise by the Moore-Pensore pseude-inverse of the matrix $$$Z$$$ and the measured local signal time-course $$$s_n(x,y)$$$:

$$\textbf{a}_{r,c}(x,y) = \textbf{Z}^{+}s(x,y)$$

Finally, summation yields the total respiratory and cardiac amplitudes:



MR data acquisition

All measurements were performed on a 1.5T MR-scanner (Avanto, Siemens Healthcare). Four healthy volunteers were scanned during free-breathing with a time-resolved 2D ultra-fast steady-state free precession (ufSSFP) sequence3 using the following parameters: field-of-view=450x450mm2, slice thickness=12mm, TE/TR/TA=0.67ms/1.46ms/119ms, flip angle α=65°, bandwidth=2056Hz/pixel, matrix=128x128 (interpolated to 256x256), 200 coronal images, 3.33 images/s, parallel imaging GRAPPA factor=2.

Image post-processing

Respiratory motion in ufSSFP data was compensated with a non-rigid image registration algorithm4. Subsequently, FD and MP decomposition of the registered time-series were performed voxel-wise to estimate the respiratory $$$(A_r)$$$ and cardiac $$$(A_c)$$$ amplitudes and create qualitative ventilation- and perfusion-weighted images. Eventually, quantitative analysis of fractional ventilation5, perfusion6, and blood arrival time7 was performed. The analysis was accomplished for time-series of different length $$$N{\in}{140...190}$$$.


Illustrative qualitative ventilation- and perfusion-weighted images obtained in a volunteer using MP and the contemporary FD for time series containing 190 images are shown in Figs.1a-1d. The variability in $$$A_r$$$ and $$$A_c$$$ as a function of time series length measured in segmented lung area (see Figure 1e) is shown in Figure 1f. The difference in maximal and minimal amplitudes $$$A_r$$$ and $$$A_c$$$ can reach up to 32% and 24% with FD, respectively; while for the MP method the difference was 5% and 12% for $$$A_r$$$ and $$$A_c$$$, respectively. Figure 2 shows exemplar quantitative maps of fractional ventilation, perfusion and blood arrival time obtained in a healthy volunteer and generated using either MP or FD for time-series length $$$N=154$$$. Furthermore, maps of differences between the MP and FD quantitative images of fractional ventilation, perfusion and blood arrival time were shown in the bottom row in Figure 2.

Discussion and Conclusion

In this work, we presented an improved method of regional pulmonary function assessment using contrast-agent free time-resolved ufSSFP imaging. The proposed method employs a spectral analysis based on a matrix pencil approach and linearized least squares fitting. The method is shown to allow for an fully automatic and improved estimation of amplitudes of respiratory and cardiac signal modulation in the lung parenchyma. We have shown that the standard implementation of FD MRI based on the straightforward FFT is due to its increased variability in amplitude estimation not well-suited for quantitative analysis. However, the improved signal analysis based on the MP decomposition method may result in increased accuracy and reproducibility of quantitative parameters derived from the time-resolved chest acquisitions such as regional ventilation, perfusion and blood arrival time.


Part of this work was supported by Siemens Healthcare.


1. Bauman G, Puderbach M, Deimling M et al. Non-contrast-enhanced perfusion and ventilationassessment of the human lung by means of fourier decomposition in proton MRI.Magn Reson Med. 2009 Sep;62(3):656-64.

2. Lin Y, Hodgkinson P, Ernst M, Pines A. A Novel Detection-Estimation Scheme for Noisy NMR Signals: Applications to Delayed Acquisition Data. J Magn Res 1997; 128:30-41.

3. Bieri O. Ultra-fast steady state free precession and its application to in vivo (1) H morphological and functional lung imaging at 1.5Tesla. Magn Reson Med. 2013 Sep;70(3):657-63.

4. Chefd’hotel C, Hermosillo G, Faugeras O. Flows of diffeomorphismsfor multimodal image registration. In: Proceedings of the IEEE International Symposium on Biomedical Imaging (ISBI’2002), Washington, DC, USA, July 2002. pp. 753–756.

5. Zapke M, Topf HG, Zenker M et al. Magnetic resonance lung function--abreakthrough for lung imaging and functional assessment? A phantom study andclinical trial. Respir Res. 2006 Aug 6;7:106.

6. Kjørstad Å, Corteville DM, Fischer A et al. Quantitative lung perfusion evaluation using Fourier decompositionperfusion MRI. Magn Reson Med. 2014 Aug;72(2):558-62.

7. Bauman G, Eichinger M, Uecker M. High temporal resolution radial bSSFP sequence with nonlinear inverse reconstruction for the measurement of the pulmonary blood inflow time using Fourier decomposition MRI. Proceedings of the 20th annual meeting of the ISMRM, Melbourne, Australia, 2012:1955.


Figure 1. (a) Native coronal ufSSFP image. The placement of the ROI used for analysis of the time series is indicated by the yellow outline. (b) Dependence between the mean respiratory ($$$A_r$$$) and cardiac ($$$A_c$$$) amplitudes calculated using MP and FD in the ROI as a function of the time series length $$$N$$$. MP shows lower variability than FD in estimation of both amplitudes.

Figure 2. Quantitative maps of fractional ventilation and perfusion in a healthy volunteer calculated using MP and contemporary FD method for the image data set length of 154. In the bottom row the images show differences between the maps obtained using MP and FD. Please note that the scales were adapted for better visualization of the differences.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)