A multi-modality approach for enhancing 4D flow magnetic resonance imaging via sparse representation

2.1. Flow reconstruction via library-based sparse representation

The MSR reconstruction assumes that the measurement process of 4D flow MRI can be modelled as:

where $\mathit{U}\in {\mathbb{R}}^N$ represents the underlying velocity field on a fine grid ${\mathcal{G}}_{\mathrm{fine}}$ with N grid points, $\mathit{u}\in {\mathbb{R}}^n$ contains the 4D flow MRI measured velocity values at n voxel-centre locations, C is the measurement matrix with a shape of n × N which transforms U on ${\mathcal{G}}_{\mathrm{fine}}$ to the coarse MRI grid ${\mathcal{G}}_{\mathrm{MRI}}$ and $ϵ\in {\mathbb{R}}^n$ represents the measurement noise. The vectors U and u contain only one velocity component, and the measurement process modelled by (2.1) applies to each of the three velocity components acquired using 4D flow MRI. In the present study, each velocity-component field was reconstructed individually such that the optimization problem has a lower rank for enabling a sparse representation and a lower dimension to save computational cost. The smoothing effect of Cartesian 4D flow MRI due to finite k-space coverage was modelled by convolving the underlying velocity field with a truncated sinc-function kernel $\mathcal{K}$ [30] aswhere Δx, Δy and Δz correspond to the grid sizes of 4D flow MRI. Therefore, the coefficients in the measurement matrix C were assigned based on the kernel function values at the grid nodes related to each MRI voxel. The row number of the matrix C was used to index the MRI voxels, and the columns were linked to the grid points of ${\mathcal{G}}_{\mathrm{fine}}$. The coefficient at the ith row and jth column of C was given as the kernel function value based on the distance between the ith voxel centre and the jth grid point normalized by the sum of the kernel function values from all grid points related to the ith voxel centre.

Using the high-resolution flow data acquired from CFD simulations and in vitro PTV measurements, a flow library $\{{\mathit{\Psi }}_{\mathit{k}}\}$ with ${\mathit{\Psi }}_{\mathit{k}}\in {\mathbb{R}}^N$ and k = 1, 2, … , m can be constructed. The library components ${\mathit{\Psi }}_{\mathit{k}}$ can be generated as a snapshot [28] or the mode extracted from the high-resolution flow data [25,26,29,31] on the grid ${\mathcal{G}}_{\mathrm{fine}}$. If the flow library contains an extensive collection of representative examples of the probable flow structures in the same aneurysm geometries, U can be accurately expressed as a linear combination of the library components as:

where $\mathit{s}\in {\mathbb{R}}^m$ is the vector of library coefficients and Ψ contains the library components Ψ_k as its columns. To find the coefficient vector $\hat{\mathit{s}}$ such that the MSR-reconstructed velocity field ($\hat{\mathit{U}}$) obtained as $\hat{\mathit{U}}=\mathit{\Psi }\hat{\mathit{s}}$ accurately reproduces the underlying velocity field U, the LASSO regression [32] was used in the MSR flow reconstruction,where $\Vert \mathit{u}-C\mathit{\Psi }\mathit{s}{\Vert }_2$ represents the L2-norm of the fitting residuals and α > 0 is the parameter that controls the strength of the regularization term $\Vert \alpha \mathit{s}{\Vert }_1$. The hyperparameter α was optimized iteratively using fivefold cross-validation, and the optimization of (2.4) was carried out using the Python scikit-learn subroutine ‘LassoCV' [33,34]. The penalization of $\Vert \mathit{s}{\Vert }_1$ in (2.4) promotes a sparse representation, i.e. most of the coefficients in $\hat{\mathit{s}}$ are zeros. However, minimizing $\Vert \mathit{s}{\Vert }_1$ also tends to reduce the amplitude of the MSR-reconstructed field since $\hspace{0.17em}\hat{U}=\mathit{\Psi }\hat{\mathit{s}}$. Thus, the obtained flow field is rescaled to keep the same L1-norm between the measured velocity data u and the reproduced velocity measurement $\hat{\mathit{u}}=C\mathit{\Psi }\hat{\mathit{s}}$ as

The procedure of the MSR reconstruction is shown schematically in figure 1a.

• Download figure
• Open in new tab
• Download PowerPoint

