Non-invasive hemodynamic analysis for aortic regurgitation using computational fluid dynamics and deep learning

by   Derek Long, et al.
The University of Auckland

Changes in cardiovascular hemodynamics are closely related to the development of aortic regurgitation (AR), a type of valvular heart disease. Pressure gradients derived from blood flows are used to indicate AR onset and evaluate its severity. These metrics can be non-invasively obtained using four-dimensional (4D) flow magnetic resonance imaging (MRI), where accuracy is primarily dependent on spatial resolution. However, insufficient resolution often results from limitations in 4D flow MRI and complex AR hemodynamics. To address this, computational fluid dynamics simulations were transformed into synthetic 4D flow MRI data and used to train a variety of neural networks. These networks generated super resolution, full-field phase images with an upsample factor of 4. Results showed decreased velocity error, high structural similarity scores, and improved learning capabilities from previous work. Further validation was performed on two sets of in-vivo 4D flow MRI data and demonstrated success in de-noising flow images. This approach presents an opportunity to comprehensively analyse AR hemodynamics in a non-invasive manner.



page 8

page 9

page 10

page 13

page 18

page 19


4DFlowNet: Super-Resolution 4D Flow MRI using Deep Learning and Computational Fluid Dynamics

4D-flow magnetic resonance imaging (MRI) is an emerging imaging techniqu...

Accurate quantification of blood flow wall shear stress using simulation-based imaging: a synthetic, comparative study

Simulation-based imaging (SBI) is a blood flow imaging technique that op...

Post-processing techniques of 4D flow MRI: velocity and wall shear stress

As the original velocity field obtained from four-dimensional (4D) flow ...

Deep Learning Super-Resolution Enables Rapid Simultaneous Morphological and Quantitative Magnetic Resonance Imaging

Obtaining magnetic resonance images (MRI) with high resolution and gener...

Learning Super-resolution 3D Segmentation of Plant Root MRI Images from Few Examples

Analyzing plant roots is crucial to understand plant performance in diff...

Resistance-Time Co-Modulated PointNet for Temporal Super-Resolution Simulation of Blood Vessel Flows

In this paper, a novel deep learning framework is proposed for temporal ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Aortic valve regurgitation, or aortic regurgitation (AR), is a common type of valvular heart disease where the aortic valve (AV) does not close properly, causing reflux of blood from the aorta into the left ventricle (LV) Valvular-Heart-Disease. This reflux of blood is known as the regurgitant jet. Diagnosis and severity of AR is determined by LV pressure gradients, which can be calculated either invasively or non-invasively. However, the non-invasive approach is known to overestimate pressure gradients and thus may lead to considerable error in determining AR severity.

Cardiovascular four-dimensional (4D) flow magnetic resonance imaging (MRI) is a novel imaging technique to quantify full-field blood flow velocities, providing a three-dimensional (3D) velocity field across a region of interest throughout the cardiac cycle. Currently, due to the small width of the regurgitant jet and the limited spatiotemporal resolution of 4D flow MRI, it fails to accurately capture the complex hemodynamics of AR and estimate LV pressure gradients. The combination of computational fluid dynamics (CFD) and deep learning with 4D flow MRI will help to generate higher resolution images and recover hemodynamic parameters lost in current MRI images. This work will extend what has already been completed in 4DFlowNet

4DFlowNet by using a wider range of flow characteristics to mimic AR, improving the data augmentation steps, and enhancing the artificial neural network with newer architecture structures.

2 Background

2.1 Clinical Background

The heart has four chambers with two separate ones for pumping. Oxygen-poor blood from the body enters the right atrium and is pumped to the right ventricle (RV) through the tricuspid valve. The RV then pumps blood to the lungs through the pulmonary valve and gas exchange occurs. The left atrium receives this blood and pumps it to the left ventricle (LV) through the mitral valve. Finally, the LV pumps oxygenated blood through the AV to the rest of the body. AR occurs when the AV does not close properly, causing blood to flow back into the LV. This forces the heart to work harder and pump more blood to the aorta, which can cause further heart problems in the future. AR also has varying levels of intensity, with most cases being trace or mild and moderate or severe cases being more rare Valvular-Heart-Disease. Acute AR is considered a medical emergency as it can cause severe pulmonary edema and hypotension, that is, excess fluid in the lungs and low blood pressure, respectively. Patients with AR are monitored yearly with echocardiography to decide whether replacement of the AV is necessary Outcomes-after-aortic-valve. LV pressure measurements are among the most important clinical parameters to assess the heart’s performance. Currently, cardiac catheterisation is the most accurate method to record LV pressure. However, this is an invasive procedure, as it involves threading a catheter through the blood vessels to the heart.

In the clinic, pressure gradients through the AV, which gives an estimate of the LV pressure, are calculated non-invasively using 2D velocities from 2D Doppler echocardiography Valvular-Heart-Disease. Due to limited information available in 2D, this method is known to overestimate pressure gradients.

2.2 4D Flow MRI

