1 Introduction

Angiographic imaging constitutes a crucial area of medical imaging. It led, in particular, to the development of devoted acquisition protocols in Magnetic Resonance Imaging, X-ray Computed Tomography, or Ultra-Sound imaging. Within this wide spectrum of modalities and specific sequences, Magnetic Resonance Angiography (MRA) is particularly used, for instance in the case of detection and diagnosis of vascular alterations (e.g., stenoses, aneurysms, thromboses, ...), and patient follow-up after treatment (e.g., by stents, coils, ...). MRA is also an effective tool for clinical research devoted to the understanding of human macro and mesovascular physiology, since it provides non-invasive/non-ionising modalities (e.g., Phase Contrast (PC) or Time-of-Flight (TOF) MRA), available in 2D, 3D, and even 3D+time.

Due to the nature of structures visualized in MRA, namely thin, elongated, curvy vessels, the development of specific image analysis tools has been an active research area during the last decades [1]. However, contrary to other morphological structures (e.g., brain tissues and structures, with BrainWeb [2]), there does not exist any common framework to facilitate the development and validation of related image analysis methods. In particular, despite a few attempts (such as MICCAI challenges [3]), that rely on datasets equiped with manual ground truth, there does not exist any “vascular analogue” of BrainWeb, i.e., a framework for generating efficiently virtual MRA images naturally equiped with a ground-truth and/or associated to a real MRA.

Stemming on this fact, we have been working, in the context of the VIVABRAIN ProjectFootnote 1 with the final goal to tackle this issue by providing an efficient and versatile methodological framework allowing for the generation of virtual MRAs from real MRAs. The notion of virtual angiography is not new, and was already developed for CT Angiography simulation [4]. However, the case of MRA was not intensively considered, in particular via the efficient coupling of CFD and MRI simulation.

This article presents our framework. In Sect. 2, we present the proposed methodological pipeline. In Sect. 3, we focus on the most challenging issues raised by the instanciation of this pipeline, and some of our approaches for tackling them. In Sect. 4, we describe experimental results that validate and/or illustrate our main contributions.

2 Overview of the Framework

Our framework relies on a pipeline composed of five main steps, going from real MRA to virtual MRA. Four of these steps require to tackle methodological issues, and are related to three main disciplinary fields, namely computer science (image processing and analysis), mathematics (CFD), and physics (MRI simulation). The first step is devoted to MRA acquisition. The second aims at segmenting vascular volumes from such data, in order to further define mixed 2D/3D computational meshes of vascular structures. The third step consists of designing a vascular model, in particular by exp- licitly modeling anatomical and physiological features. The fourth step aims to simulate the flowing blood in the designed vascular model, leading to velocity and pressure fields. The fifth step consists of generating MRA images based on the previously computed information, in order to produce virtual MRA data.

figure a

This methodological pipeline is materialized by a software framework, namely AngioTkFootnote 2, which can embed plug-ins implementing some of the above methodological steps. The data flowchart associated to this pipeline is the following: 3D real MRAs \(\rightarrow \) 3D digital volumes \(\rightarrow \) 2D/3D computational meshes \(\rightarrow \) 3D / 3D+T scalar and vectorial fields \(\rightarrow \) 3D virtual MRAs. Providing some of them as open-data is a complementary purpose of the software framework, itself developed in an open-source way.

3 Challenges and Principal Contributions

In this section, we discuss some of the principal contributions within the proposed framework, in relation with the four aforementioned methodological steps.

3.1 Segmentation of Vascular Volumes

The very first issue raised by the instanciation of the proposed framework is related to the efficient extraction of vascular volumes from 3D MRA images, i.e., the segmentation of vessels. Many approaches have been developed so far, for this specific topic. In particular, deformable model paradigms (e.g., level-sets) have been widely investigated. These approaches were mainly motivated by the dual handling of data-driven information (namely the hypersignal of blood) and a regularization term (related to the geometry of vessels). However, the use of a 2D model (either implicit or explicit) for handling structures mainly 1-dimensional generally tends to decrease the robustness of such approaches.