For complex flow fields such as the aneurysmal flow, it may not be feasible to construct a flow library containing representative examples of all the probable global flow fields. However, local flow regions may have a lower rank and enable sparse representation as demonstrated in previous studies [28]. As shown in figure 1b, the aneurysmal flow fields were divided into subdomains whose centres were located on Cartesian grids with resolutions h_sub of two times the MRI resolutions for the BT and ICA aneurysms. Each subdomain contains a region within 2h_sub from the centre r_sub and overlaps with neighbouring subdomains. The subdomains located near the wall only contain the regions within the blood flow. The MSR flow reconstruction was performed in each subdomain, and the global velocity field was subsequently constructed as the weighted superposition of the local reconstructions as

where $\hat{U}(\mathit{r})$ represents the combined velocity value at spatial location r, K is the total number of subdomains that contain r, ${\hat{U}}_k(\mathit{r})$ means the velocity value reconstructed in the kth subdomain whose centre locates at r_sub and w_k is the weight based on the spatial distance between r and r_sub.

2.2. Multi-modality flow data acquisition and flow library construction

In vivo flow data in a BT aneurysm and an ICA aneurysm were acquired using 4D flow MRI on a 3 T scanner (Skyra; Siemens Healthcare, Erlangen, Germany). MRI acquisition of the BT aneurysm was performed at the San Francisco VA Medical Center, and the ICA aneurysm was imaged at Northwestern Memorial Hospital. Gadolinium contrast was used for the BT aneurysm, whereas no contrast was used for the ICA aneurysm. The sac diameter was 7 mm and 10 mm for the BT aneurysm and ICA aneurysm, respectively. The in vivo flow data were on Cartesian grids with a spatial resolution of 1.25 × 1.25 × 1.33 mm³ for the BT aneurysm and 1.09 × 1.09 × 1.30 mm³ for the ICA aneurysm, resulting in about 5 and 9 voxels across the diameter of the aneurysmal sac for the BT aneurysm and ICA aneurysm, respectively. The temporal resolution Δt was 40.5 ms (20 frames per cycle) and 44.8 ms (13 frames per cycle) for the BT and ICA aneurysms, respectively. The encoding velocity sensitivity (venc) was 100 cm s⁻¹ for the BT aneurysm and 120 cm s⁻¹ for the ICA aneurysm. Ethical approval of all experimental procedures and protocols was granted by the institutional review boards at Purdue University, Northwestern Memorial Hospital and San Francisco VA Medical Center.

In vitro flow measurements were done using PTV in 1 : 1 scaled silicone models of the two aneurysms with blood-mimicking fluids composed of water–glycerol–urea [35]. The flow phantom for each aneurysm was fabricated as follows: based on the segmented vessel geometry, a positive-space model was 3D printed and embedded into a tear-resistant silicone block. The model was then cut from the block and replaced by a low melting point metal (Cerrobend 158; bismuth alloy). The block was cut away and embedded in a block of optically clear polydimethylsiloxane silicone (PDMS-Slygard 184). The metal was then melted out after the PDMS was hardened. The time-dependent inflow waveform was driven by a computer-controlled gear pump and designed to match the waveform obtained from the in vivo 4D flow MRI data. The acquired particle images were processed using DaVis 10.0 (LaVision Inc.) with the shake-the-box method. The temporal resolution was 2.5 ms for the BT aneurysm and 1.5 ms for the ICA aneurysm. CFD simulations were performed using FLUENT 18.1 (ANSYS) for the two aneurysms. An unstructured tetrahedral mesh generated using HyperMesh 14.0 (Altair Engineering, Troy, MI, USA) was employed with a nominal cell size of 0.15 mm. The flow-rate waveforms obtained from the in vivo 4D flow MRI were used as the boundary conditions for the CFD models. A Newtonian fluid and rigid wall were assumed in the CFD models. The rigid wall was assumed because the previous investigation with cine-MRI showed no observable wall movement of intracranial vessels over the cardiac cycle [8], and the aneurysm disease can reduce the elasticity of the arterial wall [36]. For each aneurysm, three cardiac cycles were simulated with a temporal resolution of 1.5 ms, and the results of the last cycle were used in the present study. More details on the in vivo imaging, in vitro measurement and CFD models can be found in [15].

