diff --git a/dev/AD/index.html b/dev/AD/index.html index 5e27112..5cf19dd 100644 --- a/dev/AD/index.html +++ b/dev/AD/index.html @@ -1,5 +1,5 @@ -Automatic Differentiation · EPGsim.jl
Edit on GitHub

Automatic differentiation

This page shows how to use Automatic Differentiation in combination with an EPG simulation.

The AD package tested is ForwardDiff.jl, maybe it works with others with some minor modification to the following code.

Load package

using EPGsim, ForwardDiff, CairoMakie

Building signal function

function MESE_EPG(T2,T1,TE,ETL,delta)
+Automatic Differentiation · EPGsim.jl
Edit on GitHub
Automatic differentiation
This page shows how to use Automatic Differentiation in combination with an EPG simulation. 
The AD package tested is ForwardDiff.jl, maybe it works with others with some minor modification to the following code.
Load package
using EPGsim, ForwardDiff, CairoMakie
Building signal function
function MESE_EPG(T2,T1,TE,ETL,delta)
   T = eltype(complex(T2))
   E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
   echo_vec = Vector{Complex{eltype(T2)}}()
@@ -63,7 +63,46 @@
 lines!(ax,TE_vec,abs.(amp))
 ax = Axis(f[1,2], title = "AD of MESE signal with B1 = $(deltaB1)",xlabel="TE [ms]")
 lines!(ax,TE_vec,df)
-f
Reproducibility
This page was generated with the following version of Julia:
using InteractiveUtils
+f
Differentiation along multiple variables
If we want to obtain the derivation along T1 and T2 we need to change the EPG_MESE function. The function should take as input a vector containing T2 and T1 (here noted T2/T1) :
function MESE_EPG2(T2T1,TE,ETL,delta)
+  T2,T1 = T2T1
+  T = complex(eltype(T2))
+  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
+  echo_vec = Vector{Complex{eltype(T2)}}()
+
+  E = epgRotation(E,pi/2*delta, pi/2)
+  ##loop over refocusing-pulses
+  for i = 1:ETL
+    E = epgDephasing(E,1)
+    E = epgRelaxation(E,TE,T1,T2)
+    E  = epgRotation(E,pi*delta,0.0)
+    E  = epgDephasing(E,1)
+    push!(echo_vec,E.Fp[1])
+  end
+
+  return abs.(echo_vec)
+end
+
+j2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])
50×2 Matrix{Float64}:
+ 0.00148849    0.0
+ 0.00267849    1.01626e-6
+ 0.0036074     3.11376e-9
+ 0.00424856    1.4929e-6
+ 0.0048583     2.1588e-7
+ 0.00505738    1.54224e-6
+ 0.00547205    4.77635e-7
+ 0.0053872     1.48187e-6
+ 0.00561831    5.69176e-7
+ 0.00541662    1.50582e-6
+ ⋮            
+ 0.000653146   6.78462e-7
+ 0.000604545  -4.24151e-7
+ 0.000549214   6.53958e-7
+ 0.000506571  -4.35399e-7
+ 0.000459289   6.2483e-7
+ 0.000424776  -4.38319e-7
+ 0.000383111   5.96642e-7
+ 0.000354969  -4.40833e-7
+ 0.000320195   5.7525e-7
Here we can see that the second column corresponding to T1 is equal to 0 which is expected for a MESE sequence and the derivative along T2 gives the same results :
j2[:,1] ≈ vec(j)
true
Reproducibility
This page was generated with the following version of Julia:
using InteractiveUtils
 io = IOBuffer();
 versioninfo(io);
 split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
@@ -71,10 +110,10 @@
  "Commit e4ee485e909 (2023-07-05 09:39 UTC)"
  "Platform Info:"
  "  OS: Linux (x86_64-linux-gnu)"
- "  CPU: 2 × Intel(R) Xeon(R) Platinum 8370C CPU @ 2.80GHz"
+ "  CPU: 2 × Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz"
  "  WORD_SIZE: 64"
  "  LIBM: libopenlibm"
- "  LLVM: libLLVM-14.0.6 (ORCJIT, icelake-server)"
+ "  LLVM: libLLVM-14.0.6 (ORCJIT, skylake-avx512)"
  "  Threads: 1 on 2 virtual cores"
  "Environment:"
  "  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager"
