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

1University of Reims Champagne-Ardenne, Reims, France 2University of Strasbourg, IRMA / UMR 7501, Strasbourg, France 3University 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

Blood flow simulation

Mathematics teams

Virtual images simulation

Physics teams

Segmentation

Computer science teams

Figure : The VIVABRAIN project task loop

Phantom project goal: Validation of the CFD simulations:

Cross validation Feel++/FreeFem++ Validation with experimental data

Validation of the MRI simulations

Comparison with the initial MRI

Apply the task loop on the Phantom

• 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

1

Feel++/FreeFem++ cross validation Numerical methods Fluid simulation Results

2

Feel++ and FreeFem++ comparison with experimental data

3

MRI simulation results

4

AngioTk pipeline

5

Conclusion

• 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 ∇ · u = 0, dans Ω × I + initial conditions + boundary conditions Ω ⊂ Rd(d ≥ 2): domain u : viscosity of the fluid ; p : pressure of the fluid ; D(u) = 1

2(∇u + ∇uT) : 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: u = uin

• n Γin,

(1) σ(u, p)n = g

• n Γout.

(2) For V = {v ∈ [H1(Ω)]d | v = 0 on Γw, v = uin on Γin} and M = L2

0(Ω) the

variational formulation reads: find (u, p) ∈ V × M such that ∀q ∈ M, ∀v ∈ {v ∈ [H1(Ω)]d | v = 0 on Γw ∪ Γin}, we have:

• Ω

ρ∂u ∂t v +

• Ω

ρ(u · ∇u) · v + 2µ

• Ω

D(u) : D(v) dx −

• Γout

gn · v ds −

• Ω

p div(v) dx =

• Ω

q div(u) dx = Space discretisation: Taylor-Hood finite element [PN+1

c

(Ω(h,kgeo))]d × PN

c (Ω(h,kgeo))

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)) X(x, s; s) = x (5) where X(x, s; t) is the particle position at time t who was in x at time s. That gives: (∂tu + (u.∇)u)(tn, x) ∼ u(tn+1, x) − u(tn, X n(x)) dt (6) with X n(x) = x − u(tn, x).dt + O(dt2). We finally have: ρ dt (un+1 − un ◦ X n) − µ∆un+1 + ∇pn+1 = f n+1 (7a) div(un+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 hmin hmax havrg Tetrahedrons DOF M1 2 · 10−1 5 · 10−1 3 · 10−1 157,245 769,662 M2 2.5 · 10−1 6.25 · 10−1 3.75 · 10−1 93,655 469,008 M3 3 · 10−1 7.5 · 10−1 4.5 · 10−1 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:

20 40 60 80 100 120 140 10 20 30 40 arcLength velocity profile Pressure along z axis FreeFem M3 FreeFem M2 FreeFem M1 Feel M3 Feel M2 Feel M1 20 40 60 80 100 120 140 −50 50 100 150 200 arcLength velocity profile Velocity along z axis

Figure : Comparison Feel vs FreeFem on the M3 mesh with a constant flow (Vmin), P2P1:

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 9 / 23

Feel++/FreeFem++ cross validation Fluid simulation Results

−3 −2 −1 1 2 3 20 40 60 80 arcLength velocityprofile Velocity profile at the inlet section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1 −2 −1 1 2 50 100 150 arcLength velocity profile Velocity profile at the outlet section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1 −2 −1.5 −1 −0.5 50 100 150 arcLength velocity profile Velocity profile at the lower left section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1 0.5 1 1.5 2 2.5 3 50 100 150 arcLength velocity profile Velocity profile at the upper left section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1 −2 −1.5 −1 −0.5 20 40 60 80 100 arcLength velocity profile Velocity profile at the lower right section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1 0.5 1 1.5 2 2.5 3 50 100 150 200 arcLength velocity profile Velocity profile at the upper right section Feel M3 Feel M2 Feel M1 FreeFem M3 FreeFem M2 FreeFem M1

Figure : Feel++ vs FreeFem++ comparison on the M3 mesh with a constant flow (Vmin), 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

1

Feel++/FreeFem++ cross validation Numerical methods Fluid simulation Results

2

Feel++ and FreeFem++ comparison with experimental data

3

MRI simulation results

4

AngioTk pipeline

5

Conclusion

• 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)