The unstructured PTV and CFD data were linearly interpolated to Cartesian grids ${\mathcal{G}}_{\mathrm{fine}}$ with isotropic resolutions of 0.3 and 0.4 mm for the BT and ICA aneurysms, respectively, with more than 20 grid points across the diameter of each aneurysmal sac. Instead of the basis from POD [25,26], the velocity fields were directly used as the library components. For each CA geometry and each velocity component, the flow library was constructed as a collection of 100 time frames randomly selected from the high-resolution PTV and CFD data on ${\mathcal{G}}_{\mathrm{fine}}$. Different flow libraries were constructed with different numbers of time frames from CFD and PTV data. The flow libraries with an equal number (50 each) of PTV and CFD components were referred to as the balanced flow libraries. The flow libraries with more PTV components than CFD were PTV dominant, while the flow libraries with more components from CFD than PTV were CFD dominant.

2.3. Synthetic 4D flow MRI data generation

To evaluate MSR's performance, synthetic 4D flow MRI datasets were generated from the high-resolution PTV and CFD datasets. For each velocity component, the complex-valued signal was created as:

where M_u means the signal created for the u velocity component from the high-resolution data (u_fine) and I is the signal magnitude which was set as 1.0 inside the blood flow while 0.2 outside the blood flow. The signal of the 4D flow MRI data was generated by spatially and temporally smoothing M_u aswhere * denotes convolution, $\mathcal{K}$ is the spatial smoothing kernel function defined as (2.2) and $\mathcal{T}$ represents the temporal smoothing kernel function with Δt being the temporal resolution of the MRI data. Following a four-point reference method, one reference signal (M₀) was generated from a zero field. The phase data ψ for each velocity component were generated from the flow-sensitive and the reference signals, e.g.where ψ_u is the phase for the u velocity component, $M_0^{\ast }$ is the complex conjugate of M₀ and angle() means calculating the angle (between −π and π) of a complex value. The velocity data were obtained from the phase data as

For each CA, two synthetic 4D flow MRI datasets with the same spatio-temporal resolutions as the in vivo measurements, named CFD-synMRI and PTV-synMRI, were created from the CFD data and the PTV data, respectively. To test the robustness of MSR, normally distributed noise was added to the complex-valued signal with the standard deviation defined as

where SNR_I is the signal-to-noise ratio, $\overline{V}$ is the mean velocity magnitude and VNR is the velocity-to-noise ratio [37]. The reciprocal of the VNR is defined as the noise level of the synthetic dataset, and four noise levels (5%, 10%, 15% and 20%) were considered in the present study. To test the performance of MSR on different MRI grid resolutions, additional synthetic MRI datasets were created with grid sizes varying from two to five times the grid sizes of the CFD and PTV data on ${\mathcal{G}}_{\mathrm{fine}}$.

2.4. Methods for velocity and haemodynamic analysis

The CFD and PTV velocity fields at the corresponding time frames of the synthetic 4D flow MRI data were considered as the ground truth for CFD-synMRI and PTV-synMRI, respectively. To assess the accuracy of the synthetic MRI and the MSR-reconstructed velocity fields, velocity error was determined as the difference from the ground truth, and the velocity error level for each dataset was quantified as the root-mean-square (RMS) of the velocity error magnitudes normalized by the RMS of the velocity magnitudes from the ground truth. It should be noted that the velocity fields from the ground truth time frames were excluded from the flow library for the reconstruction of CFD-synMRI and PTV-synMRI data such that the ‘ground truth' was not embedded in the flow library.

To study the effects of the MSR reconstruction on the flow-derived haemodynamic quantities, pressure and WSS were computed from the MRI and the MSR-reconstructed velocity fields. The pressure reconstruction was carried out using the measurement-error-based weighted least-squares method [14]. A pressure of 0 Pa was assigned at the inflow locations of the aneurysmal sac as the reference pressure for both CAs. The WSS was calculated from the near-wall velocity using thin-plate spline radial basis functions [15], and the time-averaged WSS (TAWSS) was subsequently obtained by averaging the WSS magnitudes across the full cardiac cycle. The pressure and WSS evaluated from the high-resolution ground-truth velocity fields were employed as the ground truth for the error analysis of the pressure and WSS from the synthetic MRI cases.

3.1. Multi-modality velocity fields