@@ -84,4 +123,4 @@
   [b6a82cc1] EPGsim v1.0.0-DEV `~/work/EPGsim.jl/EPGsim.jl`
   [f6369f11] ForwardDiff v0.10.35
   [98b081ad] Literate v2.14.0
-  [b77e0a4c] InteractiveUtils

+  [b77e0a4c] InteractiveUtils
diff --git a/dev/API/index.html b/dev/API/index.html index 4033522..c49f69f 100644 --- a/dev/API/index.html +++ b/dev/API/index.html @@ -1,2 +1,3 @@ -API · EPGsim.jl
Edit on GitHub
EPGsim.epgDephasing — Function
epgDephasing(E::EPGStates, n=1) where T

shifts the transverse dephasing states F corresponding to n dephasing-cycles.

source
EPGsim.epgRelaxation — Method
epgRelaxation(E::EPGStates,t,T1, T2)

applies relaxation matrices to a set of EPG states.

Arguments

  • E::EPGStates
  • t::AbstractFloat - length of time interval
  • T1::AbstractFloat - T1
  • T2::AbstractFloat - T2
source
EPGsim.epgRotation — Function
epgRotation(E::EPGStates, alpha::Float64, phi::Float64=0.0)

applies Bloch-rotation (<=> RF pulse) to a set of EPG states.

Arguments

  • E::EPGStates`
  • alpha::Float64 - flip angle of the RF pulse (rad)
  • phi::Float64=0.0 - phase of the RF pulse (rad)
source
EPGsim.rfRotation — Function
rfRotation_AT(alpha, phi=0.)

returns the rotation matrix for a pulse with flip angle alpha and phase phi.

source
+API · EPGsim.jl
Edit on GitHub
EPGsim.EPGStates — Type
EPGStates{T <: Real}

Stores the EPG states in 3 vectors Fp,Fn and Z.

Constructors :

EPGStates(Fp::Vector{Complex{S}},Fn::Vector{Complex{S}},Z::Vector{Complex{S}}) where {S <: Real}
+EPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number

Fields

  • Fp::Vector{Complex{T}}
  • Fn::Vector{Complex{T}}
  • Z::Vector{Complex{T}}

Related functions

  • getStates(E::EPGStates) : extract EPG states as matrix 3xN
source
EPGsim.epgDephasing — Function
epgDephasing(E::EPGStates, n=1) where T

shifts the transverse dephasing states F corresponding to n dephasing-cycles. n can be any integer

source
EPGsim.epgRelaxation — Method
epgRelaxation(E::EPGStates,t,T1, T2)

applies relaxation matrices to a set of EPG states.

Arguments

  • E::EPGStates
  • t::AbstractFloat - length of time interval
  • T1::AbstractFloat - T1
  • T2::AbstractFloat - T2
source
EPGsim.epgRotation — Function
epgRotation(E::EPGStates, alpha::Float64, phi::Float64=0.0)

applies Bloch-rotation (<=> RF pulse) to a set of EPG states.

Arguments

  • E::EPGStates`
  • alpha::Float64 - flip angle of the RF pulse (rad)
  • phi::Float64=0.0 - phase of the RF pulse (rad)
source
EPGsim.rfRotation — Function
rfRotation(alpha, phi=0.)

returns the rotation matrix for a pulse with flip angle alpha and phase phi.

Arguments

  • alpha - flip angle (radian)
  • phi=0. - phase of the flip angle (radian)
source
diff --git a/dev/index.html b/dev/index.html index 823e545..aa1d418 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · EPGsim.jl
Edit on GitHub

EPGsim

Extended Phase Graph simulation

Introduction

EPGsim is a Julia packet for magnetic resonance imaging signal simulation based on the Extended Phase Graph (EPG) concept. The principal aspect of this package was to make it compatible with Automatic Differentiation using ForwardDiff.jl in order to compute Cramér-Rao Lower Bound metrics which is used to optimized sequence protocol.

Note

EPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version

physics concept

Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :

Installation

This package is currently not registered.

Start julia and open the package mode by entering ]. Then enter

add https://github.com/aTrotier/EPGsim.jl

This will install EPGsim and all its dependencies. If you want to develop EPGsim itself you can checkout EPGsim by calling

dev https://github.com/aTrotier/EPGsim.jl

More information on how to develop a package can be found in the Julia documentation.

Tutorial

You can find an example about simulation of a Multi-Echo Spin-Echo sequence and its derivation here

+Home · EPGsim.jl
Edit on GitHub