4D flow MRI is an emerging imaging technique that captures the temporal changes of 3D blood flow patterns within individual vascular structures 4DFLOWMRI; Itatani2017NewIT. Velocities of blood particles are encoded in the phase of the MRI signal while the anatomy is visualised from the signal’s magnitude Jones98. However, 4D flow has several limitations, such as low spatiotemporal resolution, long scan time, and low signal-to-noise (SNR) ratio JIANG2015185, which makes its clinical application to AR difficult. With a spatial resolution between 1.0 and 3.5mm Itatani2017NewIT, details on the narrowest part of the jet cannot be captured as its width is typically much smaller than 3mm Sallach2007 for mild AR. Therefore, spatial resolution is the biggest limitation in 4D flow MRI.

2.3 Deep Learning in 4D Flow MRI

Deep learning (DL) Deep-Learning has had a significant impact in many scientific sectors, and is highly relevant in the field of medical imaging Medical-Image-DL. Advances in super resolution (SR) image reconstruction SR-Image-Reconstruction to obtain high-resolution (HR) images from low-resolution (LR) observations are increasingly being adopted for MRI with a DL-based approach MRI-SR-GAN-DenseNet; Brain-MRI-SR. This approach is preferred as it not only has an advantage in spatial resolution quality over conventional SR techniques MRI-SR-CascadedDL, but also successfully denoises flow images Cerebrovascular-4DFlowMRI.

The combination of 4D flow MRI with deep learning has been explored in multiple ways to increase resolution and provide more accurate estimates of physical quantities 4DFlowNet; Physics-informed-4DFlowMRI; Ferdian-Cerebrov. However, there are several limitations involved in this approach. Primarily, these have been related to insufficient data 4DFlowNet; Ferdian-Cerebrov, which is due to the requirement of paired LR and HR MR images. This can be difficult to obtain as HR MRI takes long scanning times and is subject to motion artifacts 4DFlowNet; Ferdian-Cerebrov; Applications-DL-AA. As an alternative, CFD models have been used to simulate 4D flow MRI as ground truth HR images, which are then downsampled to LR images 4DFlowNet; Ferdian-Cerebrov. Other limitations include unstable and non-robust ANN architectures Physics-informed-4DFlowMRI, which describe the organisational structure of the ANN’s layers. The architectures plays a significant role in the performance of the DL algorithm, as well as ignoring phase/velocity aliasing error Cerebrovascular-4DFlowMRI; 4DFlowNet. The aliasing error here refers to aliasing from having a velocity encoding (VENC) velocity-encoding that is too low Aliasing-Encoding rather than other types of MRI spatial aliasing which have been explored previously and reduced DL-AntiAliasing-MRI-SR; SMORE-anti-aliasing; Brain-MRI-SR-old. Note that VENC is an MR parameter to adjust the maximum velocity corresponding to a 360 phase shift in the data.

2.4 Network Architecture

Recent development in object detection has proven significant in advancing ANN architecture Survey, and appears to be widely used in many medical imaging applications Medical-Image-DL.

2.4.1 Residual Blocks

Degradation is a common trend where very deep convolutional neural network (CNN)

CNN architecture leads to higher training error if trained for too long. This means that it is not an over-fitting problem since the training error is increasing. If layers can be constructed as identity mappings, it is expected that a deeper network should have a training error at most as large as its counterpart shallower network. This is because the shallower network is essentially a subset of its counterpart deeper network. Despite this, the implication of degradation is that CNNs have difficulty learning identity mappings.

To fix this problem of degradation, residual blocks are introduced in ResNet ResNet

. These contain a shortcut connection that skips neuron connections, seen in Figure

1, with the assumption that the output and input have the same dimensions. This allows for the desired mapping to be equal to the input, where the desired mapping is the sum of the residual block mapping and the input - if an identity mapping is desired, the residual block mapping will be equal to zero. It is hypothesised that the residual block mapping is easier to learn than the desired mapping. Furthermore, no additional parameters are needed as the identity mapping can point to the same part of memory as the input. This allows for deeper and more sophisticated CNNs, with the gradient being able to flow to earlier layers easier.

Figure 1: Residual block in ResNet, showing the use of identity mappings ResNet.

2.4.2 Dense Blocks

The weakness of shortcut connections in residual blocks is that they limit the learning capacity and representation power of the network. Dense blocks DenseNet are introduced within DenseNet in an attempt to solve this weakness, seen in Figure 2. The layers within these blocks receive feature maps from all preceding layers, known as concatenation. Transition layers are also required between dense blocks, as the input and output layers within them cannot change size.

Figure 2: Dense block with concatenated dense layers CSPNet.

Concatenation means that network layers can choose to use feature maps generated by preceding layers. This allows for diversification instead of correlation, which is an improvement on ResNet. Since earlier and later layers have weak and strong semantic information, respectively, layers are now built on features of low and high complexity. This substantially reduces the number of parameters as each layer only contains a constant number of new feature maps, defined as the growth rate DenseNet

. The vanishing gradient problem is alleviated too as the loss function can be easily propagated to earlier layers.

2.4.3 Cross Stage Partial Blocks

