Add simulation of B1+-heterogeneity

KomaMRIBase/src/datatypes/Phantom.jl

@@ -1,5 +1,5 @@
 """
-    obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion)
+    obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion, B1)
 
 The Phantom struct. Most of its field names are vectors, with each element associated with
 a property value representing a spin. This struct serves as an input for the simulation.
@@ -18,6 +18,7 @@ a property value representing a spin. This struct serves as an input for the sim
 - `Dλ2`: (`::AbstractVector{T<:Real}`) spin Dλ2 (diffusion) parameter vector
 - `Dθ`: (`::AbstractVector{T<:Real}`) spin Dθ (diffusion) parameter vector
 - `motion`: (`::AbstractMotion{T<:Real}`) motion
+- `B1`: (`::AbstractVector{Complex{T<:Real}}`) spin transmit B1 parameter vector
 
 # Returns
 - `obj`: (`::Phantom`) Phantom struct
@@ -48,6 +49,7 @@ julia> obj.ρ
     #Diff::Vector{DiffusionModel}  #Diffusion map
     #Motion
     motion::AbstractMotion{T} = NoMotion{eltype(x)}() 
+    B1::AbstractVector{Complex{T}} = Complex.(ones(eltype(x), size(x)))
 end
 
 const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String,         fieldnames(Phantom))

KomaMRIBase/test/runtests.jl

@@ -381,8 +381,9 @@ end
     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θ)
+    B1 = Complex.([0.8; 0.9; 1.0; 1.1; 1.2])
+    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θ, B1=B1)
+    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θ, B1=B1)
     @test obj1 == obj2
 
     # Test size and length definitions of a phantom
@@ -578,7 +579,8 @@ end
         Dλ1,
         Dλ2,
         Dθ,
-        simplemotion
+        simplemotion,
+        B1
     )
     rng = 1:2:5
     obs2 = Phantom(
@@ -595,6 +597,7 @@ end
         Dλ2[rng],
         Dθ[rng],
         simplemotion[rng],
+        B1[rng]
     )
     @test obs1[rng] == obs2
     @test @view(obs1[rng]) == obs2
@@ -617,13 +620,14 @@ end
         [Dλ1; Dλ1[rng]],
         [Dλ2; Dλ2[rng]],
         [Dθ; Dθ[rng]],
-        vcat(obs1.motion, obs2.motion, length(obs1), length(obs2))
+        vcat(obs1.motion, obs2.motion, length(obs1), length(obs2)),
+        [B1; B1[rng]]
     )
     @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, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ, B1=B1)
     @test c * obj1 == obc
 
     #Test brain phantom 2D

KomaMRICore/src/datatypes/Spinor.jl

@@ -149,7 +149,7 @@ IEEE Transactions on Medical Imaging, 10(1), 53-65. doi:10.1109/42.75611
 
 # Arguments
 - `φ`: (`::Real`, `[rad]`) φ angle
-- `nxy`: (`::Real`) nxy factor
+- `nxy`: (`::Complex`) nxy factor
 - `nz`: (`::Real`) nz factor
 
 # Returns

KomaMRICore/src/simulation/SimMethods/Bloch/BlochCPU.jl

@@ -3,6 +3,7 @@ struct BlochCPUPrealloc{T} <: PreallocResult{T}
     M::Mag{T}                               # Mag{T}
     Bz_old::AbstractVector{T}               # Vector{T}(Nspins x 1)
     Bz_new::AbstractVector{T}               # Vector{T}(Nspins x 1)
+    B1::AbstractVector{Complex{T}}          # Vector{Complex{T}}(Nspins x 1)
     ϕ::AbstractVector{T}                    # Vector{T}(Nspins x 1)
     Rot::Spinor{T}                          # Spinor{T}
     ΔBz::AbstractVector{T}            # Vector{T}(Nspins x 1)
@@ -13,6 +14,7 @@ Base.view(p::BlochCPUPrealloc, i::UnitRange) = begin
         p.M[i],
         p.Bz_old[i],
         p.Bz_new[i],
+        p.B1[i],
         p.ϕ[i],
         p.Rot[i],
         p.ΔBz[i]
@@ -28,12 +30,13 @@ function prealloc(sim_method::Bloch, backend::KA.CPU, obj::Phantom{T}, M::Mag{T}
         ),
         similar(obj.x),
         similar(obj.x),
+        similar(obj.x, Complex{eltype(obj.x)}),
         similar(obj.x),
         Spinor(
             similar(M.xy),
             similar(M.xy)
         ),
-        obj.Δw ./ T(2π .* γ)
+        obj.Δw ./ T(2π .* γ)  # Why is this not similar(obj.x)?
     )
 end
 
@@ -126,6 +129,7 @@ function run_spin_excitation!(
 ) where {T<:Real}
     Bz = prealloc.Bz_old
     B = prealloc.Bz_new
+    B1 = prealloc.B1
     φ = prealloc.ϕ
     α = prealloc.Rot.α
     β = prealloc.Rot.β
@@ -139,12 +143,13 @@ function run_spin_excitation!(
         x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, s.t)
         #Effective field
         @. Bz = (s.Gx * x + s.Gy * y + s.Gz * z) + ΔBz - s.Δf / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z)
-        @. B = sqrt(abs(s.B1)^2 + abs(Bz)^2)
+        @. B1 = s.B1 * p.B1 # Take B1+ transmit RF field map into account. This is a vector of complex numbers
+        @. B = sqrt(abs(B1)^2 + abs(Bz)^2)
         @. B[B == 0] = eps(T)
         #Spinor Rotation
         @. φ = T(-π * γ) * (B * s.Δt) # TODO: Use trapezoidal integration here (?),  this is just Forward Euler 
         @. α = cos(φ) - Complex{T}(im) * (Bz / B) * sin(φ)
-        @. β = -Complex{T}(im) * (s.B1 / B) * sin(φ)
+        @. β = -Complex{T}(im) * (B1 / B) * sin(φ)
         mul!(Spinor(α, β), M, Maux_xy, Maux_z)
         #Relaxation
         @. M.xy = M.xy * exp(-s.Δt / p.T2)

KomaMRICore/src/simulation/SimMethods/Bloch/BlochGPU.jl

@@ -67,6 +67,7 @@ struct BlochGPUPrealloc{T} <: PreallocResult{T}
     seq_properties::AbstractVector{SeqBlockProperties{T}}
     Bz::AbstractMatrix{T}
     B::AbstractMatrix{T}
+    B1::AbstractMatrix{Complex{T}}
     φ::AbstractMatrix{T}
     ΔT1::AbstractMatrix{T}
     ΔT2::AbstractMatrix{T}
@@ -84,6 +85,7 @@ function prealloc_block(p::BlochGPUPrealloc{T}, i::Integer) where {T<:Real}
         [seq_block],
         view(p.Bz,:,1:seq_block.length),
         view(p.B,:,1:seq_block.length),
+        view(p.B1,:,1:seq_block.length),
         view(p.φ,:,1:seq_block.length-1),
         view(p.ΔT1,:,1:seq_block.length-1),
         view(p.ΔT2,:,1:seq_block.length-1),
@@ -102,6 +104,7 @@ function prealloc(sim_method::Bloch, backend::KA.GPU, obj::Phantom{T}, M::Mag{T}
         precalc.seq_properties,
         KA.zeros(backend, T, (size(obj.x, 1), max_block_length)),
         KA.zeros(backend, T, (size(obj.x, 1), max_block_length)),
+        KA.zeros(backend, Complex{T}, (size(obj.x, 1), max_block_length)),
         KA.zeros(backend, T, (size(obj.x, 1), max_block_length-1)),
         KA.zeros(backend, T, (size(obj.x, 1), max_block_length-1)),
         KA.zeros(backend, T, (size(obj.x, 1), max_block_length-1)),
@@ -186,7 +189,8 @@ function run_spin_excitation!(
 
     #Effective Field
     pre.Bz .= (x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz') .+ pre.ΔBz .- seq.Δf' ./ T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z)
-    pre.B .= sqrt.(abs.(seq.B1') .^ 2 .+ abs.(pre.Bz) .^ 2)
+    pre.B1 .= transpose(seq.B1) .* p.B1
+    pre.B .= sqrt.(abs.(pre.B1) .^ 2 .+ abs.(pre.Bz) .^ 2)
 
     #Spinor Rotation
     pre.φ .= T(-π .* γ) .* (pre.B[:,1:end-1] .* seq.Δt') # TODO: Use trapezoidal integration here (?),  this is just Forward Euler
@@ -196,7 +200,7 @@ function run_spin_excitation!(
     pre.ΔT2 .= exp.(-seq.Δt' ./ p.T2)
 
     #Excitation
-    apply_excitation!(backend, 512)(M.xy, M.z, pre.φ, seq.B1, pre.Bz, pre.B, pre.ΔT1, pre.ΔT2, p.ρ, ndrange=size(M.xy,1))
+    apply_excitation!(backend, 512)(M.xy, M.z, pre.φ, pre.B1, pre.Bz, pre.B, pre.ΔT1, pre.ΔT2, p.ρ, ndrange=size(M.xy,1))
     KA.synchronize(backend)
 
     #Reset Spin-State (Magnetization). Only for FlowPath

KomaMRICore/src/simulation/SimMethods/Bloch/KernelFunctions.jl

@@ -32,13 +32,13 @@ using KernelAbstractions: @kernel, @Const, @index, @uniform, @groupsize, @localm
         cos_φ = cos(φ[i_g, t])
         s_α_r[i_l] = cos_φ
         if (iszero(B[i_g, t]))
-            s_α_i[i_l] = -(Bz[i_g, t] / (B[i_g, t] + eps(T))) * sin_φ
-            s_β_r[i_l] = (imag(B1[t]) / (B[i_g, t] + eps(T))) * sin_φ
-            s_β_i[i_l] = -(real(B1[t]) / (B[i_g, t] + eps(T))) * sin_φ
+            s_α_i[i_l] = -(Bz[i_g, t] / (B[i_g, t] + eps(T))) * sin_φ # why not zero?
+            s_β_r[i_l] = (imag(B1[i_g, t]) / (B[i_g, t] + eps(T))) * sin_φ
+            s_β_i[i_l] = -(real(B1[i_g, t]) / (B[i_g, t] + eps(T))) * sin_φ
         else
             s_α_i[i_l] = -(Bz[i_g, t] / B[i_g, t]) * sin_φ
-            s_β_r[i_l] = (imag(B1[t]) / B[i_g, t]) * sin_φ
-            s_β_i[i_l] = -(real(B1[t]) / B[i_g, t]) * sin_φ
+            s_β_r[i_l] = (imag(B1[i_g, t]) / B[i_g, t]) * sin_φ
+            s_β_i[i_l] = -(real(B1[i_g, t]) / B[i_g, t]) * sin_φ
         end
             s_Mxy_new_r[i_l] = 2 * (s_Mxy_i[i_l] * (s_α_r[i_l] * s_α_i[i_l] - s_β_r[i_l] * s_β_i[i_l]) +
                                 s_Mz[i_l] * (s_α_i[i_l] * s_β_i[i_l] + s_α_r[i_l] * s_β_r[i_l])) +

KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl

@@ -86,11 +86,12 @@ function run_spin_excitation!(
         #Effective field
         ΔBz = p.Δw ./ T(2π .* γ) .- s.Δf ./ T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z)
         Bz = (s.Gx .* x .+ s.Gy .* y .+ s.Gz .* z) .+ ΔBz
-        B = sqrt.(abs.(s.B1) .^ 2 .+ abs.(Bz) .^ 2)
+        B1 = s.B1 .* p.B1 # Take B1+ transmit RF field map into account. This is complex
+        B = sqrt.(abs.(B1) .^ 2 .+ abs.(Bz) .^ 2)
         B[B .== 0] .= eps(T)
         #Spinor Rotation
         φ = T(-2π .* γ) .* (B .* s.Δt) # TODO: Use trapezoidal integration here (?),  this is just Forward Euler
-        mul!(Q(φ, s.B1 ./ B, Bz ./ B), M)
+        mul!(Q(φ, B1 ./ B, Bz ./ B), M)
         #Relaxation
         M.xy .= M.xy .* exp.(-s.Δt ./ p.T2)
         M.z .= M.z .* exp.(-s.Δt ./ p.T1) .+ p.ρ .* (1 .- exp.(-s.Δt ./ p.T1))

KomaMRIPlots/src/ui/DisplayFunctions.jl

@@ -992,8 +992,8 @@ Plots a phantom map for a specific spin parameter given by `key`.
 
 # Arguments
 - `obj`: (`::Phantom`) Phantom struct
-- `key`: (`::Symbol`, opts: [`:ρ`, `:T1`, `:T2`, `:T2s`, `:x`, `:y`, `:z`]) symbol for
-    displaying different parameters of the phantom spins
+- `key`: (`::Symbol`, opts: [`:ρ`, `:T1`, `:T2`, `:T2s`, `:x`, `:y`, `:z`, `:Δw`, `:B1`]) 
+    symbol for displaying different parameters of the phantom spins
 
 # Keywords
 - `height`: (`::Integer`, `=600`) plot height
@@ -1064,8 +1064,9 @@ function plot_phantom_map(
     end
 
     path = @__DIR__
-    cmin_key = minimum(getproperty(obj, key))
-    cmax_key = maximum(getproperty(obj, key))
+    this_map = getproperty(obj, key)
+    cmin_key = minimum(real.(this_map)) # allow for complex maps
+    cmax_key = maximum(real.(this_map))
     if key == :T1 || key == :T2 || key == :T2s
         cmin_key = 0
         factor = 1e3
@@ -1113,6 +1114,11 @@ function plot_phantom_map(
         factor = 1 / (2π)
         unit = " Hz"
         colormap = "Greys"
+    elseif key == :B1
+        factor = 1
+        this_map = real.(this_map)
+        unit = ""
+        colormap = "Greys"
     else
         factor = 1
         cmin_key = 0
@@ -1161,7 +1167,7 @@ function plot_phantom_map(
                 x=(x[:,i])*1e2,
                 y=(y[:,i])*1e2,
                 mode="markers",
-                marker=attr(color=getproperty(obj,key)*factor,
+                marker=attr(color=this_map*factor,
                             showscale=colorbar,
                             colorscale=colormap,
                             colorbar=attr(ticksuffix=unit, title=string(key)),
@@ -1171,7 +1177,7 @@ function plot_phantom_map(
                             ),
                 visible=i==1,
                 showlegend=false,
-                text=round.(getproperty(obj,key)*factor,digits=4),
+                text=round.(this_map*factor,digits=4),
                 hovertemplate="x: %{x:.1f} cm<br>y: %{y:.1f} cm<br><b>$(string(key))</b>: %{text}$unit<extra></extra>"))
         end
     else # 3D
@@ -1214,7 +1220,7 @@ function plot_phantom_map(
                 y=(y[:,i])*1e2,
                 z=(z[:,i])*1e2,
                 mode="markers",
-                marker=attr(color=getproperty(obj,key)*factor,
+                marker=attr(color=this_map*factor,
                             showscale=colorbar,
                             colorscale=colormap,
                             colorbar=attr(ticksuffix=unit, title=string(key)),
@@ -1224,7 +1230,7 @@ function plot_phantom_map(
                             ),
                 visible=i==1,
                 showlegend=false,
-                text=round.(getproperty(obj,key)*factor,digits=4),
+                text=round.(this_map*factor,digits=4),
                 hovertemplate="x: %{x:.1f} cm<br>y: %{y:.1f} cm<br>z: %{z:.1f} cm<br><b>$(string(key))</b>: %{text}$unit<extra></extra>"))
         end
     end

KomaMRIPlots/test/GUI_PlotlyJS_backend_test.jl

@@ -30,6 +30,11 @@
             @test true                #If the previous line fails the test will fail
         end
 
+        @testset "plot_phantom_map_B1" begin
+            plot_phantom_map(ph, :B1) #Plotting the phantom's rho map
+            @test true                #If the previous line fails the test will fail
+        end
+
         @testset "plot_phantom_map_2dview" begin
             plot_phantom_map(ph, :ρ, view_2d=true) #Plotting the phantom's rho map
             @test true                #If the previous line fails the test will fail
@@ -65,6 +70,11 @@
             @test true                #If the previous line fails the test will fail
         end
 
+        @testset "plot_motion_phantom_map_B1" begin
+            plot_phantom_map(ph, :B1) #Plotting the phantom's rho map
+            @test true                #If the previous line fails the test will fail
+        end
+
         @testset "plot_motion_phantom_map_2dview" begin
             plot_phantom_map(ph, :ρ, view_2d=true) #Plotting the phantom's rho map
             @test true                #If the previous line fails the test will fail

src/KomaUI.jl

@@ -40,7 +40,7 @@ function KomaUI(; darkmode=true, frame=true, phantom_mode="2D", sim=Dict{String,
     Observables.clear(img_ui)
 
     # For phantom sub-buttons
-    fieldnames_obj = [fieldnames(Phantom)[5:end-3]...]
+    fieldnames_obj = [fieldnames(Phantom)[[5:end-5;end]]...] # TODO should this be more explicit on which symbols to include rather than using indices?
     widgets_button_obj = button.(string.(fieldnames_obj))
 
     # Setup the Blink window

src/ui/ExportMATFunctions.jl

@@ -40,8 +40,8 @@ end
 
 function export_2_mat_phantom(phantom, matfolder; matfilename="phantom.mat")
     phantom_dict = Dict("name" => phantom.name,
-                "columns" => ["x", "y", "z", "rho", "T1", "T2", "T2s", "delta_omega"],
-                "data" => hcat(phantom.x, phantom.y, phantom.z, phantom.ρ, phantom.T1, phantom.T2, phantom.T2s, phantom.Δw))
+                "columns" => ["x", "y", "z", "rho", "T1", "T2", "T2s", "delta_omega", "B1"],
+                "data" => hcat(phantom.x, phantom.y, phantom.z, phantom.ρ, phantom.T1, phantom.T2, phantom.T2s, phantom.Δw, phantom.B1))
     matwrite(joinpath(matfolder, matfilename), Dict("phantom" => phantom_dict))
 end