The velocity fields at peak systole for all three modalities and the synthetic MRI data are presented in figure 2 for both aneurysms. For the BT aneurysm, the basilar artery's flow swirled in the aneurysmal sac and exited mainly through the posterior cerebral arteries (PCAs). For the ICA aneurysm, the flow entered from the ICA, swirled in the aneurysmal sac and exited through the distal ICA and middle cerebral artery (MCA). In the BT aneurysmal sac, the maximum descending velocity was 0.4 m s⁻¹ from MRI and CFD, while it was 0.3 m s⁻¹ from PTV. In the ICA aneurysm, a stronger inflow towards the aneurysmal sac was observed from CFD than from the other two modalities, with the maximum velocity around the ‘neck' of the aneurysmal sac being 0.55 m s⁻¹ from CFD, 0.35 m s⁻¹ from MRI and 0.3 m s⁻¹ from PTV. Compared with the corresponding CFD/PTV data, the flow structures observed from the synthetic MRI were similar but were under-resolved in the near-wall regions and smaller vessels.

• Download figure
• Open in new tab
• Download PowerPoint

3.2. Reconstruction of synthetic 4D flow MRI

The velocity error level for each synthetic MRI dataset is presented in figure 3a as a function of the noise levels. Even without noise, the error level of synthetic MRI velocity data was more than 0.27 owing to the smoothing effect. As the noise level increased from 0% to 20%, the velocity error level increased by 0.06 to 0.1. The error levels of the MSR-reconstructed velocity fields were lower than 0.15 for all the cases, as shown in figure 3a, leading to an error reduction of more than 70% for each case. The increase in the MSR-reconstructed velocity's error level due to the increase in the noise level was less than 0.04. The velocity error levels as functions of the spatial resolution of the synthetic MRI data are presented in figure 3b. As the MRI resolution increased from two to five times the MSR-reconstructed resolution h_MSR, the velocity error level of the synthetic MRI increased by more than 0.1 because of the greater spatial smoothing effect, while the error level of the MSR-reconstructed velocity fields increased by only 0.03 and did not exceed 0.13.

• Download figure
• Open in new tab
• Download PowerPoint

To investigate the effect of the flow library on the MSR's performance, flow libraries consisting of different numbers of components from CFD and PTV data were used for the reconstruction of the synthetic MRI data. The MSR-reconstructed velocity error level as a function of the library composition is shown in figure 4a for the synthetic datasets with the in vivo MRI resolution and 10% noise. The library composition is represented by the fraction of library components from the CFD data, with the rest of the PTV data components, and the total number of components was 100 for all the flow libraries. For the reconstruction of CFD-synMRI datasets, the velocity error levels were greater than 0.4 if the flow library did not contain any components from CFD data. By contrast, the velocity error levels were lower than 0.1 if there were at least 10 library components from CFD. Similarly, the reconstructions of PTV-synMRI datasets resulted in velocity error levels greater than 0.45 with libraries containing only CFD components, while the velocity error level was less than 0.18 with at least 10 PTV components. For the MSR-reconstructed flow fields, the relative contribution from the CFD library components was determined as the L1-norm of the coefficients in $\hat{s}$ corresponding to the CFD library components normalized by the total L1-norm of $\hat{s}$, with the rest of the contribution from PTV library components. Figure 4b presents the relative contribution from the CFD library components for the reconstructions of synthetic MRI datasets as a function of the library composition. With the flow libraries containing a mix of CFD and PTV components, the reconstructions of CFD-synMRI datasets yielded relative contributions from CFD of more than 0.7, while the reconstructions of PTV-synMRI datasets yielded less than 0.3 on the CFD components. Figure 4c shows the sparsity of the MSR reconstruction defined as the number of zero-valued coefficients divided by the total number of coefficients in $\hat{s}$ as a function of the library composition. The reconstructions of PTV-synMRI datasets with flow libraries containing only CFD components yielded sparsity less than 0.8, and the reconstructions of CFD-synMRI datasets with flow libraries containing only PTV components had a sparsity less than 0.78. The reconstructions from other cases showed relative sparsity of above 0.85.

• Download figure
• Open in new tab
• Download PowerPoint