EPGsim

Extended Phase Graph simulation

Introduction

EPGsim is a Julia packet for magnetic resonance imaging signal simulation based on the Extended Phase Graph (EPG) concept. The principal aspect of this package was to make it compatible with Automatic Differentiation using ForwardDiff.jl in order to compute Cramér-Rao Lower Bound metrics which is used to optimized sequence protocol.

Note

EPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version

EPG concept

Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :

Installation

This package is currently not registered.

Start julia and open the package mode by entering ]. Then enter

add https://github.com/aTrotier/EPGsim.jl

This will install EPGsim and all its dependencies. If you want to develop EPGsim itself you can checkout EPGsim by calling

dev https://github.com/aTrotier/EPGsim.jl

More information on how to develop a package can be found in the Julia documentation.

Tutorial

You can find an example about simulation of a Multi-Echo Spin-Echo sequence and its derivation here

diff --git a/dev/regular/index.html b/dev/regular/index.html index 170bda9..767cc58 100644 --- a/dev/regular/index.html +++ b/dev/regular/index.html @@ -8,9 +8,9 @@ E = EPGStates(T.([0.5+0.5im,1]),T.([0.5-0.5im,0]),T.([1,0]))
EPGStates{Float32}(ComplexF32[0.5f0 + 0.5f0im, 1.0f0 + 0.0f0im], ComplexF32[0.5f0 - 0.5f0im, 0.0f0 + 0.0f0im], ComplexF32[1.0f0 + 0.0f0im, 0.0f0 + 0.0f0im])
Warning

the F+[1] and F-[1] states should be complex conjugate and imag(Z[1])=0

EPG simulation

3 functions are used to simulate a sequence :

They take an EPGStates struct as first parameter.

E = EPGStates()
 E = epgRotation(E,deg2rad(60),0)
 E = epgDephasing(E,1)
-E = epgRotation(E,deg2rad(60),deg2rad(117))
EPGStates{Float64}(ComplexF64[0.38581714244240034 + 0.19658365292562008im, 0.0 - 0.649519052838329im], ComplexF64[0.38581714244240034 - 0.19658365292562008im, 0.17515731730550912 + 0.12725924011378179im], ComplexF64[0.2500000000000001 + 0.0im, 0.17024643740233 + 0.3341274465706379im])

Accessing states

States can seen directly as a vector :

E.Fp
2-element Vector{ComplexF64}:
+E = epgRotation(E,deg2rad(60),deg2rad(117))
EPGStates{Float64}(ComplexF64[0.38581714244240034 + 0.19658365292562008im, 0.0 - 0.649519052838329im], ComplexF64[0.38581714244240034 - 0.19658365292562008im, 0.17515731730550912 + 0.12725924011378179im], ComplexF64[0.2500000000000001 + 0.0im, 0.17024643740233 + 0.3341274465706379im])
Note

Currently, all the EPGstates are stored and used for calculation. The states equal or really close to zero are not deleted

Accessing states

States can seen directly as a vector :

E.Fp
2-element Vector{ComplexF64}:
  0.38581714244240034 + 0.19658365292562008im
                  0.0 - 0.649519052838329im

or by elements :

E.Fp[2]
0.0 - 0.649519052838329im

getStates is also available to create a 3xN matrix where 3 corresponds to Fp,Fn,Z and N is the number of states.

getStates(E)
3×2 Matrix{ComplexF64}:
  0.385817+0.196584im       0.0-0.649519im
  0.385817-0.196584im  0.175157+0.127259im
-     0.25+0.0im       0.170246+0.334127im
+ 0.25+0.0im 0.170246+0.334127im diff --git a/dev/search/index.html b/dev/search/index.html index 0a57c4a..029035d 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -Search · EPGsim.jl