Considering the relevance of the dual handling of two, data- and feature-based terms, we propose to express the problem of vessel segmentation as a variational problem. As for deformable models, this leads to minimize a two-term energy; however, this formulation no longer consists of mapping a geometric model onto the sought vascular structures. In particular, we developed a Chan-Vese formulation [5], in which the regularization term is directly driven by a vesselness measure.

In this context, we investigated two kinds of vesselness measures, namely Frangi Vesselness, and an alternative morphological feature, namely RORPOFootnote 3 [6], proved to present complementary properties in terms of accuracy and precision. In addition, such an approach can be enhanced with directional information in the regularization term. Results of filtering with RORPO and segmentation with our Chan-Vese approach are provided in Sect. 4.1.

3.2 Modeling of Vascular Networks

The digital vascular volumes can be processed to design 2D/3D computational mesh models. Many works have been devoted to this issue, and we did not investigate alternative approaches. In particular, we rely on freely available existing toolsFootnote 4. We focus here on the issue of accurately modeling the physiological hypotheses associated to the considered vascular networks. To illustrate our modus operandi, we consider the case study of the cerebral venous network. At the macroscopic scale, it is composed by input veins (7–11) draining the blood into sinuses (2,3) until their confluence (4). The blood then passes into lateral sinuses (5,6) and reaches an extracranial area, composed of the internal jugular veins (1).

figure b

We assume that: (i) the blood density is constant; (ii) the flow is incompressible and isothermal; (iii) the Newtonian model is used for blood flow. The latter assumption consists of neglecting shear thinning and viscoelastic effects [7]. The cranium is considered as a rigid closed box, and the brain tissues and cerebrospinal fluid contain mainly water, thus constituting an incompressible tissue. We then suppose that vessel walls are rigid. Patient-specific blood flow data are usually not collected from routine clinical examinations. Hence, scarce literature data are used for the velocity magnitude in cerebral veins [8]. The above hypotheses and the governing parameters computed from the literature (Reynold: 75 \(\rightarrow \) 1055; Stokes: 1.10 \(\rightarrow \) 3.84; Strouhal: 0.013 \(\rightarrow \) 0.030) lead to model the venous blood flow dynamics by the Navier-Stokes (NS) equations