An important aspect of design is to optimise computational cost and memory allocation. The Cross Stage Partial Network (CSPNet) CSPNet aims to enhance learning capacity with reduction in computation, essentially optimising DenseNet in terms of speed and performance. The partial dense block is introduced, which partitions feature maps of the base input layer into two parts and then merges it through a cross-stage hierarchy at the transition layer. This is seen in Figure 3. For example, instead of taking the 30 feature maps at the base layer into the dense block, only a portion, say 10, are taken, with the rest taken straight to the output transition layer. This lets the gradient flow propagate through different paths, alleviating the vanishing gradient problem as well. The concatenation between the partitioned part not in the dense block and the transition layer allows for back propagation to jump over the dense block and flow easily back to earlier layers.

Figure 3: Cross stage partial block, showing the partitioning of the base input layer CSPNet.

There are several advantages of using partial dense blocks as part of the network architecture. It seems to strengthen learning ability as although lightweight versions of ResNet and DenseNet do not perform well, applying partial dense blocks to these networks reduce computational bottlenecks and improve performance CSPNet. Reducing computational bottlenecks refers to the skipping of a large number of base layer features, allowing for very minimal computational cost in the block. This is especially relevant in later layers, as computation can be evenly distributed at each layer and so effectively utilise arithmetic units. Finally, it also seems to reduce memory costs through cross-channel pooling, which is the reduction of channel size, compared to normal pooling that keeps the channel size (on top of reducing feature map size).

3 Methods

3.1 Data Generation

Modelling of the AV was done using Ansys 2021 R1, which has two main options for CFD simulations; CFX and Fluent. CFX was chosen as it is better suited for more simple, low Mach number flows, such as the flow through the AV. Modelling was an iterative process to determine the best parameters and options that would increase efficiency and accuracy.

3.1.1 Geometries

The most important aspects of the real-life AR behaviour are the velocity and pressure changes. With this in mind, the design of the problem geometry focused on simplicity over replicating a real-life AV, allowing for geometries to be created easily. The basic geometry is shaped like a cylinder with a constricted section part way along the length, resembling a Venturi tube. The constricted section represents the gap in the AV present for AR to occur. Figure 4 demonstrates a sketch of what the basic geometry looks like.

Figure 4: Basic sketch of the first 10 geometries.
Figure 5: Sketches of the angled (top) and offset (bottom) constricted sections, respectively.

In total, 20 different geometries were generated. The first set of 10 geometries had variations in inlet velocity, inlet radius, and constricted section radius, while the other set of 10 geometries had different shapes. These geometries were designed to capture maximum jet velocities between 2.5 and 5.0 ms. The second set of geometries were based on the third geometry, as this geometry was reasonably small with minimal computational time. This set had diagonal and off-centre constricted sections, as well combinations of both, to model eccentric jets. Figure 5 demonstrates these differences in shape. The parameters for all geometries have laid out in Table 1.

Basic Geometries111From left to right: the geometry number, the maximum inlet velocity in mms, the inlet radius in mm, and the constricted section radius in mm. Angled/Offset Geometries222From left to right: the geometry number, the angle between the constricted section and the direction normal to the inlet surface in degrees, the offset between the constricted section and the centre of the cylinder in mm, and the direction that the constricted section is angled.
No. R R No. Direction
1 300 5.00 1.00 11 20.0 0.00 Upward
2 150 5.00 1.00 12 40.0 0.00 Upward
3 500 5.00 1.50 13 0.0 1.50 -
4 100 5.00 0.75 14 0.0 3.00 -
5 100 5.00 0.60 15 20.0 1.50 Upward
6 450 8.00 2.00 16 40.0 1.50 Upward
7 450 6.00 2.00 17 20.0 3.00 Upward
8 100 8.00 1.00 18 40.0 3.00 Upward
9 150 10.0 2.00 19 30.0 2.25 Sideways
10 100 10.0 1.50 20 30.0 2.25 Downward
Table 1: Input parameters for all geometries

3.1.2 Boundary Conditions

The relevant boundary conditions are those relating to the inlet, outlet, and inner wall. The inlet boundary conditions take the most work to define as the velocity varies both in space and time. For blood flowing through the AV, the flow is expected to be fully developed, that is, the velocity is zero at the inner walls due to friction and at its highest in the centre. This can be achieved by extending the length that the blood needs to travel from the inlet to the constricted section, allowing the flow to fully develop before reaching the AV. However, extending the length means the geometry becomes larger, hence increasing the computational time. To compensate for this, the velocity profile at the inlet was defined with a parabolic shape and the upstream length was set to 20mm Madhavan2016. This means the flow will start out more developed than with a uniform profile, and become fully developed before reaching the AV. The velocity was also time-dependent and represented the diastole (where regurgitation occurs). An example profile is displayed in Figure 6, which shows how the velocity rapidly increases initially, then slowly decreases in magnitude before falling off quickly.

Figure 6: Velocity profile of aortic regurgitation at low intensity (left) and high intensity (right) over time Sallach2007. The large shape above the horizontal line represents the velocity through the AV during diastole.

The definition for the velocity profile was made using various expressions in Ansys. The following equation defines the velocity profile at the inlet in Cartesian coordinates:


where is the maximum velocity at the centre of the inlet and R is the inlet radius. The remaining boundary conditions were for the outlet and walls. The outlet was defined as an opening with zero pressure difference and the walls were defined as non-permeable with no-slip boundary conditions.

3.2 Data Preparation

To start, the raw CFD simulations and geometries, which each had 71 timeframes, were sampled onto a uniform Cartesian grid to be used as HR images. This was completed with linear interpolation, and took from 10 minutes to 3 hours per time frame depending on the size of the geometry and the CFD data. To separate the fluid and non-fluid regions for better data processing and result quantification, binary masks were generated using k-Nearest-Neighbours

knearestneighbours. The main difference between these synthetic HR images and MR images relates to the VENC and the amount of noise – HR images are noise-free whereas MR contain considerable noise. To account for this, the HR images were downsampled using the same method as in 4DFlowNet to simulate LR MR images, shown in Figure 7. This gives the paired LR and HR synthetic images used in network training.

Figure 7: Downsampling process used in 4DFlowNet involving six steps 4DFlowNet. This process is performed in 3D, but for visualisation purposes is shown in 2D.

3.3 Data Augmentation

To augment the data set, similar techniques were used as in 4DFlowNet for each time frame. VENC values were randomly chosen from a set of velocities between 0.3ms and 6.0ms, spaced by 0.3ms

, for each velocity component. Aliasing was mostly avoided by choosing a VENC larger than the peak velocity. However, since velocity jets cannot be estimated beforehand for AR (which may cause phase aliasing), a VENC lower than the maximum velocity was chosen with a 10% probability, randomly selecting between 0.3ms

or 0.6ms lower. Some of these produced patches are shown in Figure 8

. Constant intensity values between 60 and 240 were randomly chosen for the magnitude image, and noise levels were added depending on the signal-to-noise ratio, which were randomly chosen between 14 and 17 decibels.

Figure 8: Examples of aliased patches against their non-aliased counterparts. 2D slices of the velocity in the direction were taken from the 3D patch for visualisation purposes, along the width (left) and length (right) of the geometry. Scale is in metres per second.

Since there were a limited number of geometries, further augmentation came in patch generation. From each time frame, 10 patches of 12-voxel cubes from the LR image were selected randomly with a minimum fluid region of 20%. These patches act as random translations, so no extra translation steps were taken. On top of this, for each patch generated another randomly rotated version of the patch was also created. This resulted in 20 patches generated from each time frame and thus 1420 patches per geometry.

3.4 Training and Validation

To investigate the effect of additional geometries regarding SR image quality, a subset of the data consisting of patches from only five geometries was compared against a the entire data set consisting of patches from all geometries. The validation set in both cases was the same, consisting only of patches from a single geometry.

To investigate the effect of aliasing, duplicates of the two training and validation sets were generated, but with a 10% probability of having a VENC lower than the maximum velocity in any time frame. Networks trained using these aliased data sets were validated against the previously generated validation set without any aliasing, as well as a newly generated validation set with full aliasing, that is, with each time frame having a VENC lower than the maximum velocity.

3.5 Network Architecture

The simulated pairs of LR and HR synthetic images were used to train a similar deep residual network structure to the one in 4DFlowNet, shown in Figure 9. 4DFlowNet consists of several residual blocks surrounding a central upsampling layer, with the preceding blocks in the LR space pre-processing and acting as denoisers for the input while the following blocks in the HR space refine the output. In 4DFlowNet, LR patches of 16-voxel cubes were used as input and SR patches of 32-voxel cubes were generated as output, with an upsample factor of 2.

Figure 9: 4DFlowNet architecture, with the highlighted residual block in the middle 4DFlowNet.

A couple changes were made to the above 4DFlowNet architecture to provide higher resolution images with improved accuracy. Firstly, the upsample factor was increased to 4 and the sizes of the input and output patches were changed to 12-voxel and 48-voxel cubes, respectively. The smaller patches account for smaller vessel sizes in the cardiovascular space around the AV Ferdian-Cerebrov. Secondly, the dense and cross stage partial (CSP) blocks in DenseNet and CSPNet, respectively, were experimented with by using them in place of the 12 residuals blocks in the original 4DFlowNet architecture. The growth rate DenseNet

, defined as the number of feature maps in each convolutional layer, of the dense and CSP blocks was set to 16, a quarter of the number of channels in each convolutional layer from the original residual blocks. The adapted 4DFlowNet architecture with residuals blocks (4DFlowNet-Res), with dense blocks (4DFlowNet-Dense), and with CSP blocks (4DFlowNet-CSP) had 3.34, 2.55, and 2.08 million parameters, respectively. These modified networks were implemented with TensorFlow 2.0

Tensorflow20 and trained using an Adam optimiser kingma2017adam, with an initial learning rate of and decay rate of

after every 14 epochs. Batch sizes of 16 were used, with training completed 200 epochs.

3.6 Loss Function

The network was optimised by minimising the mean squared error (MSE) between the paired HR images and the SR images generated from the corresponding input LR ones. The voxel-wise loss was calculated as the mean sum of squared differences between each velocity component:


where is the total number of voxels in the geometry, is the predicted SR velocity, and is the actual HR velocity, for .

The MSE of fluid and non-fluid regions were calculated as separate terms due to the imbalance and irregularity of these regions within a specific patch. This gives the total loss to be:


where and are the voxel-wise loss for the fluid and non-fluid regions, respectively.

The original loss function in 4DFlowNet contained a weighted velocity gradient term to smoothen the gradient between neighbouring vectors

4DFlowNet. This was taken away from the above loss function as improvements were observed in near-wall velocity estimates with its removal Ferdian-Cerebrov.

3.7 Evaluation Metric

The relative speed error (RE), the relative difference between the SR velocity magnitude (speed) compared to the actual HR speed on the validation set, was used to measure network performance and save model checkpoints. This was only calculated in fluid regions to avoid zero division error, as well as adding a small number () to the denominator. Furthermore, since many speed values in the HR images were quite small, this could risk significantly over-penalising the model. Thus, an arctangent approach MAAPE was adopted, giving the following equation for relative speed error:


where is the total number of voxels in the fluid domain, and are the predicted SR and actual HR velocities, respectively, for all , and ‘arctan’ is the arctangent function, defined for all real values from negative infinity to infinity with for .

In addition to the RE, network performance was also evaluated using the root mean squared error (RMSE) and the structural similarity (SSIM) metric SSIM in all three Cartesian velocity components. These were compared against the baseline 4DFlowNet model that had been trained with an upsample factor of 2.

4 Results

Training was performed using a Tesla V100 GPU with 32GB memory with networks being trained for 200 epochs. Improvements in relative speed error (RE) plateaued around the 100 epoch mark for 4DFlowNet-Res while still improving for 4DFlowNet-CSP and 4DFlowNet-Dense up till the very last epoch. This can be seen in Figure 10. The time taken was dependent on the type of network; 4DFlowNet-Res, 4DFlowNet-CSP, and 4DFlowNet-Dense took approximately 163, 168, and 255 hours, respectively. Note that these times were for the networks trained using all geometries. For the networks trained using only five geometries, denoted as 4DFlowNet-Res5, 4DFlowNet-CSP5, and 4DFlowNet-Dense5, the times taken were approximately 38, 40, and 63 hours, respectively. There was no significant difference in training time between networks trained with and without a portion of aliased data, denoted by ‘-A’.

Figure 10: Relative error across all 200 epochs during training for each network. Networks trained with five geometries were not included as the general trend seen was the same.

Networks were tested on one complete geometry consisting of 71 timeframes with no aliasing, and the same complete geometry with full aliasing. These predictions were required to be patch-based since patches were used as the input and output for each network. The complete geometry was reconstructed by stitching together multiple SR velocity field patches, which was done with a stride of (

) in each Cartesian direction where n is the arbitrary patch size. To reduce the data to the centre of the patch, four voxels were stripped from each patch side.

4.1 Synthetic MR Images

Figure 11: Predictions on an LR patch in the constricted region from the synthetic 4D flow MRI phase image for different networks, with (top) and without (bottom) aliasing error. A 2D slice, along the width of the geometry, of the velocity magnitude is shown from the 3D patch for visualisation purposes. Scale is in metres per second.
Figure 12: Predictions on an LR patch, focused on the constricted section, from the synthetic 4D flow MRI phase image for different networks, with (top) and without (bottom) aliasing error. A 2D slice, along the length of the geometry, of the velocity magnitude is shown from the 3D patch for visualisation purposes. Scale is in metres per second.

SR images were analysed visually and quantitatively to better understand how each model was performing. Figures 11 and 12 are visual examples of the prediction for the different networks in the constricted section at the peak flow. These display the effectiveness of each network in reducing noise, with the predictions looking quite similar to the ground truth. Furthermore, the networks trained with a proportion of aliased data seem to be performing better than networks without, especially for data with aliasing error.

The values for each evaluation metric were collected and compared in Tables

2 and 3 to quantify the performance of every model. Again, these values were taken from the time frame with peak flow, with peak velocities of 2.186, 0.355, and 0.349 ms for velocity components , , and , respectively.

Table 2: Summary of prediction errors and evaluation metrics for different networks, when predicting on non-aliased data.
Network RMSE111RMSE ± s.d. for velocity component , in metres per second. RMSE111RMSE ± s.d. for velocity component , in metres per second. RMSE111RMSE ± s.d. for velocity component , in metres per second. SSIM222Includes all three Cartesian velocity components in the form (, , ). RE333Or equivalently, MAAPE.
4DFlowNet 0.05744 ± 0.05739 0.01979 ± 0.01977 0.01767 ± 0.01765 (0.832, 0.708, 0.705) 0.543
4DFlowNet-Res 0.03613 ± 0.03602 0.01379 ± 0.01378 0.01353 ± 0.01349 (0.923, 0.843, 0.820) 0.330
4DFlowNet-CSP 0.03535 ± 0.03458 0.01419 ± 0.01416 0.01539 ± 0.01527 (0.915, 0.728, 0.766) 0.316
4DFlowNet-Dense 0.03722 ± 0.03585 0.01382 ± 0.01377 0.01386 ± 0.0138 (0.908, 0.794, 0.789) 0.296
4DFlowNet-Res-A 0.03685 ± 0.03632 0.01363 ± 0.01355 0.01382 ± 0.01357 (0.924, 0.829, 0.824) 0.340
4DFlowNet-CSP-A 0.03663 ± 0.03662 0.01348 ± 0.01348 0.01245 ± 0.01242 (0.939, 0.812, 0.803) 0.330
4DFlowNet-Dense-A 0.0350 ± 0.03434 0.01488 ± 0.01482 0.01629 ± 0.01627 (0.912, 0.747, 0.705) 0.301
4DFlowNet-Res5 0.03752 ± 0.03751 0.01596 ± 0.0158 0.01634 ± 0.01616 (0.934, 0.796, 0.825) 0.344
4DFlowNet-CSP5 0.03987 ± 0.03973 0.01519 ± 0.01515 0.01452 ± 0.0142 (0.935, 0.812, 0.804) 0.338
4DFlowNet-Dense5 0.04069 ± 0.04064 0.0166 ± 0.01655 0.01499 ± 0.01471 (0.929, 0.764, 0.799) 0.352
4DFlowNet-Res5-A 0.04146 ± 0.04146 0.01961 ± 0.01958 0.01886 ± 0.01886 (0.924, 0.815, 0.817) 0.354
4DFlowNet-CSP5-A 0.04384 ± 0.04383 0.02194 ± 0.02191 0.02113 ± 0.02107 (0.929, 0.803, 0.811) 0.362
4DFlowNet-Dense5-A 0.04023 ± 0.04023 0.01938 ± 0.01932 0.01691 ± 0.01687 (0.928, 0.807, 0.784) 0.364
00footnotetext: Note: The best values for each metric are in bold. The best values for networks trained with only five geometries have also been coloured red.
Table 3: Summary of prediction errors and evaluation metrics for different networks, when predicting on aliased data.
Network RMSE111RMSE ± s.d. for velocity component , in metres per second. RMSE111RMSE ± s.d. for velocity component , in metres per second. RMSE111RMSE ± s.d. for velocity component , in metres per second. SSIM222Includes all three Cartesian velocity components in the form (, , ). RE333Or equivalently, MAAPE.
4DFlowNet 0.06419 ± 0.0641 0.02917 ± 0.02915 0.02716 ± 0.02715 (0.855, 0.724, 0.814) 0.512
4DFlowNet-Res 0.03611 ± 0.03606 0.01423 ± 0.01423 0.01381 ± 0.01377 (0.927, 0.858, 0.836) 0.334
4DFlowNet-CSP 0.03642 ± 0.03603 0.01608 ± 0.01606 0.01638 ± 0.01625 (0.917, 0.795, 0.778) 0.309
4DFlowNet-Dense 0.03555 ± 0.03435 0.01529 ± 0.01527 0.01594 ± 0.0159 (0.917, 0.812, 0.798) 0.309
4DFlowNet-Res-A 0.03547 ± 0.03495 0.01483 ± 0.01475 0.01395 ± 0.01381 (0.925, 0.845, 0.848) 0.342
4DFlowNet-CSP-A 0.0352 ± 0.0352 0.01578 ± 0.01577 0.01464 ± 0.01463 (0.942, 0.845, 0.849) 0.313
4DFlowNet-Dense-A 0.03402 ± 0.03434 0.01635 ± 0.01482 0.01639 ± 0.01627 (0.922, 0.785, 0.757) 0.284
4DFlowNet-Res5 0.04681 ± 0.04679 0.01764 ± 0.01757 0.01695 ± 0.01687 (0.865, 0.899, 0.582) 0.337
4DFlowNet-CSP5 0.05287 ± 0.05286 0.01743 ± 0.01743 0.01636 ± 0.01619 (0.912, 0.921, 0.509) 0.341
4DFlowNet-Dense5 0.05077 ± 0.05076 0.01777 ± 0.01767 0.01713 ± 0.01706 (0.902, 0.910, 0.587) 0.353
4DFlowNet-Res5-A 0.04366 ± 0.04365 0.01748 ± 0.01748 0.01584 ± 0.01584 (0.901, 0.902, 0.522) 0.344
4DFlowNet-CSP5-A 0.04663 ± 0.04653 0.01796 ± 0.01795 0.01637 ± 0.01636 (0.914, 0.925, 0.610) 0.350
4DFlowNet-Dense5-A 0.04735 ± 0.04729 0.01935 ± 0.01934 0.01619 ± 0.01618 (0.895, 0.902, 0.506) 0.338
00footnotetext: Note: The best values for each metric are in bold. The best values for networks trained with only five geometries have also been coloured red.

The metrics were plotted in Figure 13. Briefly, 4DFlowNet-Dense-A and 4DFlowNet-CSP-A seem to perform the best with the lowest RMSE error and largest SSIM in the principle flow direction (), respectively, on both aliased and non-aliased data. However, there does appear to be considerable variation in these results between different networks, depending heavily on the size of the data set and whether aliasing is present.

Figure 13: Comparison of networks for different metrics and results, varied by the data set size and whether aliasing was used. Values were normalised to be within 0 and 1, representing the worst and best values, respectively, for a specific metric relative to other networks.

Regarding the RMSE in the principle flow direction (RMSE), all networks perform better, on average, than the base 4DFlowNet with less variation in RMSE. For the smaller data set with five geometries, the residual-based (Res) versions seem to have the smallest RMSE. There is also noticeable improvement in error between networks trained on the smaller and larger datasets when predicting on non-aliased data. However, when predicting on aliased data, the improvement is significantly more apparent. The increase in dataset size improves the RMSE in the other two flow directions for all networks too, with the CSP versions somewhat better than other versions.

Regarding the SSIM in the principle flow direction (SSIM), all networks also perform better, on average, than the base 4DFlowNet. However, the SSIM seems to worsen for a larger dataset, in general, when predicting on non-aliased data. On the other hand, when predicting on aliased data, there does appear to be slight improvement in SSIM across all networks. For the SSIM in the other two flow directions, the values are considerably more varied and worse than in the principle direction. Finally, all networks have a substantially lower RE than the base, with the Dense versions having the lowest RE.

Figure 14: Regression and Bland-Altman plots for prediction on non-aliased (left) and aliased (right) data with 4DFlowNet-CSP-A. These plots are for each of the velocity components and magnitude (, , , and , respectively, from top to bottom) between SR and synthetic HR images.
Figure 15: Regression and Bland-Altman plots for prediction on non-aliased (left) and aliased (right) data with 4DFlowNet-CSP. These plots are for each of the velocity components and magnitude (, , , and , respectively, from top to bottom) between SR and synthetic HR images.

The regression plots in Figures 14 and 15 show that there is exceptional correlation between the SR and synthetic HR images. The regression slopes are very close to one, in the principle flow direction and velocity magnitude plots, for the two CSP networks. Moreover, offset values are also essentially zero in all examples shown. Training with all non-aliased images appears to have an effect when predicting the higher velocity values in aliased images, with these values being slightly underestimated, as seen on the right in Figure 15 and confirmed by the Bland-Altman plots too. These plots also indicate minimal bias as the deviations appear constant, uniform, and are all less than 0.08ms. Note that the general trends and comments mentioned here were present in all other networks as well, so further plots were not shown.

4.2 Actual 4D Flow MRI Data

Two sets of in-vivo 4D flow MRI data (CH34 and CH37) were acquired and used to show how the network would perform, seen in Figure 16. These were only LR images, with no HR images available to compare the predicted SR images against. However, the SR images produced were compared against the baseline 4DFlowNet to help better understand the predictions. The SR images seem to have effectively removed noise from the LR image, with the new 4DFlowNet-CSP networks appearing to remove slightly more noise throughout the LR image than the baseline 4DFlowNet. Other networks performed similarly, with predicted SR images almost identical to the ones shown.

The top set of data (CH34) was also segmented and visualised with ParaView paraview, shown in Figure 17. Image stitching appears to have been performed correctly, with velocities within the fastest section preserved well.

Figure 16: Predicted SR images on actual 4D flow MRI data. Two sets of 4D flow MRI data were used (top and bottom), with the original LR images on the left and the SR images in the other three columns, titled by the network used for prediction. These are 2D slices of the 3D image, showing the velocity magnitude taken at peak flow. Scale is in metres per second.
Figure 17: Predicted SR images visualised in ParaView. The columns show the velocity magnitude, its vector field, and the streamline reconstruction within the largest magnitde section. Scale is in metres per second.

5 Discussion

A noteworthy consideration when analysing the results is that the focus should be more towards metrics and values obtained in the principle flow direction . Since the peak velocities in the other flow directions, and , were almost 10 times smaller than that in the principle direction, networks would have difficulty differentiating between velocity and noise. This is evident in the and RMSE values, which were over half of their corresponding peak velocities. The and SSIM values were also much worse, potentially exhibiting strange behaviour in these low velocity fields too SSIM_limitations.

5.1 Synthetic 4D Flow MRI Data

The main limitation of previous work has been related to insufficient data and flow characteristics 4DFlowNet; Cerebrovascular-4DFlowMRI; Ferdian-Cerebrov. To understand the effect of additional flow characteristics in the dataset, a wider range of them were used in this study. There was definitely noticeable improvement in the RMSE across all networks tested, particularly when predicting on aliased data. This indicates that with more geometries and hence flow characteristics, the network seems to generalise better and is more robust against data it has not seen.

In terms of the synthetic LR and HR image pairs, the downsampling process and patch-based approach seemed to work effectively. The noisy LR patches allowed for the network to learn noise removal while also enabling greater generalisation to unknown flow characteristics or geometries 4DFlowNet. However, a small portion of these incorrect non-fluid velocity areas around the fluid domain appear to have been included in network training, seen in Figures 11 and 12. These were generated during the linear interpolation process when obtaining the HR images from the CFD data. This suggests that the binary mask for separating between the fluid and non-fluid regions, created using k-Nearest-Neighbours, needs improvement. An obvious consequence is that the network may learn incorrect flow characteristics and predict less accurately.

Incorporating aliased patches into the training data seemed to work effectively as well, improving the RMSE across networks trained with all geometries. This was done by choosing VENC values lower than the maximum velocity, within a particular time frame, with a 10% probability. Note that this probability was chosen arbitrarily. For the VENC, if it is set too high, visualization of the jet may not be obtained and be inaccurate, as well as having poorer SNR. On the other hand, if it is set too low, flow characteristics may be lost and a mosaic pattern will be shown MRI_venc. This means that for the time frames with aliasing, the velocities lower than the VENC would have been captured significantly more clearly, with only a small portion of high velocity characteristics being lost. For these time frames, the purpose would be to help the network better learn the flow characteristics in lower velocity fields, improving robustness and generalisability. However, this hypothesis has not been proven yet and more experiments on aliasing will be required to fully understand its effect on network performance. This would involve varying its probability of occurring as well as choosing VENC values considerably lower than the maximum velocity.

5.2 Network Architecture

A major drawback of residual blocks is that they suffer from limited learning ability DenseNet. This was seen in the results, in which the Res networks did not show much improvement even after adding many more geometries or introducing aliased images into the training data. Furthermore, the RE for Res networks seemed to plateau at around the 100th epoch, whereas the RE for the other networks seemed to continue improving, although at a decreasing rate, even up till the last epoch. These observations can be seen in Figure 10. Despite this, the Res networks still performed well, with RMSE, SSIM, and RE values similar to the other networks, as well as having the fastest training and prediction times. This suggests that, when data is insufficient or very limited, Res networks may work best.

With sufficient or abundant data, CSP or dense networks may be preferred. The learning ability in these types of networks have a significantly higher ceiling than Res networks DenseNet; CSPNet

, with the only difference between these two network structures being the training and prediction times. CSP networks had training and prediction times almost as fast as Res networks, whereas dense networks were almost 1.6 times slower. Although training times may not be a crucial problem, clinicians and patients alike may require and desire fast prediction times. This leads to CSP networks being preferred. Otherwise, if training and prediction times are not significant constraints, then dense networks may be the optimal choice. Furthermore, the growth rate for the dense and CSP networks was set quite low, at a quarter of the number of feature maps in each convolutional layer within the residual blocks. This limits the learning ability of these networks, as there are significantly less parameters, so testing larger values of this hyperparameter will be beneficial and likely improve network performance. Similarly, only a quarter of feature maps were taken from the base input layer within each partial dense block within the CSP networks, so larger values of this hyperparameter will likely improve performance as well.

Despite the network architecture seeming to work quite effectively, there are still improvements that could be made on top of modifying the residual blocks. Presently, the network is not taking full advantage of the temporal aspect in 4D flow MRI. Modifying it to incorporate characteristics of recurrent neural networks

RecurrentNN may help the network understand this temporal aspect better. This could be done by using predictions of the same patches from one or two time frames prior. Additionally, including physical properties of fluids may also help the network in learning flow characteristics, as seen in Physics-informed-4DFlowMRI; Ferdian-Cerebrov. However, due to the bulky nature of the velocity data, which were 3D cubes for each velocity component, these ideas were not considered further as it would have been too costly to process given the available resources.

5.3 Limitations

Although the number of geometries and flow patterns have significantly increased from previous studies, a wider range of characteristics can still be included. Additional data will diversify the data set further and improve the model’s robustness and generalisability even more. On top of this, more testing on real 4D flow MRI data will be required to validate model performance. Currently, this validation is only done by visually analysing model predictions to see if they look reasonable and sensible. The preferred approach would to be validate quantitatively with a pair of corresponding LR and HR 4D flow MRI images, calculating metrics such as RMSE and SSIM to properly understand model performance.

Lastly, several practical limitations can also be noted. Due to the limited GPU resources, training took approximately 60 days for all networks. With either more time or increased GPU resources, more networks could be trained by using different values for hyperparameters such as the growth rate and the probability that aliasing occurs. Moreover, networks could also be trained for more epochs, or until the error plateaus, to better gauge the learning ability of different networks. Finally, memory constraints were a factor too, as an upsample factor of 4 led to 64 times more usage in disk space. This would not be feasible for much larger upsample factors, so a different data representation may be required in future work.

6 Conclusion

In this study, it was shown how deep learning and computational fluid dynamics can be combined to effectively quantify hemodynamics for aortic regurgitation. 4DFlowNet was successfully adapted and modified to produce 4D flow MRI SR images with an upsample factor of 4. The results show that by adding more geometries and hence flow characteristics into the data set, the accuracy of 4DFlowNet predictions are improved. Moreover, the comparison of different network architecture suggested that the original residual network structure limits learning ability and can be further refined.