Loading search...

    +Search · EPGsim.jl

    Loading search...

      diff --git a/dev/search_index.js b/dev/search_index.js index 2c2baf5..7057874 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"API/","page":"API","title":"API","text":"","category":"page"},{"location":"API/","page":"API","title":"API","text":"Modules = [EPGsim]","category":"page"},{"location":"API/#EPGsim.epgDephasing","page":"API","title":"EPGsim.epgDephasing","text":"epgDephasing(E::EPGStates, n=1) where T\n\nshifts the transverse dephasing states F corresponding to n dephasing-cycles.\n\n\n\n\n\n","category":"function"},{"location":"API/#EPGsim.epgRelaxation-Tuple{EPGStates, Any, Any, Any}","page":"API","title":"EPGsim.epgRelaxation","text":"epgRelaxation(E::EPGStates,t,T1, T2)\n\napplies relaxation matrices to a set of EPG states.\n\nArguments\n\nE::EPGStates\nt::AbstractFloat - length of time interval\nT1::AbstractFloat - T1\nT2::AbstractFloat - T2\n\n\n\n\n\n","category":"method"},{"location":"API/#EPGsim.epgRotation","page":"API","title":"EPGsim.epgRotation","text":"epgRotation(E::EPGStates, alpha::Float64, phi::Float64=0.0)\n\napplies Bloch-rotation (<=> RF pulse) to a set of EPG states.\n\nArguments\n\nE::EPGStates`\nalpha::Float64 - flip angle of the RF pulse (rad)\nphi::Float64=0.0 - phase of the RF pulse (rad)\n\n\n\n\n\n","category":"function"},{"location":"API/#EPGsim.rfRotation","page":"API","title":"EPGsim.rfRotation","text":"rfRotation_AT(alpha, phi=0.)\n\nreturns the rotation matrix for a pulse with flip angle alpha and phase phi.\n\n\n\n\n\n","category":"function"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"EPG implementation that mimics the regular implementation from Julien Lamy in Sycomore","category":"page"},{"location":"regular/#Short-description","page":"Regular EPG ","title":"Short description","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"Regular implementation use a constant positive or negative gradient dephasing. We use a vector Fp, Fn and Z to store the states.","category":"page"},{"location":"regular/#Initialization","page":"Regular EPG ","title":"Initialization","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"EPG states are stored as a structure :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"mutable struct EPGStates{T <: Real} \n Fp::Vector{Complex{T}}\n Fn::Vector{Complex{T}}\n Z::Vector{Complex{T}}\nend","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"which can be initialized with default parameters Fp = 0, Fn = 0 and Z = 1 states using :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"using EPGsim\nE = EPGStates()","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or by :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E = EPGStates(0,0,1)","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"which convert any numbers of the same types in Vector{ComplexF64}","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or directly by passing Vector{Complex{T}} where {T <: Real} which means it can accept a complex{dual} type :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"T = ComplexF32\nE = EPGStates(T.([0.5+0.5im,1]),T.([0.5-0.5im,0]),T.([1,0]))","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"warning: Warning\nthe F+[1] and F-[1] states should be complex conjugate and imag(Z[1])=0 ","category":"page"},{"location":"regular/#EPG-simulation","page":"Regular EPG ","title":"EPG simulation","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"3 functions are used to simulate a sequence :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"epgDephasing\nepgRelaxation\nepgRotation","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"They take an EPGStates struct as first parameter.","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E = EPGStates()\nE = epgRotation(E,deg2rad(60),0)\nE = epgDephasing(E,1)\nE = epgRotation(E,deg2rad(60),deg2rad(117))","category":"page"},{"location":"regular/#Accessing-states","page":"Regular EPG ","title":"Accessing states","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"States can seen directly as a vector :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E.Fp","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or by elements :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E.Fp[2]","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"getStates is also available to create a 3xN matrix where 3 corresponds to Fp,Fn,Z and N is the number of states.","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"getStates(E)","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = EPGsim","category":"page"},{"location":"#EPGsim","page":"Home","title":"EPGsim","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Extended Phase Graph simulation","category":"page"},{"location":"#Introduction","page":"Home","title":"Introduction","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"EPGsim is a Julia packet for magnetic resonance imaging signal simulation based on the Extended Phase Graph (EPG) concept. The principal aspect of this package was to make it compatible with Automatic Differentiation using ForwardDiff.jl in order to compute Cramér-Rao Lower Bound metrics which is used to optimized sequence protocol.","category":"page"},{"location":"","page":"Home","title":"Home","text":"note: Note\nEPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version","category":"page"},{"location":"#physics-concept","page":"Home","title":"physics concept","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :","category":"page"},{"location":"","page":"Home","title":"Home","text":"Lecture-04A: Definition of the Extended Phase Graph Basis\nLecture-04B: Sequence Operations in the Extended Phase Graph Domain\nLecture-04C: Examples using Extended Phase Graphs","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package is currently not registered.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Start julia and open the package mode by entering ]. Then enter","category":"page"},{"location":"","page":"Home","title":"Home","text":"add https://github.com/aTrotier/EPGsim.jl","category":"page"},{"location":"","page":"Home","title":"Home","text":"This will install EPGsim and all its dependencies. If you want to develop EPGsim itself you can checkout EPGsim by calling","category":"page"},{"location":"","page":"Home","title":"Home","text":"dev https://github.com/aTrotier/EPGsim.jl","category":"page"},{"location":"","page":"Home","title":"Home","text":"More information on how to develop a package can be found in the Julia documentation.","category":"page"},{"location":"#Tutorial","page":"Home","title":"Tutorial","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"You can find an example about simulation of a Multi-Echo Spin-Echo sequence and its derivation here","category":"page"},{"location":"AD/#Automatic-differentiation","page":"Automatic Differentiation","title":"Automatic differentiation","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"This page shows how to use Automatic Differentiation in combination with an EPG simulation. ","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"The AD package tested is ForwardDiff.jl, maybe it works with others with some minor modification to the following code.","category":"page"},{"location":"AD/#Load-package","page":"Automatic Differentiation","title":"Load package","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"using EPGsim, ForwardDiff, CairoMakie","category":"page"},{"location":"AD/#Building-signal-function","page":"Automatic Differentiation","title":"Building signal function","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"function MESE_EPG(T2,T1,TE,ETL,delta)\n T = eltype(complex(T2))\n E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])\n echo_vec = Vector{Complex{eltype(T2)}}()\n\n E = epgRotation(E,pi/2*delta, pi/2)\n # loop over refocusing-pulses\n for i = 1:ETL\n E = epgDephasing(E,1)\n E = epgRelaxation(E,TE,T1,T2)\n E = epgRotation(E,pi*delta,0.0)\n E = epgDephasing(E,1)\n push!(echo_vec,E.Fp[1])\n end\n\n return abs.(echo_vec)\nend;","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"warning: Specific types with AD\nForwardDiff use a specific type : Dual <: Real. The target function must be written generically enough to accept numbers of type T<:Real as input (or arrays of these numbers).We also need to create an EPGStates that is of that type. We need to force it to be complex :T = eltype(complex(T2))\nE = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])","category":"page"},{"location":"AD/#Define-parameters-for-simulation-and-run-it","page":"Automatic Differentiation","title":"Define parameters for simulation and run it","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"T2 = 60.0\nT1 = 1000.0\nTE = 7\nETL = 50\ndeltaB1 = 1\n\nTE_vec = range(TE,TE*ETL,ETL)\n\namp = MESE_EPG(T2,T1,TE,ETL,deltaB1)\nlines(TE_vec,abs.(amp),axis =(;title = \"MESE Signal\", xlabel=\"TE [ms]\"))","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"As expected, we get a standard T2 decaying exponential curve :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"S(TE) = exp(-TET_2)","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"We can analytically derive the equation according to T_2 :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"fracpartial Spartial T_2 = fracTET_2^2 exp(-TET_2)","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"which give the following curves:","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"df = TE_vec .* exp.(-TE_vec./T2)./(T2^2) \n\nlines(TE_vec,abs.(df),axis =(;title = \"dS/dT2\", xlabel=\"TE [ms]\"))","category":"page"},{"location":"AD/#Find-the-derivative-with-Automatic-Differentiation","page":"Automatic Differentiation","title":"Find the derivative with Automatic Differentiation","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Because we want to obtain the derivate at multiple time points (TE), we will use ForwardDiff.jacobian :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"j = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Let's compare it to the analytical equation :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"f=Figure()\nax = Axis(f[1,1],title =\"Analytic vs Automatic Differentiation\")\n\nlines!(ax,TE_vec,abs.(df),label = \"Analytic Differentiation\",linewidth=3)\nlines!(ax,TE_vec,abs.(vec(j)),label = \"Automatic Differentiation\",linestyle=:dash,linewidth=3)\naxislegend(ax)\nf","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Of course, in that case we don't really need the AD possibility. But if we reduce the B1+ value the equation becomes complicated enough and might lead to error during derivation if we don't use AD.","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"deltaB1 = 0.8\n\namp = MESE_EPG(T2,T1,TE,ETL,deltaB1)\nj = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])\n\nf = Figure()\nax = Axis(f[1,1], title = \"MESE signal with B1 = $(deltaB1)\",xlabel=\"TE [ms]\")\nlines!(ax,TE_vec,abs.(amp))\nax = Axis(f[1,2], title = \"AD of MESE signal with B1 = $(deltaB1)\",xlabel=\"TE [ms]\")\nlines!(ax,TE_vec,df)\nf","category":"page"},{"location":"AD/#Reproducibility","page":"Automatic Differentiation","title":"Reproducibility","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"This page was generated with the following version of Julia:","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"using InteractiveUtils\nio = IOBuffer();\nversioninfo(io);\nsplit(String(take!(io)), '\\n')","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"And with the following package versions","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"import Pkg; Pkg.status()","category":"page"}] +[{"location":"API/","page":"API","title":"API","text":"","category":"page"},{"location":"API/","page":"API","title":"API","text":"Modules = [EPGsim]","category":"page"},{"location":"API/#EPGsim.EPGStates","page":"API","title":"EPGsim.EPGStates","text":"EPGStates{T <: Real}\n\nStores the EPG states in 3 vectors Fp,Fn and Z.\n\nConstructors :\n\nEPGStates(Fp::Vector{Complex{S}},Fn::Vector{Complex{S}},Z::Vector{Complex{S}}) where {S <: Real}\nEPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number\n\nFields\n\nFp::Vector{Complex{T}}\nFn::Vector{Complex{T}}\nZ::Vector{Complex{T}}\n\nRelated functions\n\ngetStates(E::EPGStates) : extract EPG states as matrix 3xN\n\n\n\n\n\n","category":"type"},{"location":"API/#EPGsim.EPGStates-Union{Tuple{}, Tuple{T}, Tuple{T, T}, Tuple{T, T, T}} where T<:Number","page":"API","title":"EPGsim.EPGStates","text":"getStates(E::EPGStates)\n\nExtract EPG states as matrix 3xN\n\n\n\n\n\n","category":"method"},{"location":"API/#EPGsim.epgDephasing","page":"API","title":"EPGsim.epgDephasing","text":"epgDephasing(E::EPGStates, n=1) where T\n\nshifts the transverse dephasing states F corresponding to n dephasing-cycles. n can be any integer\n\n\n\n\n\n","category":"function"},{"location":"API/#EPGsim.epgRelaxation-Tuple{EPGStates, Any, Any, Any}","page":"API","title":"EPGsim.epgRelaxation","text":"epgRelaxation(E::EPGStates,t,T1, T2)\n\napplies relaxation matrices to a set of EPG states.\n\nArguments\n\nE::EPGStates\nt::AbstractFloat - length of time interval\nT1::AbstractFloat - T1\nT2::AbstractFloat - T2\n\n\n\n\n\n","category":"method"},{"location":"API/#EPGsim.epgRotation","page":"API","title":"EPGsim.epgRotation","text":"epgRotation(E::EPGStates, alpha::Float64, phi::Float64=0.0)\n\napplies Bloch-rotation (<=> RF pulse) to a set of EPG states.\n\nArguments\n\nE::EPGStates`\nalpha::Float64 - flip angle of the RF pulse (rad)\nphi::Float64=0.0 - phase of the RF pulse (rad)\n\n\n\n\n\n","category":"function"},{"location":"API/#EPGsim.rfRotation","page":"API","title":"EPGsim.rfRotation","text":"rfRotation(alpha, phi=0.)\n\nreturns the rotation matrix for a pulse with flip angle alpha and phase phi.\n\nArguments\n\nalpha - flip angle (radian)\nphi=0. - phase of the flip angle (radian)\n\n\n\n\n\n","category":"function"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"EPG implementation that mimics the regular implementation from Julien Lamy in Sycomore","category":"page"},{"location":"regular/#Short-description","page":"Regular EPG ","title":"Short description","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"Regular implementation use a constant positive or negative gradient dephasing. We use a vector Fp, Fn and Z to store the states.","category":"page"},{"location":"regular/#Initialization","page":"Regular EPG ","title":"Initialization","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"EPG states are stored as a structure :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"mutable struct EPGStates{T <: Real} \n Fp::Vector{Complex{T}}\n Fn::Vector{Complex{T}}\n Z::Vector{Complex{T}}\nend","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"which can be initialized with default parameters Fp = 0, Fn = 0 and Z = 1 states using :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"using EPGsim\nE = EPGStates()","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or by :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E = EPGStates(0,0,1)","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"which convert any numbers of the same types in Vector{ComplexF64}","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or directly by passing Vector{Complex{T}} where {T <: Real} which means it can accept a complex{dual} type :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"T = ComplexF32\nE = EPGStates(T.([0.5+0.5im,1]),T.([0.5-0.5im,0]),T.([1,0]))","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"warning: Warning\nthe F+[1] and F-[1] states should be complex conjugate and imag(Z[1])=0 ","category":"page"},{"location":"regular/#EPG-simulation","page":"Regular EPG ","title":"EPG simulation","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"3 functions are used to simulate a sequence :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"epgDephasing\nepgRelaxation\nepgRotation","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"They take an EPGStates struct as first parameter.","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E = EPGStates()\nE = epgRotation(E,deg2rad(60),0)\nE = epgDephasing(E,1)\nE = epgRotation(E,deg2rad(60),deg2rad(117))","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"note: Note\nCurrently, all the EPGstates are stored and used for calculation. The states equal or really close to zero are not deleted","category":"page"},{"location":"regular/#Accessing-states","page":"Regular EPG ","title":"Accessing states","text":"","category":"section"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"States can seen directly as a vector :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E.Fp","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"or by elements :","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"E.Fp[2]","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"getStates is also available to create a 3xN matrix where 3 corresponds to Fp,Fn,Z and N is the number of states.","category":"page"},{"location":"regular/","page":"Regular EPG ","title":"Regular EPG ","text":"getStates(E)","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = EPGsim","category":"page"},{"location":"#EPGsim","page":"Home","title":"EPGsim","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Extended Phase Graph simulation","category":"page"},{"location":"#Introduction","page":"Home","title":"Introduction","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"EPGsim is a Julia packet for magnetic resonance imaging signal simulation based on the Extended Phase Graph (EPG) concept. The principal aspect of this package was to make it compatible with Automatic Differentiation using ForwardDiff.jl in order to compute Cramér-Rao Lower Bound metrics which is used to optimized sequence protocol.","category":"page"},{"location":"","page":"Home","title":"Home","text":"note: Note\nEPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version","category":"page"},{"location":"#EPG-concept","page":"Home","title":"EPG concept","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :","category":"page"},{"location":"","page":"Home","title":"Home","text":"Lecture-04A: Definition of the Extended Phase Graph Basis\nLecture-04B: Sequence Operations in the Extended Phase Graph Domain\nLecture-04C: Examples using Extended Phase Graphs","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package is currently not registered.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Start julia and open the package mode by entering ]. Then enter","category":"page"},{"location":"","page":"Home","title":"Home","text":"add https://github.com/aTrotier/EPGsim.jl","category":"page"},{"location":"","page":"Home","title":"Home","text":"This will install EPGsim and all its dependencies. If you want to develop EPGsim itself you can checkout EPGsim by calling","category":"page"},{"location":"","page":"Home","title":"Home","text":"dev https://github.com/aTrotier/EPGsim.jl","category":"page"},{"location":"","page":"Home","title":"Home","text":"More information on how to develop a package can be found in the Julia documentation.","category":"page"},{"location":"#Tutorial","page":"Home","title":"Tutorial","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"You can find an example about simulation of a Multi-Echo Spin-Echo sequence and its derivation here","category":"page"},{"location":"AD/#Automatic-differentiation","page":"Automatic Differentiation","title":"Automatic differentiation","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"This page shows how to use Automatic Differentiation in combination with an EPG simulation. ","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"The AD package tested is ForwardDiff.jl, maybe it works with others with some minor modification to the following code.","category":"page"},{"location":"AD/#Load-package","page":"Automatic Differentiation","title":"Load package","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"using EPGsim, ForwardDiff, CairoMakie","category":"page"},{"location":"AD/#Building-signal-function","page":"Automatic Differentiation","title":"Building signal function","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"function MESE_EPG(T2,T1,TE,ETL,delta)\n T = eltype(complex(T2))\n E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])\n echo_vec = Vector{Complex{eltype(T2)}}()\n\n E = epgRotation(E,pi/2*delta, pi/2)\n # loop over refocusing-pulses\n for i = 1:ETL\n E = epgDephasing(E,1)\n E = epgRelaxation(E,TE,T1,T2)\n E = epgRotation(E,pi*delta,0.0)\n E = epgDephasing(E,1)\n push!(echo_vec,E.Fp[1])\n end\n\n return abs.(echo_vec)\nend;","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"warning: Specific types with AD\nForwardDiff use a specific type : Dual <: Real. The target function must be written generically enough to accept numbers of type T<:Real as input (or arrays of these numbers).We also need to create an EPGStates that is of that type. We need to force it to be complex :T = eltype(complex(T2))\nE = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])","category":"page"},{"location":"AD/#Define-parameters-for-simulation-and-run-it","page":"Automatic Differentiation","title":"Define parameters for simulation and run it","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"T2 = 60.0\nT1 = 1000.0\nTE = 7\nETL = 50\ndeltaB1 = 1\n\nTE_vec = range(TE,TE*ETL,ETL)\n\namp = MESE_EPG(T2,T1,TE,ETL,deltaB1)\nlines(TE_vec,abs.(amp),axis =(;title = \"MESE Signal\", xlabel=\"TE [ms]\"))","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"As expected, we get a standard T2 decaying exponential curve :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"S(TE) = exp(-TET_2)","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"We can analytically derive the equation according to T_2 :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"fracpartial Spartial T_2 = fracTET_2^2 exp(-TET_2)","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"which give the following curves:","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"df = TE_vec .* exp.(-TE_vec./T2)./(T2^2) \n\nlines(TE_vec,abs.(df),axis =(;title = \"dS/dT2\", xlabel=\"TE [ms]\"))","category":"page"},{"location":"AD/#Find-the-derivative-with-Automatic-Differentiation","page":"Automatic Differentiation","title":"Find the derivative with Automatic Differentiation","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Because we want to obtain the derivate at multiple time points (TE), we will use ForwardDiff.jacobian :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"j = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Let's compare it to the analytical equation :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"f=Figure()\nax = Axis(f[1,1],title =\"Analytic vs Automatic Differentiation\")\n\nlines!(ax,TE_vec,abs.(df),label = \"Analytic Differentiation\",linewidth=3)\nlines!(ax,TE_vec,abs.(vec(j)),label = \"Automatic Differentiation\",linestyle=:dash,linewidth=3)\naxislegend(ax)\nf","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Of course, in that case we don't really need the AD possibility. But if we reduce the B1+ value the equation becomes complicated enough and might lead to error during derivation if we don't use AD.","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"deltaB1 = 0.8\n\namp = MESE_EPG(T2,T1,TE,ETL,deltaB1)\nj = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])\n\nf = Figure()\nax = Axis(f[1,1], title = \"MESE signal with B1 = $(deltaB1)\",xlabel=\"TE [ms]\")\nlines!(ax,TE_vec,abs.(amp))\nax = Axis(f[1,2], title = \"AD of MESE signal with B1 = $(deltaB1)\",xlabel=\"TE [ms]\")\nlines!(ax,TE_vec,df)\nf","category":"page"},{"location":"AD/#Differentiation-along-multiple-variables","page":"Automatic Differentiation","title":"Differentiation along multiple variables","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"If we want to obtain the derivation along T1 and T2 we need to change the EPG_MESE function. The function should take as input a vector containing T2 and T1 (here noted T2/T1) :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"function MESE_EPG2(T2T1,TE,ETL,delta)\n T2,T1 = T2T1\n T = complex(eltype(T2))\n E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])\n echo_vec = Vector{Complex{eltype(T2)}}()\n\n E = epgRotation(E,pi/2*delta, pi/2)\n ##loop over refocusing-pulses\n for i = 1:ETL\n E = epgDephasing(E,1)\n E = epgRelaxation(E,TE,T1,T2)\n E = epgRotation(E,pi*delta,0.0)\n E = epgDephasing(E,1)\n push!(echo_vec,E.Fp[1])\n end\n\n return abs.(echo_vec)\nend\n\nj2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"Here we can see that the second column corresponding to T1 is equal to 0 which is expected for a MESE sequence and the derivative along T2 gives the same results :","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"j2[:,1] ≈ vec(j)","category":"page"},{"location":"AD/#Reproducibility","page":"Automatic Differentiation","title":"Reproducibility","text":"","category":"section"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"This page was generated with the following version of Julia:","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"using InteractiveUtils\nio = IOBuffer();\nversioninfo(io);\nsplit(String(take!(io)), '\\n')","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"And with the following package versions","category":"page"},{"location":"AD/","page":"Automatic Differentiation","title":"Automatic Differentiation","text":"import Pkg; Pkg.status()","category":"page"}] }