$$\begin{aligned} \left\{ \begin{array}{rcll} \rho (\partial _t \mathbf {u} + (\mathbf u \cdot \nabla ) \mathbf u ) - \mu \varDelta \mathbf u + \nabla p &{} = &{} \mathbf f &{} \text {in}\; \varOmega \times [0,T] \\ \text{ div }~\mathbf u &{} = &{} 0 &{} \text {in}\; \varOmega \times [0,T] \\ \mathbf u \vert _{t=0} &{} = &{} \mathbf {u_0} &{} \text {in}\; \varOmega \\ \end{array} \right. \end{aligned}$$
(1)

where \(\mathbf u \): fluid velocity, p: pressure, \(\rho \): density, \(\mu \): dynamic viscosity, \(\mathbf f \): an applied body force. For boundary conditions, we use a constant velocity of small magnitude at the inflow (coming from microcirculation), and free traction boundary condition for the outflow. The no-slip condition is imposed on the wall boundaries, assumed to be rigid.

3.3 Simulation of Flowing Blood

In order to simulate flowing blood in complex networks, reduced order models are often considered [9]. However, this is not sufficient to accurately capture subtle behaviours, in particular for MRA simulation. Then, we aim to generalize the 3D CFD approaches generally considered for vessel samples [10] to the case of more complete networks.

The two numerical methods we use to compute an approximate solution to the NS equations are finite element methods. Both are available as open-source (FreeFem++/Feel++)Footnote 5. FreeFem++ uses characteristics method to handle the non-linear term, while Feel++ uses the Oseen linearization. Both methods were validatedFootnote 6 with analytical solutions, e.g., Ethier-Steinman. The finite element spaces used to discretize the velocity and the pressure are the Inf-Sup stable Taylor-Hood finite elements P2P1. In order to speed-up the resolution, this method is coupled with an iterative Uzawa conjuguate gradient algorithm with a Cahouet-Chabart preconditionner. For Feel++ simulations, we used a second order finite difference scheme to approximate the time derivative and a second order extrapolation of the nonlinear convective term. The resolution in the latter is handled by a LU solver. The discretized NS equations were supplemented for both Feel++ and Freefem++ simulations by a no-slip boundary conditions on the lateral surface of the computational domain, since we are dealing with viscous fluid. At the inflow and the outflow, we imposed a Dirichlet and a Neumann boundary conditions.

3.4 MRA Simulation

From the 3D geometry of the vascular network, and the velocity field obtained from the NS equations resolution, MRA can finally be simulated. The sample to be imaged is divided into equal subvolumes called isochromats, assumed to possess uniform physical properties (spin relaxation times T1, T2, T2*, equilibrium magnetization \(\mathbf {M_{0}}, \ldots \)). The measurable signal emitted by one isochromat is obtained by numerically solving the Bloch equations, that give the temporal evolution of tissue magnetization for one isochromat

$$\begin{aligned} \frac{d\mathbf {M}}{dt} = \gamma \mathbf {M} \times \mathbf {B} - \mathbf {\hat{R}}(\mathbf {M} - \mathbf {M_{0}}) \end{aligned}$$
(2)

with \(\mathbf {M}\) the magnetization vector of the tissue, \(\gamma \) the gyromagnetic ratio of hydrogen, \(\mathbf {B}\) the external magnetic field, \(\mathbf {\hat{R}}\) the relaxation matrix containing T1 and T2. The whole MR signal is obtained by summing the contribution of each isochromat over the sample.

To handle movement of spins over the time, we consider the Lagrangian approach. This requires to determine each individual spin trajectory. While solving Bloch equations, we need not use a different treatment for static tissues and flowing particles. We vary the position of the spin over time, which changes the field value seen by the particle

$$\begin{aligned} \mathbf {B}(\mathbf r ,t) = [\mathbf {G}(t).\mathbf {r}(t) + \varDelta B(\mathbf r ,t) ] .\mathbf {e_{z}} + \mathbf {B_{1}}(\mathbf r ,t) \end{aligned}$$
(3)

where the magnetic field term \(\mathbf {B}\) contains all the MR sequence elements (gradients and RF pulses): \(\mathbf {G}\) the gradient sequence, \(\mathbf {r}\) the isochromat position, \(\varDelta B\) the field inhomogeneities, \(\mathbf {B_{1}}\) the RF pulses sequence.

Specific software packages have been designed for MRI simulation, including JEMRISFootnote 7, which is an advanced MRI simulator software written in C++, open-source and freely modifiable. Natively, Bloch equation solving in JEMRIS is only dedicated to simulate static tissues; indeed, JEMRIS only allows to specify one trajectory for all spins. Then, we added a specific class to the C++ code to allow users to specify a different trajectory for each spin, thus implementing an MRA simulation functionality [11].

4 Experiments and Results

4.1 Validation of Image Processing Methods

In terms of image processing and analysis, two validations were carried out. The first is devoted to RORPO by comparing the efficiency of the associated vesselness feature vs. the gold standard of Frangi vesselness, in terms of image filtering (namely, vessel enhancement). Some qualitative results for filtering of a 3D MRA of brain vasculature are illustrated in Fig. 1. The efficiency of RORPO vs. Frangi vesselness, was also emphasised by a comparative threshold-based segmentation study, whose results are provided in the ROC curves, on the right (mean value and standard deviation).

figure c

The second validations are devoted to the Chan-Vese segmentation method embedding a vesselness measure as regularization term. In order to quantitatively assess the robustness of this approach on complex vascular networks, we applied this method on the DRIVE database [12], that is composed of 2D images (our approach is valid both for 2D and 3D), and comes with accurate ground truths. The results obtained on this database present an accuracy score of 0.9434, compared to 0.9442 for the state of the art [12] and 0.9479 for Human observer. These results, very close to the state of the art, prove the relevance of the approach, in particular for a method that was neither designed nor parametrically tuned for this application.

Fig. 1.
figure 1

Left: TOF MRA volume rendering (from low intensities, in red, to high intensities, in yellow). Center: RORPO filtering. Right: Frangi Vesselness. (Color figure online)

Full size image
Fig. 2.
figure 2

Physical phantom, designed as a double bifurcation fluid circuit, to reproduce acceleration, deceleration and recirculation all of which can be encountered in real vascular networks.

Full size image

4.2 Validation of CFD Methods

Various experiments were performed to assess the relevance of the proposed CFD methods. First, the numerical results were compared to those actually measured from MRA images of a physical phantom (Fig. 2). A pulsatile flow, obtained with MRA measurements, was used to impose the velocity at the inlet of the phantom. In order to carry out a cross-validation of the two finite elements methods, we compared velocity (and pressure) profiles at six different radial sections previously defined, and along the centerline. Two of these six profiles are presented on the right. Most of the results are equivalent in this Feel++–Freefem++ comparison. The input and output flows are very similar, both between the two simulations and between simulations and experiment. Branches flow are equally similar for the two simulations, while MRA measurements are slightly higher. These differences can be explained by the fact that MRA images were segmented beforehand, and this process tends to overestimate the real conditions of the experiment. In a more qualitative way, we also run a simplified NS problem on a realistic cerebral venous network with both Feel++ and Freefem++. The results very close from each other are illustrated in Fig. 3.

figure d
Fig. 3.
figure 3

Blood flow simulation in a cerebral venous network: streamlines and velocity fields.

Full size image

4.3 Validation of MRA Simulation Methods

MRI simulations were performed using JEMRIS with constant flows computed from numerical simulations and pulsatile ones, and on rigid and nonrigid phantoms. Only a part of these experiments is reported here; see [11, 13] for complementary details.

Fig. 4.
figure 4

Virtual MRA obtained from the velocity field computed by the CFD methods, and the MRA simulation implemented in JEMRIS. (Phantom geometry superimposed, in black.)

Full size image
figure e

In Fig. 4, we show a simulated phase contrast image of constant flow obtained with JEMRIS. About 19 000 spin trajectories were calculated on a time interval of 5 s from numerical data calculated for constant flow with Feel++ and Freefem++. We used a phase contrast sequence with resolution 128, matrix \(180\,\times \,30\), TE of 8 ms, TR of 100 ms, Venc of 400 mm/s and Nex = 1. Computations took 15 min using 20 CPUs. We find the initial phantom geometry as expected; we can notice a lack of spins in the lower branch due to the low number of flow particles travelling along that path. The velocities measured in each branch are on par, magnitude-wise, with the initial data provided by Feel++ and Freefem++. The low resolution used for the MRI simulation leads to an averaging of the velocities in the voxel, and consequently to an underestimation of the peak velocity in each branch of the phantom (partial volume effect). It is worth mentioning that the proposed framework, although preliminary, is already efficient enough to allow for simulating subtle MRA details, and in particular to reproduce specific acquisition artifacts. This is, for instance, the case for vascular signal shift due to TE delay (first row: real data; second row: simulated data).

5 Conclusion

Beyond specific contributions to the state of the art in vessel segmentation, blood flow simulation and MRI simulation, with the purpose of providing methods both fairly validated and openly available to the whole community, the purpose of this work was to provide a full methodological framework allowing to generate virtual MRAs from real ones, in a versatile and efficient paradigm. A software pipeline is associated to this framework, and will hopefully constitute a strong basis for facilitating research activities of the scientific community on angiographic imaging, and MRA in particular.

From a methodological point of view, we plan to further embed some solutions for generating vascular atlases [14] in this framework, thus allowing to go a step further by parameterizing the vessel geometric and topological properties. The generation of vascular atlases will however require to tackle open issues, in particular with regard to vessel registration. This will constitute our immediate further works.