Matrix pencil decomposition of time-resolved proton MRI for robust and improved assessment of lung ventilation and perfusion

Grzegorz Bauman^{1,2} and Oliver Bieri^{1,2}

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) method^{2} 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:

$$s_n=\sum^{M}_{p=1}a_{p}e^{(i\omega_p-R_p)}t_n+\eta_n=\sum^{M}_{p=1}a_{p}z^n_p+\eta_n$$

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:

$$S_r(n)=\sum^{u}_{p=1}a_{r,p}z^n_{r,p},{\quad}S_c(n)=\sum^{v}_{p=1}a_{c,p}z^n_{c,p}$$

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:

$$\textbf{Z}_r=\begin{bmatrix}z_{r,1}&z_{r_2}&{\cdots}&z_{r,u}\\z_{r,1}^2&z_{r_2}^2&{\cdots}&z_{r,u}^2\\{\vdots}&\vdots&{\ddots}&\vdots\\z_{r,1}^N&z_{r_2}^N&{\cdots}&z_{r,u}^N\end{bmatrix},{\quad}\textbf{Z}_c=\begin{bmatrix}z_{c,1}&z_{c_2}&{\cdots}&z_{c,u}\\z_{c,1}^2&z_{c_2}^2&{\cdots}&z_{c,u}^2\\{\vdots}&\vdots&{\ddots}&\vdots\\z_{c,1}^N&z_{c_2}^N&{\cdots}&z_{c,u}^N\end{bmatrix}$$

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:

$$A_r(x,y)=\sum_{p=1}^{u}a_{r,p}(x,y),{\quad}A_c(x,y)=\sum_{p=1}^{v}a_{c,p}(x,y)$$

**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) sequence^{3} using the following parameters: field-of-view=450x450mm^{2}, 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 algorithm^{4}. 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 ventilation^{5}, perfusion^{6}, and blood arrival time^{7} 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.

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)

1606