Phantom project Alexandre Ancel 2 Alexandre Fortin 1 Simon Garnotel 3 - PowerPoint PPT Presentation

Phantom project Alexandre Ancel 2 Alexandre Fortin 1 Simon Garnotel 3 Olivia Miraucourt 1 Stéphanie Salmon 1 Ranine Tarabay 2 1 University of Reims Champagne-Ardenne, Reims, France 2 University of Strasbourg, IRMA / UMR 7501, Strasbourg, France 3 University of Picardie Jules Verne, BioFlowImage laboratory, Amiens, France August 25th 2015 A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 1 / 23

RMA Medical teams Phantom project goal: Validation of the CFD simulations: Virtual images Segmentation Cross validation simulation Computer science Feel++/FreeFem++ Physics teams teams Validation with experimental data Blood flow simulation Validation of the MRI Mathematics teams simulations Comparison with the initial MRI Apply the task loop on the Phantom Figure : The VIVABRAIN project task loop A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 2 / 23

Phantom What is a phantom ? A device designed to reproduce some features of flows, compatible with the MRI. Figure : Physical phantom for cerebral arteries. A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 3 / 23

Feel++/FreeFem++ cross validation Outline Feel++/FreeFem++ cross validation 1 Numerical methods Fluid simulation Results Feel++ and FreeFem++ comparison with experimental data 2 MRI simulation results 3 AngioTk pipeline 4 Conclusion 5 A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 4 / 23

Feel++/FreeFem++ cross validation Numerical methods Mathematical model for blood flow simulations Blood : homogeneous, incompressible fluid, with “standard” Newtonian behaviour, Mathematical model: unsteady Navier-Stokes equations: ρ∂ u ∂ t − 2 ∇ · ( µ D ( u )) + ρ ( u · ∇ ) u + ∇ p = f , dans Ω × I 0 , dans Ω × I ∇ · u = + initial conditions + boundary conditions Ω ⊂ R d ( d ≥ 2 ) : domain u : viscosity of the fluid ; p : pressure of the fluid ; D ( u ) = 1 2 ( ∇ u + ∇ u T ) : deformation tensor ; ρ and µ density and dynamic viscosity of the fluid. A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 5 / 23