To assess the accuracy of the flow-derived haemodynamic quantities, Bland–Altman analysis was performed to compare the pressure and WSS evaluated from the synthetic MRI and from the MSR-reconstructed data with the ground truth as shown in figure 5. In addition, the statistical distributions of these haemodynamic parameters are compared in figure 5 using probability density functions (PDFs) with the mean values indicated using vertical lines. Compared with the synthetic MRI data created from the CFD and PTV data, the haemodynamic quantities derived from the MSR-reconstructed fields were more consistent with the ground truth, as suggested by the lower bias and STDs of the differences in the Bland–Altman plots. As suggested in figure 5, the WSS derived from synthetic MRI data also showed greater underestimation with greater mean WSS values, while the WSS from the MSR-reconstructed data did not have this correlation. Moreover, the PDFs of the WSS and pressure from MSR-reconstructed fields closely resembled those of the ground truth. The mean of the WSS from the synthetic MRI data was underestimated by 40–60% compared with the ground truth.

• Download figure
• Open in new tab
• Download PowerPoint

3.3. Reconstruction of in vivo 4D flow MRI

The MSR flow reconstruction was applied to the in vivo 4D flow MRI data acquired for the same two CAs using the balanced flow libraries with 50 CFD components and 50 PTV components. The flow fields at peak systole, pressure fields at peak systole and TAWSS obtained from the in vivo MRI and from the MSR-reconstructed flow fields are presented in figure 6a,b for the BT and ICA aneurysms, respectively. For both aneurysms, the MSR-reconstructed flow structures in the aneurysmal sacs were qualitatively similar to the in vivo MRI data, while the flows in the near-wall region and small vessels were better resolved by the MSR reconstruction. Compared with the in vivo MRI, the MSR-reconstructed flow data yielded more distinct low-pressure regions around the centre of the aneurysmal sac, as clearly shown on the lateral view. For the BT aneurysm, the minimum relative pressure observed around the core of the aneurysmal sac was 0 Pa and −25 Pa for the in vivo MRI and MSR-reconstructed flow data, respectively. For the ICA aneurysm, the relative pressure was −5 Pa and −10 Pa around the centre of the aneurysmal sac for the in vivo MRI and MSR-reconstructed flow data, respectively. The TAWSS obtained from the MSR-reconstructed fields were higher than the TAWSS from in vivo MRI. The TAWSS of the BT aneurysm calculated from the MSR-reconstructed data also had different spatial distributions compared with the TAWSS from the MRI data, as suggested by the anterior and superior views. For example, regions with TAWSS as high as 4 Pa were observed from the anterior and superior views of the MSR-reconstructed data for the BT aneurysm, which were absent from the in vivo MRI data. For the ICA aneurysm, the high TAWSS region (greater than 3 Pa) around the neck of the aneurysmal sac was smaller for the in vivo MRI data than for the MSR-reconstructed data.

• Download figure
• Open in new tab
• Download PowerPoint

The time-dependent median and the interquartile range of velocity magnitudes, WSS and pressure obtained from the MRI and MSR-reconstructed data are presented in figure 7 at all time points within the cardiac cycle. As suggested by the median values, the velocity and WSS magnitudes of the MSR-reconstructed fields were greater at around peak systole (at t/T ≈ 0.7 for the BT aneurysm and t/T ≈ 0.4 for the ICA aneurysm), while preserved similar waveforms with MRI. The MSR reconstruction increased the mean velocity magnitudes by 13% and 6% compared with the MRI data for the BT and ICA aneurysms, respectively, while the mean WSS were increased by 60% and 51% for the BT and ICA aneurysms, respectively. As a reference, the mean WSS from CFD was 39% and 61% higher than the MRI data for the BT and ICA aneurysms, respectively, while the mean WSS from PTV was 47% and 60% higher than the MRI data for the BT and ICA aneurysms, respectively. As shown in figure 7, the MSR reconstruction altered the pressure waveforms and reduced the median pressure at peak systole by 30 Pa and 5 Pa for the BT and ICA aneurysms, respectively. The reduction in the median pressure was due to the improvement in the spatial resolution by the MSR reconstruction, which allowed the pressure estimation to better resolve the low-pressure region around the centre of the aneurysmal sac corresponding to the core of the vortical flow structure, as shown in figure 6. The contributions of the CFD and PTV library components were also obtained for the reconstructions of the in vivo 4D flow MRI data. The two-dimensional histograms in figure 8 show the distributions of the relative contributions of the CFD library components in all subdomains as a function of time. The relative contribution from CFD library components varied from 20% to 70% for the BT aneurysm, while it was 40–80% for the ICA aneurysm. The sparsity of the reconstruction was 0.88 for the BT aneurysm and 0.8 for the ICA aneurysm.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint