. Author manuscript; available in PMC: 2026 Jan 6.

Published in final edited form as: J Biomech Eng. 2024 Sep 1;146(9):091006. doi: 10.1115/1.4065165

Unsupervised denoising and super-resolution of vascular flow data by physics-informed machine learning

PMCID: PMC12768327 NIHMSID: NIHMS2126686 PMID: 38529728

Abstract

Purpose:

We present an unsupervised deep learning method to perform flow denoising and super-resolution without high resolution labels. We demonstrate the ability of this model to reconstruct 3D stenosis and aneurysm flows, with varying geometries, orientations, and boundary conditions.

Methods:

Ground truth data was generated using computational fluid dynamics, and then corrupted with multiplicative Gaussian noise. Autoencoders were used to compress the representations of the flow domain geometry and the (possibly noisy and low-resolution) flow field. These representations were used to condition a physics-informed neural network. A physics-based loss was implemented to train the model to recover lost information from the noisy input by transforming the flow to a solution of the Navier-Stokes equations.

Results:

Our experiments achieved mean squared errors in the true flow reconstruction of $O (1.0 e - 4)$ , and root mean squared residuals of $O (1.0 e - 2)$ for the momentum and continuity equations. Our method yielded correlation coefficients of 0.971 for the hidden pressure field and 0.82 for the derived wall shear stress field.

Conclusion:

By performing point-wise predictions of the flow, the model was able to robustly denoise and super-resolve the field to 20x the input resolution.

1. Introduction

A variety of fluid flows are, by nature, difficult to observe. Blood flow is no exception for several reasons including fluidopacity, being obscured by surrounding tissues and the presence of physiological factors. Computational fluid dynamics (CFD) has proven to be successful in a wide range of cardiovascular hemodynamics modeling applications, with a growing base of software to perform such modeling [1, 2]. However, one of the biggest challenges with CFD, particularly in non-idealized geometries where complicated flow structures exist, is uncertainty of boundary conditions. While boundary conditions can, in theory, be specified directly from measurements, more often assumptions and intuition play a large role. This uncertainty can be attenuated in applications where rough measurement of the flow in the domain interior is also available. Namely, in applications where magnetic resonance velocimetry (MRV) can be performed [3], there is the potential to use such data to calibrate boundary conditions (or other uncertain inputs) to improve modeling performance [4,5]. Nonetheless, the best way to merge data- and theory-based models is largely an open research question, particularly given recent advances in machine learning.

As an alternative to traditional CFD, machine learning techniques, and specifically deep learning, have emerged as a promising framework to tackle a variety of fluid modeling challenges [6]. For instance, various machine learning algorithms have been explored to enhance the spatial and temporal resolutions of fluid flows. Some of these rely on ridge-regression [7] or Lasso regression [8] applied to a known high resolution dataset to extract key features of such flows. Others rely on combining CFD simulations with data using adjoint-based methods [9], or with low resolution simulations using a differentiable solver and neural networks [10]. The main drawbacks of these methods include the uncertainty of boundary conditions or cost of running a CFD simulation, particularly at higher resolution. Deep learning physics-informed methods offer an alternative framework to merge data- and theory-based models. Several recent reviews highlight the potential for physics-informed neural networks in the modeling of fluid mechanics, dynamical systems and cardiovascular biomechanics [11–13].

Two main approaches exist when training a machine learning model on a theory-based (i.e., partial differential equation, PDE) model. A first approach relies on numerical discretization to impose PDE constraints as a loss function on the output of a convolutional neural network (CNN). For example, Subramaniam et al. [14] successfully applied this method for turbulence enrichment to transform a low-resolution turbulent field to a high-resolution counterpart whilst imposing the continuity and the pressure Poisson equations. Despite its computation efficiency, the discretized stencils often limit the geometrical complexity of the applications that can be readily considered (e.g., precluding cardiovascular domains). Recently, Gao et al. [15] addressed this problem by solving elliptic equations to map irregular grids to regular grids, though the geometric complexity of the domain remains constrained by the method and the approach is limited to fixed super-resolution scales.

A second approach relies on the differentiable nature of neural networks to enforce the PDE. Raissi et al. [16] set the inputs of a neural network as the Cartesian coordinates $x_{i}$ of a point in the flow domain, and the output as being the velocity components $u_{i}$ and the pressure field $p$ at that point. By doing so, the residual of the PDE can be directly computed by taking the differentials of $u_{i}$ and $p$ with respect to $x_{i}$ . The residuals are then used to train the model to satisfy the Navier-Stokes equations (NSEs) continuously over the domain. This continuous evaluation opens the door to more complex geometries, which could be particularly useful for image-based vascular geometries. Raissi et al. [16] coined the term Physics-Informed Neural Networks (PINN), to describe such networks. For instance, Sun and Wang [5] introduced an idealized stenosis parameter in the physics-informed neural network input to train the model over various 2D ideal stenosis geometries with varying viscosity. This method does not leverage any image data but assumes exact knowledge of boundary conditions and is limited to parameterized 2D stenosis geometries. In a similar vein, Oldenburg et al. [17] combined a variational autoencoder [18] with a physics-informed neural network to solve for a distribution of more complex 2D stenosis geometries. This model assumes a constant Reynolds number ( $R e$ ). Some authors have also been successful in reconstructing fluid flows in a variety of complex geometries without such constraints, by using CFD simulations as the training data. Maleki et al. [19] introduced a geometry autoencoder that was used by Ranade et al. [20] to develop a fluid flow neural network solver. Rutkowski et al. [21] trained a CNN to enhance cerebrovascular 4D flow MRI. A drawback is that these methods assume the presence of exact ground truth data and do not ensure physical viability of the model’s output.

There is substantial potential to leverage physics-informed learning to improve MRV data. Namely, the temporal and spatial resolutions of MRV are generally poor, measurements are noisy, and no direct measurement of blood pressure is available [22,23]. While CFD simulations based on coarse MRV data have been considered, this approach has been shown to poorly estimate derived quantities such as wall shear stress, a critical hemodynamic parameter [24]. Alternatively, physics-informed neural networks have successfully been used to super-resolve noisy flow inputs in both the discretized [15,25] and differentiable [26,27] training frameworks. Fathi et al. [27] used 4D flow visualisation images, transformed into the complex plane, to study the use of physics-informed neural networks to super-resolve and denoise the input without segmentation labels. The authors exhibited the performance of their method on example flows such as an aneurysm flow, and revealed that the input can be efficiently super-resolved and denoised. In addition, Kissas et al. [28] recovered the absolute pressure wave in a bifurcating arterial network from noisy clinical MRV data using one-dimensional models of pulsatile flow to constrain a physics-informed machine learning model. In both articles, the authors demonstrate the performance of their method on an example cardiovascular flow, and reveal that a single PDE solution can be learnt efficiently. However, these models were trained for a specific geometry, which limits their practical utility.

Herein, we are motivated by the goal of denoising and super-resolution of MRV. However, there is typically no ground truth available with MRV data. Thus, for the purposes of testing our proposed framework, we herein apply our method to synthetic flow data derived from CFD whereby a range of ground truth can be generated for evaluating accuracy. And since, broadly, the most important variables governing local cardiovascular hemodynamics are geometry and inflow/outflow boundary conditions, we develop a method capable of predicting flow and pressure fields when vascular geometry and/or boundary conditions are varied. The contributions of our work herein are threefold: (1) We develop a methodology to train a single physics-informed machine learning model that can perform denoising of low resolution and noisy 3D flow images over a spectrum of geometries and $R e$ . (2) We demonstrate the ability of our model to recover the hidden pressure gradient field and reasonable estimation of derived hemodynamic parameters such as the wall shear stress that are highly sensitive to resolution and noise. (3) We demonstrate the ability of our model to be queried at arbitrary resolution, effectively performing super-resolution of the flow without high-resolution labels.

2. Methods

A super-resolution model broadly seeks to learn a mapping from a noisy and low resolution field ${\tilde{u}}_{i}$ to a high resolution field $u_{i}$ . Here, $u_{i}$ represent components of the fluid velocity in a domain of interest $Ω$ , where $i = 1,2, 3$ signify the 3 coordinate directions. Ideally $u_{i}$ represents the “true” flow field, i.e., the unique hypothetical flow field resulting from the exact domain geometry, boundary conditions and fluid properties.

We assume that the flow field $u_{i}$ should satisfy the NSEs for an incompressible, Newtonian fluid. Herein, we will also assume steady flow with no body forces. The resulting NSEs can be non-dimensionalized as to be parameterized solely by the Reynolds number $R e$ . Namely, define $U$ , and $L$ , as the respective velocity and length scales of the flow. We can then define the following non-dimensional variables:

\{\begin{cases} u_{i}^{'} = \frac{u_{i}}{U} \\ x_{i}^{'} = \frac{x_{i}}{L} \\ p^{'} = \frac{p}{ρ U^{2}} \end{cases}

(1)

where $x_{i}$ denotes coordinate position ( $i = 1,2, 3$ ) and $p$ is the fluid pressure. The momentum equation is then given by

u_{j}^{'} \frac{\partial u_{i}^{'}}{\partial x_{j}^{'}} = - \frac{\partial p^{'}}{\partial x_{i}^{'}} + \frac{1}{R e} \frac{\partial^{2} u_{i}^{'}}{\partial x_{k}^{'} \partial x_{k}^{'}}

(2)

and the continuity equation is

\frac{\partial u_{i}^{'}}{\partial x_{i}^{'}} = 0

(3)

where $R e = \frac{U L}{v}$ , with $ν$ representing the fluid kinematic viscosity assumed constant. To ease notation burden, we remove the prime values in subsequent expressions of these variables.

2.1. Dataset

In this section we present the data generation and data processing used for our experiments. Consistent with the nominal properties of blood, we set $μ = 0.04 g r . (c m . s)^{- 1}, ρ = 1.06 g r . {c m}^{- 3}$ for all generated flows. Note that we express the spatial coordinates as $x_{i}$ and $x, y, z$ interchangeably. The former is used to write the equations compactly using index notation, while the latter is favored for figures and analysis.

2.1.1. Flow data

We used Simvascular [1] to perform CFD simulations of idealized stenotic and aneurysmal geometries. These geometries can be parameterized with a radius $r = r (y)$ described by the following equation:

r = R \cdot [1 + \frac{\pm α}{2} \cdot (1 + c o s (2 π \cdot \frac{y - y_{c}}{L}))]

(4)

where $α$ is the obstruction fraction of the native vessel radius $R, y_{c}$ is the axial position of the center of the stenosis/aneurysm, and $L$ its length. For a stenosis we use $- α$ to reduce the cross-sectional area, while we use $+ α$ to increase the cross-sectional area for an aneurysm. We set a total length of 20 cm and base radius $R$ of 2 or 3 cm. For the stenosis flows, we varied the inflow in the range [40, 140] cm³.s⁻¹, and the radius obstruction in the range [0.2, 0.5], corresponding to reductions of area of 36 to 75%. For aneurysm flows, we varied the inflow in the range [35, 75] cm³.s⁻¹, and the obstruction in the range [0.5, 1.5]. The inflow is smaller for the aneurysm models in order to generate flows that exhibit greater recirculation patterns. Overall, the flows ranged with $R e \in [460,2000]$ , using each flow’s maximum velocity and diameter taken at the middle of the stenosis/aneurysm. For data generation, the inflow boundary condition was prescribed a parabolic profile, and a zero pressure condition was prescribed at the outlet (since the vessel was modeled as inelastic). Based on these parameters, we generated 315 different stenosis flows, and 315 aneurysm flows, with minimum mesh sizes of 75,000 nodes.

2.1.2. Generating low resolution image data

The CFD simulations result in flow field data output on an unstructured mesh. To be more representative of measured flow data, we downsample the unstructured mesh data onto a 3D image (structured) grid that forms a bounding box over the domain. The image grid is 64³ and downsampling is achieved by tri-linear interpolation using the Python SciPy griddata function. We preserve the geometric aspect ratio, which, because the domain is narrow in two dimensions, results in relatively sparse 3D data where only ~25,000 voxels contain (nonzero) flow information. This sparsity is typical in measured data. To evaluate the error induced by this downsampling, we interpolated the CFD positions from the 64³ bounding box and compare the CFD and reconstructed fields. The interpolation was computed to introduce a mean squared error (MSE) of 3.1e − 2 between the two fields.

2.1.3. Transformations

We perform translations, rotations, and inversions (reversal) of the downsampled, voxelized flows to evaluate the ability of the model to perform independently of the vessel orientation, position, and flow direction. We introduce variations in these flows by random: translations in the range ±15 pixels (on a volume of size 64³), rotations in the range [−90°, 90°], and inversion of the flow direction. For each experiment, we use the hold-one out validation method, where we select and transform: 80% of the stenosis and aneurysm flows for training, the other 20% of the flows for validation.

2.1.4. Signed Distance Fields

The geometry is segmented by labeling image voxels as either inside or outside the flow domain. For simplicity, the voxels were labeled as outside if they are zero-valued in the three velocity components, or otherwise labeled inside. Voxels composing the boundaries of the flow domain were, in the case of the vessel wall obtained by finding the edge between voxels labeled outside and those labeled inside, and in the case of inlets and outlets determined by the intersection of the fluid domain with the bounding box. The geometry information is further encoded by computing a signed distance field, whereby each voxel first gets a value equal to its distance to the closest boundary. Namely, this value is computed as the Euclidean distance between the center position of the voxel of interest and the center position of the closest boundary voxel. If the voxel is within the background, the distance is made negative, and otherwise the distance remains unaltered. This step enables richer information than a segmentation map, as points far and close from the boundary are encoded differently. An example of a signed distance field is shown in Figure 1.

Fig. 1. — Overview of the model architecture. The noisy flow volume is passed through the flow CNN autoencoder from which the latent representation $γ^{f l o w}$ is extracted. The signed distance field is built from the flow volume and passed through the geometry CNN autoencoder to generate the geometry latent representation $γ^{geo}$ . Both $γ^{flow}$ and $γ^{geo}$ are combined, with the space coordinates $x, y, z$ of a point of interest to predict the velocity vector and pressure scalar at that point. (Here the $s$ and $a$ subscripts described in the text are removed.)

2.1.5. Normalization

We normalize the input variables through non-dimensionalisation, using the velocity and length scales of the flow. The velocity scale is set to the maximal magnitude of all the generated flows. By scaling all flows with a shared velocity scale, we ensure that our model always solves the same NSE, parameterized by a constant $R e$ . This makes the task less complex for the model. Indeed, if each sample were to have its own normalization scheme, then the model would have to learn a solution for each NSE. Hence, the velocity values of the normalized flows range in [−1, 1]. Similarly, we normalize the signed distance fields between [0, 1]. Hence, for points outside the domain, the further the point is from the flow domain the closer its value is to 0. Doing this reduces the impact of non-fluid regions on the learning process, and ensures that the network focuses more on reconstructing the flow and boundary regions.

2.1.6. Noise

Multiplicative noise can be generated during coherent imaging processes acquisition [29], such as in MRV. Hence, in order to simulate noisy flow data from an MRV acquisition, we add multiplicative Gaussian noise in a similar fashion to [15] and [30]:

{\tilde{u}}_{i} = u_{i} * (\vec{I} + 𝒩 (\vec{0}, \vec{I}))

(5)

$𝒩$ represents the normal distribution of mean $\vec{0}$ and variance $\vec{I}$ , given by identity so that noise is added independently for each voxel and velocity component. Even though the unit variance may seem excessive, it is explicitly selected to test the ability of the model to remove a significant level of noise. After the noise addition, the flow values ranged in [−4, 4]. To quantify the amount of noise introduced in the 3D flow images, we used the signal to noise ratio defined as:

S N R = 10 {l o g}_{10} \frac{σ_{u}^{2}}{σ^{2}}

(6)

where $σ_{u}$ is the standard deviation of the flow field.

2.2. Machine Learning Framework

2.2.1. Overview

Our model comprises three main components described below and depicted in Figure 1.

The geometry autoencoder takes in the signed distance field of each sample and reduces it to a latent space representation of 60 features. For a given input sample $s$ , we label the vector representation of the geometry as $γ_{s}^{g e o}$ .
The flow autoencoder takes in the velocity field as input, and reduces each sample into a latent space representation of 60 features. For a given input sample $s$ , we label the latent vector representation of the flow as $γ_{s}^{flow}$ .
The PINN takes as inputs the $x_{a}, y_{a}, z_{a}$ spatial coordinates of a point $a$ from a sample $s$ , alongside the latent space representations of both the geometry $γ_{s}^{geo}$ and the flow $γ_{s}^{flow}$ . It then outputs a prediction for the flow values ${\hat{u}}_{i}^{s} (a)$ and pressure ${\hat{p}}^{s} (a)$ at that point.

Following a standard approach in machine learning for image analysis, we processed the image data in 64³ patches. Each patch represents a region of 20 cm³ in physical space. These patches can then be stitched together to cover the entire domain as needed; however this recombination process is secondary to the process described herein.

The flow and geometry autoencoders are trained individually first. The weights of these models are then frozen when training the PINN, and the decoder is discarded. We describe each of these individual components below.

2.2.2. Geometry Autoencoder

The geometry autoencoder conditions the PINN on the flow geometry. This enables a single model to be used for different geometries (e.g., vessel shape, orientation, and position).

We train an autoencoder to map the 3D signed distance field of each sample to a latent space. The rationale is that the latent space vector of the geometry will encode the necessary representational power to condition the PINN on the sample’s specific geometry. Since the input resembles an image, we use a CNN to train the autoencoder. This is consistent with approaches and models from other works [19, 20, 31, 32].

We perform a wide search on hyper-parameters and model architectures, resulting in parameters shown in Table 1. The interpolation scale factor corresponds to the resizing of the feature maps. Between each encoder layer, the maps are scaled down by a factor of 2 in each direction using trilinear interpolation. Conversely, the maps are then scaled up by a factor of 2 between each decoder layer, using using trilinear interpolation, to recover the original input size. The model is fully convolutional, with no linear layers in the latent space. This provided better results, and stems from the fact that throughout every layer in the model we benefit from the local features of the kernel.

Table 1.

Hyper-parameters of the geometry autoencoder.

Number of channels per layer	Kernel Size	Padding	Stride	Interpolation scale factor	Activation

Enc: 1, 3, 6, 12, 24, 48, 96, 60 Dec: 96, 48, 24, 12, 6, 1	3×3×3	1	1	2	ReLU

Open in a new tab

We train the model with the mean average error (MAE) loss, defined as:

ℒ_{geo} = \frac{1}{N^{3} \cdot B} \sum_{s = 1}^{B} \sum_{a \in Ω_{i m}^{l}} |d_{a}^{s} - {\hat{d}}_{a}^{s}|

(7)

where $N$ is the number voxels in each direction $x, y, z$ of the input sample, $B$ is the batch size, $d_{a}^{s}$ is the distance value (normalized signed Euclidean distance to the boundary) at point $a$ in sample $s$ and ${\hat{d}}_{a}^{s}$ is the prediction. $Ω_{i m}^{l}$ is the whole domain of the low resolution input image, which includes background, boundary and flow points. This loss provided better reconstruction results than the mean squared error (MSE) for this task. We use a learning rate of 10⁻⁴ and add a scheduler with exponential decay, which stabilises the training as the loss decreases. Every 20 epochs, the learning rate is multiplied by 0.99.

2.2.3. Flow Autoencoder

We introduce the flow samples in the model through the flow autoencoder. The role of this model is to compress the flow information into a feature vector $γ_{s}^{flow}$ to condition the PINN on the image-based flow data. Similar to the geometry encoder, we use the same CNN structure for the flow autoencoder. The only difference lies in the input and output dimensions. For the flow data, the input and output is the velocity vector field – thus having dimension of $[3, N, N, N]$ . The hyper-parameters of the CNN are supplied in Table 2.

Table 2.

Hyper-parameters of the flow autoencoder.

Number of channels per layer	Kernel Size	Padding	Stride	Interpolation scale factor	Activation

Enc: 3, 3, 6, 12, 24, 48, 96, 60 Dec: 96, 48, 24, 12, 6, 3	3×3×3	1	1	2	ReLU

Open in a new tab

We train the flow autoencoder with the mean squared error loss, defined as:

f l o w = \frac{1}{3 \cdot N^{3} \cdot B} \sum_{s = 1}^{B} \sum_{a \in Ω_{i m}^{l}} \sum_{i = 1}^{3} {({\tilde{u}}_{i}^{s} (a) - {\hat{\tilde{u}}}_{i}^{s} (a))}^{2} .

(8)

where ${\tilde{u}}_{i}^{s} (a)$ is the low resolution velocity vector in direction $x_{i}$ at point $a$ of sample $s$ , and ${\hat{\tilde{u}}}_{i}^{s} (a)$ is the prediction. Note, we explored training the model with the mean average error and cosine similarity losses. The cosine similarity was expected to provide better results, as in [33], since it considers the directional difference between two vectors. Indeed, the authors developed a direction-sensitive loss function to learn flow super-resolution that outperformed the mean squared and average errors. Nonetheless, we observed that the MSE provided the best results. We used the same learning rate and scheduler as for the geometry autoencoder.

2.2.4. Physics-Informed Neural Network

The PINN combines the latent representations of the flow and geometry autoencoders to provide pointwise solutions to the NSE for each sample. The latent representations condition the model on both the geometry and flow image data. For each point $a$ of a given input sample $s$ , we input the geometry and flow feature vectors $γ_{s}^{geo}$ and $γ_{s}^{flow}$ respectively. We define the concatenation:

γ_{s} = [γ_{s}^{geo}, γ_{s}^{flow}]

which acts as the conditioning of the PDE and the model. We then concatenate $γ_{s}$ with the normalized $x, y, z$ coordinates to generate the input vector for our PINN. For each point $a$ from a sample $s$ , the input vector hence reads:

c_{a}^{s} = (x_{a}, y_{a}, z_{a}, γ_{s})

of length 123, with 3 features for the position, 60 for the geometry and 60 for the noisy, low resolution image. The input data is a vector, hence we use a multi-layer perceptron model to build the PINN. We assume that the ordering and locality of the features have no particular information, hence a CNN is not deemed appropriate. The output of the PINN is a vector of length 4 composed of the velocity vector and the scalar pressure: $u_{1}, u_{2}, u_{3}, p$ . The hyper-parameters of the PINN are summarized in Figure 1. We use the hyperbolic tangent (tanh) as the activation function (instead of the ReLU) since the activation must be twice differentiable to compute the Laplacian of the velocity field in the NSE.

The role of the PINN is to predict the velocity vector and the scalar pressure at a given point $a \in Ω$ , for a given sample $s$ . The optimization task can be expressed as:

\tilde{Θ} = \underset{Θ}{a r g m i n} \sum_{s = 1}^{n} \sum_{a \in Ω} ‖P I N N (c_{a}^{s}) - u_{i}^{s} (a)‖

(9)

where $u_{i}^{s} (a)$ is the $i$ th component of the ground truth velocity at point $a$ of sample $s, Θ$ are the trainable parameters of the PINN, $n$ corresponds to the number of samples in the training set, and $‖ \cdot ‖$ is a norm operator. Notice that the domain here is the continuous domain $Ω$ , and not that of the low resolution data $Ω^{l}$ .

We set-up our training process as a multi-objective optimization task. We want to achieve two goals: (1) reconstruct the flow input to perform data assimilation, and (2) ensure the output is a solution to the NSE, thereby learning the pressure gradient field associated with the input flow. While obtaining pressure through the pressure Poisson equation would be possible, we expect that the presence of noise and coarse data would lead to large errors in spatial gradients and hence pressure calculations.

Reconstruction Loss

The reconstruction loss used for the PINN is the mean squared error, similar to Eq. (8) for the flow autoencoder. The input images correspond to the velocity field only. Therefore, the reconstruction loss only applies to the velocity terms. On the other hand, pressure is only constrained by the PDE loss and the associated gradient term. Moreover, note that we evaluate the loss at any point $a \in Ω$ , to ensure that our model can be queried continuously over the domain. Therefore, we must compute approximations of the velocity field ${\tilde{u}}_{i}$ at points $a \notin Ω^{l}$ . This interpolation process does not add new information per se; rather, it enables evaluation anywhere within $Ω$ . If training was limited to the grid locations of the low resolution input image, then the model could overfit to these locations.

To perform the interpolation, we use the torch.nn.grid sample function from Pytorch. This function performs a trilinear interpolation on the noisy and low resolution ${\tilde{u}}_{i}$ field to provide a reconstruction target for the points $a \notin Ω^{l}$ . We label this interpolation $I_{a}$ . The interpolation relies on the neighboring points of $a$ in $Ω^{l}$ . Over a batch of $B$ of samples, we select $n_{p}$ points per sample to train the model. We label the set of selected points as $ω$ . Note that the point locations vary for each sample. The reconstruction loss becomes:

ℒ_{recon} = \frac{1}{n_{p} \cdot B} \sum_{s = 1}^{B} \sum_{a \in ω} \sum_{i = 1}^{3} {[{\tilde{u}}_{i}^{s} (I_{a}) - {\hat{u}}_{i}^{s} (a)]}^{2}

(10)

Navier-Stokes Equation Loss

We apply the NSE loss to flow points $a_{p}$ such that $a_{p} \notin \partial Ω$ . We label this set of points as $ω_{p}$ . Based on the output vector ${\hat{u}}_{i}^{s} (a)$ and ${\hat{p}}^{s} (a)$ , the residuals at a specific point $a$ are obtained as:

ℛ_{cont}^{s} (a) = \frac{\partial {\hat{u}}_{i}^{s} (a)}{\partial x_{i}}

(11)

\begin{array}{r} ℛ_{mtm, i}^{s} (a) & = \frac{1}{Re} \frac{\partial^{2} {\hat{u}}_{i}^{s} (a)}{\partial x_{k} \partial x_{k}} - \frac{\partial {\hat{p}}^{s} (a)}{\partial x_{i}} - {\hat{u}}_{j}^{s} (a) \frac{\partial {\hat{u}}_{i}^{s} (a)}{\partial x_{j}} \end{array}

(12)

$ℛ_{cont}^{s} (a)$ represents the residual of the continuity equation, whilst $ℛ_{mtm, i}^{s} (a)$ that of the momentum equations. Hence, the loss is computed as the mean square error of the residuals:

ℒ_{p} = \frac{1}{n_{p} \cdot B} \sum_{s = 1}^{B} \sum_{a \in ω} [ℛ_{cont}^{s} (a)^{2} + \sum_{i = 1}^{3} (ℛ_{mtm, i}^{s} (a)^{2})]

(13)

Note that the boundary condition loss is included in $ℒ_{recon}$ . Moreover, $ℒ_{p}$ can be computed at any point in $Ω$ to ensure adherence to the physical laws expressed in the PDEs. We combine the two losses with a weight parameter $λ$ :

ℒ_{total} = ℒ_{recon} + λ ℒ_{p} .

(14)

The parameter $λ$ enables us to control the importance of the PDE loss. A large $λ$ drives the model to favor the PDE loss, whilst a smaller $λ$ favors the exact reconstruction of the (potentially noisy or coarse) input. For points on the vessel wall boundary, $a \in \partial Ω$ we have $u_{i} (a) = 0$ . We label this set of points $ω_{b}$ , and automatically set $γ = 0$ to compute only the reconstruction loss. Indeed, the no-slip vessel walls points must enforce the boundary conditions. In practice, we set $B = 48$ , and enforce the NSEs at $|ω_{p}| = 100$ flow points and the no-slip condition at $|ω_{b}| = 10$ boundary points (thus $| ω | = = 110$ in Eq. (13)).

2.2.5. Evaluation Metrics

We use five metrics to evaluate the performance of our model on the validation data.

Noiseless $ℒ_{recon}$ : the reconstruction loss of the model on the true flow; which it has never seen during training when $σ \neq 0$ .
Peak signal-to-noise ratio (PSNR): a quality measurement between the noiseless and reconstructed flow. The larger the value the better the model performance. The PSNR is defined as:
$P S N R = 10 {l o g}_{10} (\frac{m a x (u_{i})}{ℒ_{recon}})$ (15)
where $ℒ_{recon}$ must be evaluated with respect to the noiseless data, and $u_{i}$ is the ground truth flow field.
$ℒ_{p}$ : the physics loss value evaluated on flow points at the resolution of the low dimensional grid.
RMSE: the root of the mean of the squared error; which we use to evaluate the performance of the model on satisfying the continuity and momentum equations.
Pearson correlation coefficient $ρ$ to evaluate the recovery of physical parameters from our model. The absolute pressure value is non-unique up to a constant shift, thus the correlation coefficient enables us to compare the predicted and the true (CFD) pressure fields indifferent to this shift. We perform the same comparison for the wall shear stress.

3. Results

In this section we present the results on flow denoising and super-resolution for the different types of flows discussed in section 2.1. We explore the model’s denoising ability depending on the PDE loss weight $λ$ and the noise level $σ$ . We present results for the recovery of the pressure field for these different models, the wall shear stress and super-resolution performance for the best overall model. For benchmarking, we compare our three stages approach to a Baseline CNN method that consists of taking the output of the flow autoencoder as the denoised flow (as motivated below).

3.1. Denoising

We explore the impact of two specific parameters on our model performance: (1) the noise level σ and (2) the NSE loss weight $λ$ . It is important to note that for the settings where $σ > 0$ , the model never has access to the true flow. The input to the PINN, and the reconstruction targets are the noisy flow. Changing $σ$ enables us to quantify the denoising ability of our model.

We trained a geometry autoencoder for these flows according to section 2.2.2 and achieved a test loss of $ℒ_{geo} = 3.87 e - 3$ . The geometry was reasonably encoded using 60 features for conditioning the PINN. We trained a flow autoencoder for $σ = 0.0$ and $σ = 1.0$ , reaching MSE of 2.32e-5 and 3.13e-4 respectively, evaluated on the entire input images $Ω_{i m}^{l}$ . The SNR of the flows with $σ = 1.0$ was −4.51, indicating the large amount of corruption to the input. Since the MSE increases as the noise level increases we observed that the flow autoencoder is unable to learn the high frequency information that the noise adds to the flow. In that sense, the autoencoder already contributes to the denoising process, as mentioned in [34] and forms our Baseline CNN comparison method. Note that these MSE values should not be compared with $Noiseless ℒ_{recon}$ since their sparsity drives the loss closer to 0. Results for setting $λ$ to [0.1, 0.5, 1.0] are presented in Table 3.

Table 3.

Model reconstruction performance for stenosis and aneurysm flows for varying λ and σ. Each model was trained for 6000 epochs on a A6000 Nvidia GPU. The results were evaluated on the validation dataset.

	Baseline CNN		PINN
σ	${Noiseless ℒ}_{recon}$	PSNR	λ	${Noiseless ℒ}_{recon}$	PSNR

0.0	1.49e-4	38.3	0.0	1.03e-4	39.9
			0.1	5.86e-5	42.3
			0.5	7.61e-5	41.2
			1.0	9.70e-5	40.1

1.0	3.91e-4	34.1	0.0	2.52e-4	36.0
			0.1	1.02e-4	39.9
			0.5	1.48e-4	38.3
			1.0	1.82e-4	37.4

Open in a new tab

The data in Table 3 demonstrate that for an NSE loss weight $λ > 0.0$ , increasing $λ$ effectively increases the noiseless reconstruction loss of the PINN, decreasing the PSNR. This is expected since the impact of $ℒ_{p}$ on the total loss increases, which leads the model to favor the NSE over the reconstruction of the input data. For instance, with $σ = 1.0, λ = 0.1$ gives a PSNR of 39.9 while $λ = 1.0$ gives a PSNR of 37.4. Despite these observations, the PINN outperforms the Baseline CNN independently of $λ$ , for which the PSNR is 34.1, whilst being robust to noise. Indeed, the results are of similar order of magnitudes (PSNR~ 40). This supports the overall approach, and the ability of the knowledge of the NSE to denoise the input flow. We demonstrate such denoising abilities in Figure 2 for two representative validation stenosis and aneurysm flows. The comparison between the Baseline CNN and the $λ = 0.0$ case highlights that the PINN architecture naturally provides improvements on the task. This could be due to the explicit incorporation of the geometry signed distance fields, which help the PINN to learn shared features across different examples. Additionally, comparing the $λ = 0.0$ case with $λ > 0.0$ cases reveals that the consideration of NSE loss reduces the output noise, by a factor ranging from 1.4 to 2.5, reinforcing the impact of combining the NSE constraint and data assimilation for improved recovery of the underlying flow.

Fig. 2. — Slices horizontal velocity component of input flows for four different cases of the validation dataset and their denoising with $σ = 1.0$ and $λ = 0.5$ . From left to right: the noisy input flow, the *Baseline CNN* prediction, the PINN prediction, the ground truth data. First reference in section 3.1

Focusing on the noiseless training data set-up ( $σ = 0$ ) in Table 3, and comparing the results of the Baseline CNN and the λ = 0.0 PINN, we observe the PINN architecture perhaps has greater expressive power regardless of utilizing the NSE loss. Nonetheless, introducing the NSE loss ( $λ > 0$ ) enables the model to learn more general features across examples, and generalize better, increasing the reconstruction ability of the model up to a factor of 1.8. Note, even though the input is a solution to the NSE, increasing $λ$ too much leads to a lower performance on $Noiseless ℒ_{recon}$ . This is reasonable since the imposition of the PDE constraint reduces both the weight of the reconstruction loss and the space of functions the network can take. This reduces the performance on reconstruction, even though the two objectives are ostensibly aligned, suggesting that the CFD simulation, which satisfies the PDE over the volume, might have disparities when interpreted point-wise as data. This observation may also be a consequence of the downsamplnig step from the CFD data to the 3D images.

3.2. Solving the Navier-Stokes Equations

Here we explored the ability of the model to solve the NSE and recover the pressure field.

3.2.1. PDE Loss

The multi-optimization task leads the PINN to recover the pressure gradient field necessary to evaluate the NSE. We explore the $ℒ_{p}$ component, the RMSE for each component of the PDE, and the correlation coefficient of the predicted pressure field versus the ground truth as a function of $σ$ and $λ$ in Table 4.

Table 4.

Model performance on satisfying the PDE for stenosis and aneurysm flows for varying λ and σ. Each model was trained for 6000 epochs. The results were evaluated on the validation dataset.

σ	λ	$ℒ_{p}$	RMSE_continuity	RMSE_mtm,x	RMSE_mtm,y	RMSE_mtm,z	ρ_press

0.0	0.0	90.5	1.33e-1	3.07	1.40	8.90	0.10
	0.1	2.68e-3	3.77e-2	2.19e-2	2.14e-2	1.78e-2	0.873
	0.5	8.35e-4	2.01e-2	1.24e-2	1.25e-2	1.08e-2	0.967
	1.0	5.73e-4	1.67e-2	1.02e-2	1.08e-2	8.39e-3	0.941

1.0	0.0	9.41	1.58e-1	1.58	1.04	2.41	−0.155
	0.1	2.14e-3	3.76e-2	1.80e-2	1.80e-2	8.76e-3	0.966
	0.5	1.31e-3	2.51e-2	1.71e-2	1.66e-2	1.05e-2	0.971
	1.0	5.86e-4	1.53e-2	1.10e-2	1.13e-2	1.02e-2	0.946

Open in a new tab

Table 4 highlights that the PDE loss decreases as $λ$ increases, and $σ$ decreases. These results are expected since the importance of the PDE loss is increased with $λ$ , which forces the model to optimize for solving the PDE. Specifically, when $λ = 0.0$ , the predicted pressure field is not constrained in any way. Hence, the model randomly predicts values that lead to $ℒ_{p}$ being up to 4 orders of magnitude larger than the other cases, and a weak pressure field correlation coefficient $ρ_{press}$ . As soon as $λ > 0.0$ , the pressure field becomes constrained by the NSE loss, and the RMSE of the momentum equation is robustly satisfied in all directions. For instance with $σ = 1.0$ and $λ = 1.0$ , the momentum RMSE are 1.10e-2, 1.13e-2 and 1.02e-2 for the $x, y$ and $z$ directions respectively, which are all of the same order of magnitude. The continuity equation generally has a larger RMSE, of 1.53e-2. This could be due to both the implicit definition of pressure, and the irregularities in the boundary conditions that are overcome by viscous forces in the momentum equations. Moreover, we notice that the addition of noise increases the RMSE overall. This observation is due to the increased discrepancy between the two optimization tasks. Nonetheless, the overall values achieved are in the same order of magnitude as for the noiseless data, which demonstrates the robust denoising ability of the proposed method.

3.2.2. Pressure Field Recovery

The recovery of the pressure gradient field is inherent to the PINN framework. The correlation values $ρ_{press}$ in Table 4 for all 3 values of λ were all consistently high. We observed that the best correlation coefficient was reached for $λ = 0.5$ for both $σ = 0.0$ and $σ = 1.0$ . Likely, when $λ$ is too high, the model overly favors the solution of the PDE and effectively the model “solves” the NSE for a different set of boundary conditions than that of the true underlying flow. This results from the fact that the weight of the reconstruction loss on the total loss reduces, and that the boundaries are inherently inexact due to noise and resolution.

We display examples of normalized output and CFD pressure fields in Figure 3, which reveals their high correlations and that our model predicts the overall patterns of the pressure field. Notably, in the the stenosis, the adverse pressure gradient at the stenosis center is predicted whereas for the aneurysm, the circular high pressure zone at the entry (left) is predicted.

Fig. 3. — Normalized pressure field reconstruction for a stenosis and an aneurysm with $σ = 1.0$ and $λ = 0.5$ . From left to right: slice of a horizontal velocity component of an input flow, the pressure field predicted by the PINN, the ground truth pressure field. First reference in section 3.2.2

3.3. Wall Shear Stress

Our framework enables us to compute the wall shear stress directly from the trained model. To compute the wall shear stress, we need to evaluate the gradient and divergence of the velocity field, as well as the normal vector to the geometry. The signed distance field is a level set of the geometry, hence taking its gradient yields its normal vectors $n_{i}$ . We computed gradients by using an 8 points stencil central difference method, and alternatively, we obtained the velocity gradient and divergence by leveraging automatic differentiation of the PINN. To evaluate the wall shear stress, we first computed the viscous stress vector at the wall:

σ_{i} = [μ (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}) - \frac{2 μ}{3} \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}] n_{j}

where $δ_{i j}$ is the Kronecker delta. Since the zero divergence condition is not exactly satisfied by the soft constraint during training, we preserved the divergence term in the evaluation of $σ_{i}$ . The shear stress vector at the wall was then obtained by:

τ_{i} = σ_{i} - σ_{j} n_{j} n_{i}

We evaluated the magnitude of $τ_{i}$ at the wall and compared it with the CFD simulations for the model trained with $σ = 1.0$ and $λ = 0.5$ , which has the best pressure field correlation coefficient, on the validation set and obtained $ρ_{w s s} = 0.82$ . This coefficient is comparable to deep learning methods focused at predicting the wall shear stress directly [35] and reveals the advantage of solving the NSE with a conditioned PINN, for which the recovery of the wall shear stress is a byproduct. Figure 4 displays examples of normalized wall shear stress reconstruction for a stenosis and an aneurysm flow from the validation dataset. The zones of high and low shear stress are well captured.

Fig. 4. — Normalized wall shear stress magnitude reconstruction for a stenosis and an aneurysm with $σ = 1.0$ and $λ = 0.5$ . From left to right: slice of a horizontal velocity component of an input flow, the wall shear stress predicted by the PINN, the ground truth wall shear stress. First reference in section 3.3

3.4. Super-resolution

We evaluated the ability of our model trained with $σ = 1.0$ and $λ = 0.5$ to perform super-resolution in two ways: (1) super-resolve data from the same resolution than the training data, and (2) super-resolve data from lower resolution than the training data. In both cases, the super-resolution capability of the model stems from increasing the resolution of the images fed as an input. These images represent the MRV images to be processed. By processing low-resolution, noisy flow inputs and comparing them with high-resolution, noiseless outputs, our model effectively achieves super-resolution. The model in this section remains unchanged from its initial training with 64³ input images. The two autoencoders and the PINN are not re-trained but kept frozen with their previously acquired parameters.

3.4.1. Same Resolution Data

As motivated in section 2.1.2, we evaluated the ability of our model to super-resolve the input from ~ 25,000 to ~ 75,000 degrees of freedom. This represents a super-resolution 3x or ~ 1.44x per Cartesian direction. Such an ability is required to overcome the low resolution of common measurement. In this experiment, the input images are kept as the noisy 64³ images from the validation set. For the PINN, we predict the output field on a 92³ domain instead of 64³ to achieve the super-resolution. For the Baseline CNN, we perform a trilinear upsampling of the flow autoencoder output, from 64³ to 92³. Indeed, the Baseline CNN architecture does not allow one to directly get a 92³ output image. As shown in Table 5, we observe that the performance of the PINN decreases slightly reaching a PSNR of 37.9 vs. 38.3 for the constant resolution denoising task. Moreover, the RMSE_continuity increases slightly to 2.71e-2 vs 2.51e-2. Nonetheless, the model appears to be robust, specifically when compared to the Baseline CNN for which the PSNR drops from 34.1 to 33.3. Examples of the 3x super-resolution are shown in Figure 5, revealing the finer predictions of the model while corresponding to the ground truth.

Table 5.

Super-resolution performance of the PINN trained with σ = 1.0 and λ = 0.5. The results are evaluated on the high resolution validation dataset.

		Baseline CNN		PINN
Input	Prediction	${Noiseless ℒ}_{recon}$	PSNR	${Noiseless ℒ}_{recon}$	PSNR	ρ_pres

64	92	4.65e-4	33.3	1.63e-4	37.9	0.967

Open in a new tab

Fig. 5. — Slices of horizontal velocity component of input flows from the validation dataset, and their denoising/super-resolution. From left to right: the noisy input flow, the *Baseline CNN* prediction, the PINN prediction, the high resolution ground truth data. First reference in section 3.4.1

Moreover, we present visualization of super-resolution by factors [8, 64] in Figure 6 to demonstrate the ability of the model to perform higher super-resolution levels. In these cases, the RMSE_continuity increases to 2.76e-2, and remains of similar value, since it is a property of the model.

Fig. 6. — Slices of velocity components (top: horizontal, bottom:vertical) from the validation dataset, and their denoising/super-resolution with $σ = 1.0$ and $λ = 0.5$ . From left to right: the noisy input flow, the PINN x8 prediction, the PINN x64 prediction, the 64³ ground truth data. First reference in section 3.4.1

3.4.2. Lower Resolution Data

Acquiring high resolution 3D flow data is a slow and tedious process. However, such data can be obtained in limited quantities. Hence we explore the ability of our model, to leverage the training data of resolution 64³, to super-resolve lower resolution data at inference time. Such an ability could have crucial importance to get high resolution results in reasonable times from clinical scans. To this end, we preserve the model trained with $σ = 1.0$ and $λ = 0.5$ and generate lower resolution data by down-sampling the validation data to values in [64³; 12³] using tri-linear interpolation. We then add the noise on this lower resolution data, and perform a trilinear up-sampling to 64³ before passing it to the trained model. Therefore, the input images are still of shape 64³, but have lost some information due to the down-sampling step. Figure 7 exhibits examples of the new inputs as a function of their down-sampled shape.

Fig. 7. — Super-resolution ability of the PINN with respect to the input size. Top row: the ground truth flow and pressure with $σ = 1.0$ and $λ = 0.5$ . From left to right: the flow dimension, the upsampled noisy flow, the *Baseline CNN* prediction, the PINN flow prediction, the PINN pressure prediction. First reference in section 3.4.2

Table 6 displays the reconstruction performance of the PINN based on the dimension of the input. The noiseless reconstruction specifically drops for a super-resolution factor of 151.7 with a PSNR of 33.9, while remaining robust up to a scaling factor of 18.9, with a PSNR of 36.7 vs. 38.3 on the original data. For all resolutions, the PINN is seen to outperform the Baseline CNN method significantly, yielding a $Noiseless ℒ_{recon}$ twice smaller. Moreover, the prediction of the pressure field remains above 0.957, suggesting it is well captured in all cases. These results suggest that the method could be used to train a model on high resolution data to super-resolve low dimensional data efficiently. This ability could have particular importance in a clinical setting, where rare high resolution data could be leveraged to increase the resolution of lower resolution images. We exhibit a visual example of super-resolution in Figure 7. Independently of the input dimension, the difference in the predictions of the PINN are hardly visible, while the improvement of the baseline with input image size is observed when comparing the top and final rows.

Table 6.

Super-resolution performance of the PINN trained with σ = 1.0 and λ = 0.5 for input flows of lower resolution. The results are evaluated on the 64³ validation dataset.

Input	Scaling factor	Baseline CNN		PINN
Input	Scaling factor	${Noiseless ℒ}_{recon}$	PSNR	${Noiseless ℒ}_{recon}$	PSNR	ρ_pres

12	151.7	7.16e-4	31.4	4.06e-4	33.9	0.957
24	18.9	4.65e-4	33.3	2.15e-4	36.7	0.964
32	8.0	4.16e-4	33.8	1.89e-4	37.2	0.968
48	2.4	4.09e-4	33.9	1.64e-4	37.9	0.969
64	1	3.91e-4	34.1	1.48e-4	38.3	0.971

Open in a new tab

As a side note, since the flow autoencoder is fully convolutional, we also explored inputting the lower resolution 3D image directly into the network to generate the flow conditioning vector $γ_{flow}$ . However, the latter method proved to be less effective. We hypothesize that the model performs better on up-sampled data since the flow autoencoder features have been trained on that resolution. We link these observations to usual practices with widespread image processing CNNs such as the EfficientNets [36], for which the images are expected to be re-sampled to the training input size despite their fully convolutional nature.

3.5. Model Generalization

The performance of the model has been explored with average values. However, to explore the generalization ability of the model, we plot the RMSE_continuity as a function of the $R e$ and geometry, for both the models trained with $σ = 0.0$ and $σ = 1.0$ at $λ = 0.5$ in Figure 8.

Fig. 8. — RMSE_continuity as a function of $R e$ for stenosis and aneurysm flows with $λ = 0.5$ . From left to right: results for the noiseless setting $σ = 0.0$ , results for the noisy setting $σ = 1.0$ . First reference in section 3.5

We draw two main observations from Figure 8. Firstly, for the $σ = 0.0$ it appears that the RMSE_continuity is approximately independent of the geometry, whilst it slightly increases with $R e$ . For $σ = 1.0$ , the observations are not geometry independent. For the stenosis, the RMSE_continuity seems to increase linearly with $R e > 1000$ . For the aneurysm, the RMSE_continuity seems independent of $R e$ . The rationale behind these observations and the differences between the noiseless and noisy set-up stems from the multiplicative Gaussian noise. Indeed, a high $R e$ in a stenosis flow is linked to a high velocity at the center of the stenosis. Since the noise is proportional to the velocity, these flow sections will be more highly impacted by the noise. Consequently, the model performance on these regions will decrease. We do not observe such behaviours in aneurysms since the flow is not accelerated since the radius of the vessel increases and the flow is incompressible.

4. Discussion and Conclusions

We herein proposed and tested a physics-informed machine learning method for flow denoising and super-resolution. In particular, the proposed a method trains a single model that can be applied to a variety of flow and geometry conditions. We demonstrated the potential of this methodology to (1) denoise 3D stenosis and aneurysm flows of varying geometries and flow conditions, and (2) recover the hidden pressure field and predict the wall shear stress, and (3) super-resolve coarse and noisy flow image data to 20x resolution robustly. The model produces physically viable solutions, even without the exact knowledge of boundary conditions. This approach illustrates potential for broader applicability in flow measurement denoising and super-resolution.

To tackle the task of building a single model for varying geometries and flow conditions, we demonstrated that a latent representation of both the flow field and the geometry are efficient for learning a conditioned PDE solution. In this framework, our model performs reduced order modeling in a space that is unintelligible to humans. However, the latent space can be decoded back to the flow space through the PINN. Introducing the physics-informed loss with the mean squared error in the training process enabled the model to perform data assimilation whilst optimizing its weights to satisfy the conditioned NSE pointwise. The mean squared error of both the physics loss and noiseless flow reconstruction loss were observed to be of order of magnitude ≤ 1.5e − 3 and ≤ 1.5e − 4 respectively for 100% multiplicative Gaussian noise. These values confirm the model’s ability to effectively perform a multi-optimization task, and thereby of the ability of the NSE to perform denoising. In the process, the model is able to partially recover the associated pressure gradient field and wall shear stress. Increasing the weight parameter of the physics loss leads to an output that is more physically consistent but may deviate from the input data.

In the post-processing workflow of MRV images, the position of our method can be discussed. Firstly, the efficacy of our processing step is contingent upon the accurate labeling of flow voxels, necessitating a preliminary labelling stage. In this stage, readings from voxels outside the vascular domain are excluded and assigned a zero value to ensure that the subsequent analysis is confined to regions of actual flow. Secondly, the comparison of the results with $σ = 0.0$ and $σ = 1.0$ reveals that as the noise level in the input data decreases, the PINN recovery of the ground truth flow increases. This indicates that our method could be effectively integrated with existing denoising techniques, potentially enhancing the quality of image reconstruction from raw data. We acknowledge that our method is trained exclusively on data with multiplicative Gaussian noise, and does not consider other types noises that may be more aligned with MRV corruptions. However, we theorize that it could address errors such as phase wrapping and aliasing since it utilizes physics-informed principles that discern and rectify physically inconsistent patterns. Moreover, its capability for super-resolution, grounded in maintaining physical consistency, may allow for improved image detail and the potential mitigation of intra-voxel dephasing. Nonetheless, all of this needs to be tested. Although we do not apply our method to actual MRV measurements, nor do we consider the range of geometries or boundary conditions observable in cardiovascular flows, we do believe the range of flow conditions tested is sufficiently varied for proof of concept of the proposed methodology, and this preliminary success motivates application of this method to more complex and varied flow measurements.

Through data transformation, we indirectly impose rotational and translational equivariance in the model. Developing further research into equivariant neural networks could provide a methodology that inherently performs efficiently on such transformations. Nonetheless, since the gravity force acts in a single direction, rotations along all axes may not be sensible. Additionally, future works should focus on improving the robustness of the model to higher $R e$ flows. Similarly, exploring different activation functions could help improve performance on higher velocity flows. E.g., since tanh tends to saturate, the gradients on large values throughout the model could saturate as well, impacting the training process negatively. Finally, introducing unsteady data and uncertainty quantification in the model output are interesting research directions to build upon the method.

Acknowledgements

We thank Arjun Narayanan and Numi Sveinsson for insightful discussions. SCS acknowledges funding from the American Heart Association, Award No. 18EIA33900046.

Footnotes

• The authors do not have conflicts of interest relevant to this manuscript.

• Data: The script to generate the data using Simvascular and Python will be provided in https://github.com/Theosau/superresolution/tree/encoders_geometry_flow/data/stenosis_sv upon publication.

Code availability:

• The code will be made available at https://github.com/Theosau/superresolution/tree/encoders_geometry_flow upon publication.

References

[1].Updegrove A, Wilson NM, Merkow J, Lan H, Marsden AL, and Shadden SC, 2016. “SimVascular: An open source pipeline for cardiovascular simulation”. Annals of Biomedical Engineering, 45(3), Dec., pp. 525–541. [DOI] [PMC free article] [PubMed] [Google Scholar]
[2].Izzo R, Steinman D, Manini S, and Antiga L, 2018. “The Vascular Modeling Toolkit: A Python Library for the Analysis of Tubular Structures in Medical Images”. Journal of Open Source Software, 3(25), May, p. 745. [Google Scholar]
[3].Markl M, Frydrychowicz A, Kozerke S, Hope M, and Wieben O, 2012. “4D flow MRI”. Journal of Magnetic Resonance Imaging, 36(5), pp. 1015–1036. [DOI] [PubMed] [Google Scholar]
[4].Xiao H, Wu J-L, Wang J-X, Sun R, and Roy C, 2016. “Quantifying and reducing model-form uncertainties in reynolds-averaged navier–stokes simulations: A data-driven, physics-informed bayesian approach”. Journal of Computational Physics, 324, Nov., pp. 115–136. [Google Scholar]
[5].Sun L, Gao H, Pan S, and Wang J-X, 2019. “Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data”. Computer Methods in Applied Mechanics and Engineering, 361, 11, p. 112732. [Google Scholar]
[6].Brunton SL, Noack BR, and Koumoutsakos P, 2020. “Machine learning for fluid mechanics”. Annual Review of Fluid Mechanics, 52(1), Jan., pp. 477–508. [Google Scholar]
[7].Bakhshinejad A, Baghaie A, Vali A, Saloner D, Rayz VL, and D’Souza RM, 2017. “Merging computational fluid dynamics and 4d flow MRI using proper orthogonal decomposition and ridge regression”. Journal of Biomechanics, 58, June, pp. 162–173. [DOI] [PMC free article] [PubMed] [Google Scholar]
[8].Fathi MF, Bakhshinejad A, Baghaie A, Saloner D, Sacho RH, Rayz VL, and D’Souza RM, 2018. “Denoising and spatial resolution enhancement of 4d flow MRI using proper orthogonal decomposition and lasso regularization”. Computerized Medical Imaging and Graphics, 70, Dec., pp. 165–172. [DOI] [PubMed] [Google Scholar]
[9].Töger J, Zahr MJ, Aristokleous N, Bloch KM, Carlsson M, and Persson P-O, 2020. “Blood flow imaging by optimal matching of computational fluid dynamics to 4d-flow data”. Magnetic Resonance in Medicine, 84(4), Apr., pp. 2231–2245. [DOI] [PubMed] [Google Scholar]
[10].De Avila Belbute-Peres F, Economon T, and Kolter Z, 2020. “Combining differentiable PDE solvers and graph neural networks for fluid flow prediction”. In Proceedings of the 37th International Conference on Machine Learning, III HD and Singh A, eds., Vol. 119 of Proceedings of Machine Learning Research, PMLR, pp. 2402–2411. [Google Scholar]
[11].Cai S, Mao Z, Wang Z, Yin M, and Karniadakis GE, 2021. “Physics-informed neural networks (PINNs) for fluid mechanics: a review”. Acta Mechanica Sinica, 37(12), Dec., pp. 1727–1738. [Google Scholar]
[12].Wang R, 2021. “Physics-guided deep learning for dynamical systems: A survey”. ArXiv. [DOI] [PMC free article] [PubMed] [Google Scholar]
[13].Arzani A, Wang J-X, Sacks MS, and Shadden SC, 2022. “Machine learning for cardiovascular biomechanics modeling: Challenges and beyond”. Annals of Biomedical Engineering, 50(6), pp. 615–627. [DOI] [PubMed] [Google Scholar]
[14].Subramaniam A, Wong ML, Borker RD, Nimmagadda S, and Lele SK, 2020. “Turbulence enrichment using physics-informed generative adversarial networks”. arXiv. [Google Scholar]
[15].Gao H, Sun L, and Wang J-X, 2021. “Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels”. Physics of Fluids, 33(7), p. 073603. [Google Scholar]
[16].Raissi M, Perdikaris P, and Karniadakis G, 2019. “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”. Journal of Computational Physics, 378, Feb., pp. 686–707. [Google Scholar]
[17].Oldenburg J, Borowski F, Öner A, Schmitz K-P, and Stiehm M, 2022. Geometry aware physics informed neural network surrogate for solving navier-stokes equation (GAPINN), Mar. [Google Scholar]
[18].Kingma DP, and Welling M, 2022. “Auto-encoding variational bayes”. arXiv. [Google Scholar]
[19].Maleki A, Heyse J, Ranade R, He H, Kasimbeg P, and Pathak J, 2021. “Geometry encoding for numerical simulations”. arXiv. [Google Scholar]
[20].Ranade R, Hill C, Ghule L, and Pathak J, 2022. “A composable machine-learning approach for steady-state simulations on high-resolution grids”. arXiv. [Google Scholar]
[21].Rutkowski D, Roldán-Alzate A, and Johnson K, 2021. “Enhancement of cerebrovascular 4d flow mri velocity fields using machine learning and computational fluid dynamics simulation data”. Scientific Reports, 11, 05. [DOI] [PMC free article] [PubMed] [Google Scholar]
[22].Callaghan FM, and Grieve SM, 2016. “Spatial resolution and velocity field improvement of 4d-flow MRI”. Magnetic Resonance in Medicine, 78(5), Nov., pp. 1959–1968. [DOI] [PubMed] [Google Scholar]
[23].Zhuang B, Sirajuddin A, Zhao S, and Lu M, 2021. “The role of 4d flow MRI for clinical applications in cardiovascular disease: current status and future perspectives”. Quantitative Imaging in Medicine and Surgery, 11(9), Sept., pp. 4193–4210. [DOI] [PMC free article] [PubMed] [Google Scholar]
[24].Cherry M, Khatir Z, Khan A, and Bissell M, 2022. “The impact of 4d-flow MRI spatial resolution on patientspecific CFD simulations of the thoracic aorta”. Scientific Reports, 12(1), Sept. [DOI] [PMC free article] [PubMed] [Google Scholar]
[25].Ferdian E, Suinesiaputra A, Dubowitz DJ, Zhao D, Wang A, Cowan B, and Young AA, 2020. “4dflownet: Super-resolution 4d flow MRI using deep learning and computational fluid dynamics”. Frontiers in Physics, 8, May. [Google Scholar]
[26].Jiang C, Esmaeilzadeh S, Azizzadenesheli K, Kashinath K, Mustafa M, Tchelepi HA, Marcus P, Prabhat M, and Anandkumar A, 2020. “MESHFREEFLOWNET: a physics-constrained deep continuous space-time super-resolution framework”. In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. [Google Scholar]
[27].Fathi MF, Perez-Raya I, Baghaie A, Berg P, Janiga G, Arzani A, and D’Souza RM, 2020. “Super-resolution and denoising of 4d-flow MRI using physics-informed deep neural nets”. Computer Methods and Programs in Biomedicine, 197, Dec., p. 105729. [DOI] [PubMed] [Google Scholar]
[28].Kissas G, Yang Y, Hwuang E, Witschey WR, Detre JA, and Perdikaris P, 2020. “Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow MRI data using physics-informed neural networks”. Computer Methods in Applied Mechanics and Engineering, 358, Jan., p. 112623. [Google Scholar]
[29].Goodman JW, 1976. “Some fundamental properties of speckle∗”. J. Opt. Soc. Am., 66(11), Nov, pp. 1145–1150. [Google Scholar]
[30].Frank S, Lee J, Lantz J, Ebbers T, and Shadden S, 2021. “Cardiac kinetic energy and viscous dissipation rate from radial flow data”. Frontiers in Physiology, 12, 09. [DOI] [PMC free article] [PubMed] [Google Scholar]
[31].Huang H, Hu X, Zhao Y, Makkie M, Dong Q, Zhao S, Li K, and Liu T, 2017. “Modeling task fmri data via deep convolutional autoencoder”. IEEE Transactions on Medical Imaging, PP, 06, pp. 1–1. [DOI] [PubMed] [Google Scholar]
[32].Cheng Z, Sun H, Takeuchi M, and Katto J, 2018. “Deep convolutional autoencoder-based lossy image compression”. CoRR, abs/1804.09535. [Google Scholar]
[33].Shit S, Zimmermann J, Ezhov I, Paetzold JC, Sanches AF, Pirkl C, and Menze BH, 2022. “SRflow: Deep learning based super-resolution of 4d-flow MRI data”. Frontiers in Artificial Intelligence, 5, Aug. [DOI] [PMC free article] [PubMed] [Google Scholar]
[34].López MM, Frederick JM, and Ventura J, 2021. “Evaluation of MRI denoising methods using unsupervised learning”. Frontiers in Artificial Intelligence, 4, June. [DOI] [PMC free article] [PubMed] [Google Scholar]
[35].Ferdian E, Dubowitz DJ, Mauger CA, Wang A, and Young AA, 2022. “WSSNet: Aortic wall shear stress estimation using deep learning on 4d flow MRI”. Frontiers in Cardiovascular Medicine, 8, Jan. [DOI] [PMC free article] [PubMed] [Google Scholar]
[36].Tan M, and Le Q, 2019. “EfficientNet: Rethinking model scaling for convolutional neural networks”. In Proceedings of the 36th International Conference on Machine Learning, Chaudhuri K and Salakhutdinov R, eds., Vol. 97 of Proceedings of Machine Learning Research, pp. 6105–6114. [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Data Availability Statement

• The code will be made available at https://github.com/Theosau/superresolution/tree/encoders_geometry_flow upon publication.

PERMALINK

Unsupervised denoising and super-resolution of vascular flow data by physics-informed machine learning

Théophile Sautory

Shawn C Shadden

Abstract

Purpose:

Methods:

Results:

Conclusion:

1. Introduction

2. Methods

2.1. Dataset

2.1.1. Flow data

2.1.2. Generating low resolution image data

2.1.3. Transformations

2.1.4. Signed Distance Fields

Fig. 1.

2.1.5. Normalization

2.1.6. Noise

2.2. Machine Learning Framework

2.2.1. Overview

2.2.2. Geometry Autoencoder

Table 1.

2.2.3. Flow Autoencoder

Table 2.

2.2.4. Physics-Informed Neural Network

Reconstruction Loss

Navier-Stokes Equation Loss

2.2.5. Evaluation Metrics

3. Results

3.1. Denoising

Table 3.

Fig. 2.

3.2. Solving the Navier-Stokes Equations

3.2.1. PDE Loss

Table 4.

3.2.2. Pressure Field Recovery

Fig. 3.

3.3. Wall Shear Stress

Fig. 4.

3.4. Super-resolution

3.4.1. Same Resolution Data

Table 5.

Fig. 5.

Fig. 6.

3.4.2. Lower Resolution Data

Fig. 7.

Table 6.

3.5. Model Generalization

Fig. 8.

4. Discussion and Conclusions

Acknowledgements

Footnotes

Code availability:

References

Associated Data

Data Availability Statement

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases