Coil sensitivities: proof of concept

KomaMRIBase/src/datatypes/Phantom.jl

@@ -47,12 +47,8 @@ julia> obj.ρ
     T2s::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000
     #Off-resonance related
     Δw::AbstractVector{T} = zeros(eltype(x), size(x))
-    #χ::Vector{SusceptibilityModel}
-    #Diffusion
-    Dλ1::AbstractVector{T} = zeros(eltype(x), size(x))
-    Dλ2::AbstractVector{T} = zeros(eltype(x), size(x))
-    Dθ::AbstractVector{T}  = zeros(eltype(x), size(x))
-    #Diff::Vector{DiffusionModel}  #Diffusion map
+    #EXPERIMENTAL: Coils
+    coil_sens::AbstractVector{T} = ones(eltype(x), length(x))
     #Motion
     motion::MotionModel{T} = NoMotion{eltype(x)}()
 end
@@ -123,7 +119,7 @@ end
 """
     obj = heart_phantom(...)
 
-Heart-like LV 2D phantom. The variable `circumferential_strain` and `radial_strain` are for streching (if positive) 
+Heart-like LV 2D phantom. The variable `circumferential_strain` and `radial_strain` are for streching (if positive)
 or contraction (if negative). `rotation_angle` is for rotation.
 
 # Arguments
@@ -159,11 +155,11 @@ function heart_phantom(
     ring = ⚪(R) .- ⚪(r)
     ρ = ⚪(r) .+ 0.9 * ring #proton density
     # Diffusion tensor model
-    D = 2e-9 #Diffusion of free water m2/s
-    D1, D2 = D, D / 20
-    Dλ1 = D1 * ⚪(R) #Diffusion map
-    Dλ2 = D1 * ⚪(r) .+ D2 * ring #Diffusion map
-    Dθ = atan.(x, -y) .* ring #Diffusion map
+    # D = 2e-9 #Diffusion of free water m2/s
+    # D1, D2 = D, D / 20
+    # Dλ1 = D1 * ⚪(R) #Diffusion map
+    # Dλ2 = D1 * ⚪(r) .+ D2 * ring #Diffusion map
+    # Dθ = atan.(x, -y) .* ring #Diffusion map
     T1 = (1400 * ⚪(r) .+ 1026 * ring) * 1e-3 #Myocardial T1
     T2 = (308 * ⚪(r) .+ 42 * ring) * 1e-3 #T2 map [s]
     # Generating Phantoms
@@ -174,9 +170,9 @@ function heart_phantom(
         ρ=ρ[ρ .!= 0],
         T1=T1[ρ .!= 0],
         T2=T2[ρ .!= 0],
-        Dλ1=Dλ1[ρ .!= 0],
-        Dλ2=Dλ2[ρ .!= 0],
-        Dθ=Dθ[ρ .!= 0],
+        # Dλ1=Dλ1[ρ .!= 0],
+        # Dλ2=Dλ2[ρ .!= 0],
+        # Dθ=Dθ[ρ .!= 0],
         motion=SimpleMotion(
             PeriodicHeartBeat(;
                 period=period,
@@ -225,7 +221,7 @@ julia> plot_phantom_map(obj, :ρ)
 ```
 """
 function brain_phantom2D(; axis="axial", ss=4, us=1)
-    # check and filter input    
+    # check and filter input
     ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, ss, us)
 
     # Get data from .mat file
@@ -322,7 +318,7 @@ end
     obj = brain_phantom3D(; ss=4, us=1)
 
 Creates a three-dimentional brain Phantom struct.
-Default ss=4 sample spacing is 2 mm. Original file (ss=1) sample spacing is .5 mm. 
+Default ss=4 sample spacing is 2 mm. Original file (ss=1) sample spacing is .5 mm.
 
 # References
 - B. Aubert-Broche, D.L. Collins, A.C. Evans: "A new improved version of the realistic
@@ -349,7 +345,7 @@ julia> plot_phantom_map(obj, :ρ)
 ```
 """
 function brain_phantom3D(; ss=4, us=1, start_end=[160, 200])
-    # check and filter input    
+    # check and filter input
     ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(3, ss, us)
 
     # Get data from .mat file
@@ -473,7 +469,7 @@ julia> pelvis_phantom2D(obj, :ρ)
 ```
 """
 function pelvis_phantom2D(; ss=4, us=1)
-    # check and filter input    
+    # check and filter input
     ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, ss, us)
 
     # Get data from .mat file
@@ -556,7 +552,7 @@ julia> ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(3, 4, [2, 2, 2])
 ```
 """
 function check_phantom_arguments(nd, ss, us)
-    # check for valid input    
+    # check for valid input
     ssz = -9999
     usz = -9999
     if length(us) > 1 || prod(us) > 1

KomaMRIBase/test/runtests.jl

@@ -371,11 +371,8 @@ end
     T2 = [0.09; 0.05; 0.04; 0.07; 0.005]
     T2s = [0.1; 0.06; 0.05; 0.08; 0.015]
     Δw = [-2e-6; -1e-6; 0.0; 1e-6; 2e-6]
-    Dλ1 = [-4e-6; -2e-6; 0.0; 2e-6; 4e-6]
-    Dλ2 = [-6e-6; -3e-6; 0.0; 3e-6; 6e-6]
-    Dθ = [-8e-6; -4e-6; 0.0; 4e-6; 8e-6]
-    obj1 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ)
-    obj2 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ)
+    obj1 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw)
+    obj2 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw)
     @test obj1 == obj2
 
     # Test size and length definitions of a phantom
@@ -558,7 +555,7 @@ end
     arbitrarymotion = @suppress ArbitraryMotion(t_start, t_end, 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt))
 
     # Test phantom subset
-    obs1 = Phantom(
+    obs1 = Phantom(;
         name,
         x,
         y,
@@ -568,26 +565,20 @@ end
         T2,
         T2s,
         Δw,
-        Dλ1,
-        Dλ2,
-        Dθ,
-        simplemotion
+        motion=simplemotion
     )
     rng = 1:2:5
-    obs2 = Phantom(
+    obs2 = Phantom(;
         name,
-        x[rng],
-        y[rng],
-        z[rng],
-        ρ[rng],
-        T1[rng],
-        T2[rng],
-        T2s[rng],
-        Δw[rng],
-        Dλ1[rng],
-        Dλ2[rng],
-        Dθ[rng],
-        simplemotion[rng],
+        x=x[rng],
+        y=y[rng],
+        z=z[rng],
+        ρ=ρ[rng],
+        T1=T1[rng],
+        T2=T2[rng],
+        T2s=T2s[rng],
+        Δw=Δw[rng],
+        motion=simplemotion[rng],
     )
     @test obs1[rng] == obs2
     @test @view(obs1[rng]) == obs2
@@ -598,26 +589,23 @@ end
     # @test @view(obs1[rng]) == obs2
 
     # Test addition of phantoms
-    oba = Phantom(
+    oba = Phantom(;
         name,
-        [x; x[rng]],
-        [y; y[rng]],
-        [z; z[rng]],
-        [ρ; ρ[rng]],
-        [T1; T1[rng]],
-        [T2; T2[rng]],
-        [T2s; T2s[rng]],
-        [Δw; Δw[rng]],
-        [Dλ1; Dλ1[rng]],
-        [Dλ2; Dλ2[rng]],
-        [Dθ; Dθ[rng]],
-        [obs1.motion; obs2.motion]
+        x=[x; x[rng]],
+        y=[y; y[rng]],
+        z=[z; z[rng]],
+        ρ=[ρ; ρ[rng]],
+        T1=[T1; T1[rng]],
+        T2=[T2; T2[rng]],
+        T2s=[T2s; T2s[rng]],
+        Δw=[Δw; Δw[rng]],
+        motion=[obs1.motion; obs2.motion]
     )
     @test obs1 + obs2 == oba
 
     # Test scalar multiplication of a phantom
     c = 7
-    obc = Phantom(name=name, x=x, y=y, z=z, ρ=c*ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ)
+    obc = Phantom(name=name, x=x, y=y, z=z, ρ=c*ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw)
     @test c * obj1 == obc
 
     #Test brain phantom 2D

KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl

@@ -63,7 +63,7 @@ function run_spin_precession!(
     M.xy .= Mxy[:, end]
     M.z  .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1))
     #Acquired signal
-    sig .= transpose(sum(Mxy[:, findall(seq.ADC)]; dims=1)) #<--- TODO: add coil sensitivities
+    sig .= transpose(sum(p.coil_sens .* Mxy[:, findall(seq.ADC)]; dims=1)) #<--- TODO: add coil sensitivities
 
     return nothing
 end

KomaMRIPlots/src/ui/DisplayFunctions.jl

@@ -1337,7 +1337,7 @@ function plot_phantom_map(
         l.width = width
     end
 
-    return Plot(trace, l, frames)
+    return plot_koma(trace, l, frames)
 end
 
 """

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 KomaMRI = "6a340f8b-2cdf-4c04-99be-4953d9b66d0a"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 PlutoSliderServer = "2fc8631c-6f24-4c5b-bca7-cbb509c42db4"
 

examples/3.tutorials/lit-06-Coils.jl

@@ -0,0 +1,63 @@
+# # Experimental: Simulating with realistic coils
+
+using KomaMRI, MAT # hide
+sys = Scanner() # hide
+
+# Loading coil sensitivities generated from electromagnetic simulation software.
+# To plot the coil sensitivities we will use the Phantom object:
+ph_file = joinpath(dirname(pathof(KomaMRI)), "../examples/2.phantoms/sphere_fields.mat")
+sphere = matread(ph_file)
+Δx = 1e-3                   # 5mm
+N = size(sphere["b1m"])  # Number of spins
+FOV = (N .- 1) .* Δx    # Field of view
+xr = -FOV[1]/2:Δx:FOV[1]/2  # x spin coordinates vector
+ρ = abs.(sphere["b1m"][:]) .> 0
+x = [x for (x, y, z) in Iterators.product(xr, xr, xr)][ρ .!= 0]
+y = [y for (x, y, z) in Iterators.product(xr, xr, xr)][ρ .!= 0]
+z = [z for (x, y, z) in Iterators.product(xr, xr, xr)][ρ .!= 0]
+coil_sens = 1.0 * sphere["b1m"][:][ρ .!= 0] ./ maximum(abs.(sphere["b1m"][:][ρ .!= 0]))
+ρ = 1.0 * ρ[ρ .!= 0]
+obj = Phantom(; x, y, z, ρ, coil_sens)
+p1 = plot_phantom_map(obj, :coil_sens ; height=400, width=400, darkmode=true)
+
+#md savefig(p1, "../assets/6-phantom1.html") # hide
+#jl display(p1)
+
+#md # ```@raw html
+#md # <center><object type="text/html" data="../../assets/6-phantom1.html" style="width:50%; height:420px;"></object></center>
+#md # ```
+
+# Now we will interpolate the coils into a brain phantom (this will be done internally):
+
+using Interpolations
+coil_sens = sphere["b1m"] ./ maximum(abs.(sphere["b1m"][:]))
+obj = brain_phantom2D()
+obj.coil_sens .= LinearInterpolation((xr,xr,xr), coil_sens, extrapolation_bc=0).(obj.x, obj.y, obj.z)
+p2 = plot_phantom_map(obj, :coil_sens ; height=400, width=400, darkmode=true)
+#md savefig(p2, "../assets/6-phantom2.html") # hide
+#jl display(p2)
+
+#md # ```@raw html
+#md # <center><object type="text/html" data="../../assets/6-phantom2.html" style="width:50%; height:420px;"></object></center>
+#md # ```
+
+# Then, we will load an EPI sequence.
+
+seq_file = joinpath(dirname(pathof(KomaMRI)), "../examples/5.koma_paper/comparison_accuracy/sequences/EPI/epi_100x100_TE100_FOV230.seq")
+seq = read_seq(seq_file)
+# And simulate:
+
+raw = simulate(obj, seq, sys) # hide
+acq = AcquisitionData(raw) # hide
+acq.traj[1].circular = false # hide
+Nx, Ny = raw.params["reconSize"][1:2] # hide
+reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx, Ny)) # hide
+image = reconstruction(acq, reconParams) # hide
+slice_abs = abs.(image[:, :, 1]) # hide
+p3 = plot_image(slice_abs; height=400) # hide
+#md savefig(p3, "../assets/6-recon.html") # hide
+#jl display(p3)
+
+#md # ```@raw html
+#md # <center><object type="text/html" data="../../assets/6-recon.html" style="width:65%; height:420px;"></object></center>
+#md # ```