Feel++/FreeFem++ cross validation Numerical methods Feel++ FE method: Oseen scheme Let’s consider a Dirichlet condition at the inflow and a Neumann BC at the outflow: on Γ in , (1) u = u in σ ( u , p ) n = g on Γ out . (2) For V = { v ∈ [ H 1 (Ω)] d | v = 0 on Γ w , v = u in on Γ in } M = L 2 and 0 (Ω) the variational formulation reads: find ( u , p ) ∈ V × M such that ∀ q ∈ M , ∀ v ∈ { v ∈ [ H 1 (Ω)] d | v = 0 on Γ w ∪ Γ in } , we have: � ρ∂ u � � � � ∂ t v + ρ ( u · ∇ u ) · v + 2 µ D ( u ) : D ( v ) dx − gn · v ds − p div ( v ) dx = Ω Ω Ω Γ out Ω � q div ( u ) dx = Ω (Ω ( h , k geo ) )] d × P N Space discretisation: Taylor-Hood finite element [ P N + 1 c (Ω ( h , k geo ) ) c Time discretisation: Finite difference order 2 Convective term treatment: Extrapolation of order 2 A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 6 / 23

Feel++/FreeFem++ cross validation Numerical methods FreeFem++ FE method: Method of characteristics For every particle, we write: � dX dt ( x , s ; t ) = u ( t , X ( x , s ; t )) (5) X ( x , s ; s ) = x where X ( x , s ; t ) is the particle position at time t who was in x at time s . That gives: ( ∂ t u + ( u . ∇ ) u )( t n , x ) ∼ u ( t n + 1 , x ) − u ( t n , X n ( x )) (6) dt with X n ( x ) = x − u ( t n , x ) . dt + O ( dt 2 ) . We finally have: ρ dt ( u n + 1 − u n ◦ X n ) − µ ∆ u n + 1 + ∇ p n + 1 = f n + 1 (7a) div ( u n + 1 ) = 0 (7b) This method is implemented using the convect operator of FreeFem++. A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 7 / 23

Feel++/FreeFem++ cross validation Fluid simulation Results Benchmark setup Figure : Radial slices where the velocity profiles are plotted Geometry Tetrahedrons DOF h min h max h avrg 2 · 10 − 1 5 · 10 − 1 3 · 10 − 1 M1 157 , 245 769 , 662 2 . 5 · 10 − 1 6 . 25 · 10 − 1 3 . 75 · 10 − 1 M2 93 , 655 469 , 008 3 · 10 − 1 7 . 5 · 10 − 1 4 . 5 · 10 − 1 M3 60 , 349 307 , 510 Table : The characteristics of the three types of geometries A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 8 / 23

Feel++/FreeFem++ cross validation Fluid simulation Results Constant Poiseuille flow Figure : Velocity magnitude, M3 mesh, constant flow Velocity and pressure magnitude along the centerline: Pressure along z axis Velocity along z axis 40 200 FreeFem M3 FreeFem M2 FreeFem M1 150 30 Feel M3 velocity profile velocity profile Feel M2 100 Feel M1 20 50 10 0 0 − 50 0 20 40 60 80 100 120 140 0 20 40 60 80 100 120 140 arcLength arcLength A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 9 / 23 Figure : Comparison Feel vs FreeFem on the M3 mesh with a constant flow ( V min ), P2P1:

Feel++/FreeFem++ cross validation Fluid simulation Results Velocity profile at the inlet section Velocity profile at the outlet section Velocity profile at the lower left section 150 Feel M3 80 150 Feel M2 Feel M1 60 FreeFem M3 100 velocity profile velocity profile velocityprofile FreeFem M2 100 FreeFem M1 40 Feel M3 Feel M3 50 50 Feel M2 Feel M2 20 Feel M1 Feel M1 FreeFem M3 FreeFem M3 FreeFem M2 FreeFem M2 0 0 0 FreeFem M1 FreeFem M1 − 3 − 2 − 1 0 1 2 3 − 2 − 1 0 1 2 − 2 − 1 . 5 − 1 − 0 . 5 arcLength arcLength arcLength Velocity profile at the upper left section Velocity profile at the lower right section Velocity profile at the upper right section 200 100 150 80 150 velocity profile velocity profile velocity profile 60 100 100 40 Feel M3 Feel M3 Feel M3 Feel M2 Feel M2 Feel M2 50 50 Feel M1 Feel M1 Feel M1 20 FreeFem M3 FreeFem M3 FreeFem M3 FreeFem M2 FreeFem M2 FreeFem M2 0 0 0 FreeFem M1 FreeFem M1 FreeFem M1 0 0 . 5 1 1 . 5 2 2 . 5 3 − 2 − 1 . 5 − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 arcLength arcLength arcLength Figure : Feel++ vs FreeFem++ comparison on the M3 mesh with a constant flow ( V min ), P2P1: Velocity profile at the right and left sections in the upper and lower channels A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 10 / 23

Feel++ and FreeFem++ comparison with experimental data Outline Feel++/FreeFem++ cross validation 1 Numerical methods Fluid simulation Results Feel++ and FreeFem++ comparison with experimental data 2 MRI simulation results 3 AngioTk pipeline 4 Conclusion 5 A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 11 / 23

Feel++ and FreeFem++ comparison with experimental data Comparison of the numerical outputs with respect to the experimental measurements: (pulsatile flow) Velocity profile at the inlet section Velocity profile at the outlet section 2 , 500 2 , 500 MRI MRI FreeFem FreeFem Feel Feel 2 , 000 2 , 000 flow profile flow profile 1 , 500 1 , 500 1 , 000 1 , 000 500 0 100 200 300 400 500 0 100 200 300 400 500 time time Velocity profile at the lower left section Velocity profile at the upper left section 600 MRI MRI FreeFem FreeFem 500 Feel Feel 1 , 500 400 flow profile flow profile 300 1 , 000 200 100 500 0 100 200 300 400 500 0 100 200 300 400 500 time time A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 12 / 23 Figure : Comparison Feel vs FreeFem on the M3 mesh with a pulsatl flow, flow profile at the

MRI simulation results Outline Feel++/FreeFem++ cross validation 1 Numerical methods Fluid simulation Results Feel++ and FreeFem++ comparison with experimental data 2 MRI simulation results 3 AngioTk pipeline 4 Conclusion 5 A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 13 / 23

MRI simulation results Mathematical model for MRI simulations In MRI, the signal collected over time is generated by the temporal variations of the macroscopic magnetization of tissues. This signal contains all the information needed to reconstruct the final image. The most popular technique for MRI simulation is isochromat summation. The sample to be imaged is divided into equal subvolumes called isochromats, see Fig. 8. Those subvolumes are supposed to possess uniform physical properties: spin relaxation times T 1, T 2, T 2*, equilibrium magnetization M 0 and magnetic susceptibility χ . Figure : Cutting the sample into isochromats. A magnetization vector is associated to each isochromat and its evolution is monitored during the acquisition sequence. The collected MR signal corresponds to the transverse component of the magnetization. A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay Phantom project 14 / 23