100 200 300 400 500 1,000 1,500 2,000 2,500 time flow profile Velocity profile at the inlet section MRI FreeFem Feel 100 200 300 400 500 500 1,000 1,500 2,000 2,500 time flow profile Velocity profile at the outlet section MRI FreeFem Feel 100 200 300 400 500 100 200 300 400 500 600 time flow profile Velocity profile at the lower left section MRI FreeFem Feel 100 200 300 400 500 500 1,000 1,500 time flow profile Velocity profile at the upper left section MRI FreeFem Feel

Figure : Comparison Feel vs FreeFem on the M3 mesh with a pulsatl flow, flow profile at the

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 12 / 23

MRI simulation results

Outline

1

Feel++/FreeFem++ cross validation Numerical methods Fluid simulation Results

2

Feel++ and FreeFem++ comparison with experimental data

3

MRI simulation results

4

AngioTk pipeline

5

Conclusion

• 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 T1, T2, T2*, equilibrium magnetization M0 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

MRI simulation results

Mathematical model: Bloch equations: Bloch equations : temporal evolution of magnetization dM dt = γM × B − ^ R(M − M0) (8) where M is the magnetization vector of the tissue, γ is the gyromagnetic ratio of hydrogen, B is the external magnetic field ^ R the relaxation matrix containing T1 and T2. The magnetic field term B(r, t) contains all the MR sequence elements (gradients and RF pulses). Its expression is given by: B(r, t) = [G(t).r + ∆B(r, t)].ez + B1(r, t) (9) where G(t) is the gradients sequence, r is the isochromat position, ∆B(r, t) is the field inhomogeneities due to off-resonance and non uniform gradients, B1(r, t) the RF pulses sequence.

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 15 / 23

MRI simulation results

JEMRIS results

(a) Amplitude (b) Phase

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 16 / 23

AngioTk pipeline

Outline

1

Feel++/FreeFem++ cross validation Numerical methods Fluid simulation Results

2

Feel++ and FreeFem++ comparison with experimental data

3

MRI simulation results

4

AngioTk pipeline

5

Conclusion

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 17 / 23

AngioTk pipeline

RMA

Medical teams

Blood flow simulation

Mathematics teams

Virtual images simulation

Physics teams

Segmentation

Computer science teams

Figure : The VIVABRAIN project task loop

Reproduce the vivabrain loop with a collection of software: Filtering: algorithms developed in the project Segmentation: algorithms developed in the project Mesh generation: VTK, vmtk Simulation: Feel++, FreeFrem++ MRA simulation: JEMRIS

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 18 / 23

AngioTk pipeline

Input data

Input data: Images in the DICOM format

Figure : DICOM image stack Figure : Volume rendering of DICOM data

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 19 / 23

AngioTk pipeline

AngioTK Pipeline: Mesh processing

Input Data Surface extraction Mesh refinement Mesh smoothing Centerlines extraction Reconstruction Surface opening Tetrahedralization Partitioning

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 20 / 23

AngioTk pipeline

(a) Thresholded MRI (b) Segmented MRI

Figure : Comparison of the realistic geometry with MRI segmentations, the fisrt is obtained by a simple threshold and the second by the classical method of snake.

Radial Slices IN UL LL UC LC UR LR OUT Realistic mesh 5.00 3.00 2.00 3.00 2.00 3.00 2.00 4.00 Thresholded MRI 3.12 2.70 1.23 2.45 1.73 2.70 1.97 3.36 Segmented MRI 4.30 5.23 3.30 2.20 5.28 5.04

Table : Comparison of the diameters (mm) at the radial slices : IN=inlet, UP=uppper left, LL=lower left, UC=upper center, UL=upper left, UR=upper right, LR=lower right and OUT=outlet. As the UL and LL parts and respectively the UR and LR parts are sticked together in the segmented IRM, we compute the diameters on the entire left slice and respectively on the entire right slice which are theorically equal to 5.11mm.

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 21 / 23

Conclusion

Outline

1

Feel++/FreeFem++ cross validation Numerical methods Fluid simulation Results

2

Feel++ and FreeFem++ comparison with experimental data

3

MRI simulation results

4

AngioTk pipeline

5

Conclusion

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 22 / 23

Conclusion

Conclusion

Cross validation Feel++/FremFem++: comparable results qualitatively and quantitatively Feel++/FreeFem++ vs experimental data: qualitatively ok, but not quantitatively

Improve MRA acquisition

First, virtual images from simulation

Geometry OK, but still issues with velocities

AngioTk pipeline: from MRA to computational meshes

Need tuning for the simulation step

• A. Ancel, A. Fortin, S. Garnotel, O. Miraucourt, S. Salmon, R. Tarabay

Phantom project 23 / 23