Scanner object considering multiple coils

KomaMRIBase/src/datatypes/Phantom.jl

@@ -34,7 +34,7 @@ julia> obj.ρ
     x::AbstractVector{T}   = @isdefined(T) ? T[] : Float64[]
     y::AbstractVector{T}   = zeros(eltype(x), size(x))
     z::AbstractVector{T}   = zeros(eltype(x), size(x))
-    ρ::AbstractVector{T}   = ones(eltype(x), size(x))
+    ρ::Union{AbstractVector{T}, AbstractVector{Complex{T}}}   = ones(eltype(x), size(x))
     T1::AbstractVector{T}  = ones(eltype(x), size(x)) * 1_000_000
     T2::AbstractVector{T}  = ones(eltype(x), size(x)) * 1_000_000
     T2s::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000
@@ -124,7 +124,7 @@ end
         heart_rate, asymmetry
     )
 
-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.
 
 # Keywords
@@ -227,7 +227,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
@@ -323,7 +323,7 @@ end
     obj = brain_phantom3D(; ss=4, us=1, start_end=[160,200])
 
 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
@@ -350,7 +350,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 +473,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 +556,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/src/datatypes/Scanner.jl

@@ -1,3 +1,44 @@
+# Hardware limits
+@with_kw mutable struct HardwareLimits{T}
+    B0::T = 1.5
+    B1::T = 10e-6
+    Gmax::T = 60e-3
+    Smax::T = 500.0
+    ADC_Δt::T = 2e-6
+    seq_Δt::T = 1e-5
+    GR_Δt::T = 1e-5
+    RF_Δt::T = 1e-6
+    RF_ring_down_T::T = 20e-6
+    RF_dead_time_T::T = 100e-6
+    ADC_dead_time_T::T = 10e-6
+end
+
+# Gradients
+abstract type Gradients{T} end
+@with_kw mutable struct LinearXYZGradients{T} <: Gradients{T} 
+    Gx::AbstractVector{T} = zeros(T, 1)
+    Gy::AbstractVector{T} = zeros(T, 1)
+    Gz::AbstractVector{T} = zeros(T, 1)
+end
+
+# RF coils
+abstract type RFCoils{T} end
+@with_kw mutable struct UniformRFCoils{T} <: RFCoils{T} 
+    coil_sens::AbstractMatrix{Complex{T}} = complex(ones(Complex{T}, 1, 1))
+end
+
+struct ArbitraryRFCoils{T} <: RFCoils{T}
+    x::AbstractVector{T} 
+    y::AbstractVector{T} 
+    z::AbstractVector{T} 
+    coil_sens::AbstractMatrix{Complex{T}}  
+    B1⁺::AbstractMatrix{Complex{T}} 
+end
+
+struct RFCoilsSensDefinedAtPhantomPositions{T} <: RFCoils{T}
+    coil_sens::AbstractMatrix{Complex{T}}
+end
+
 """
     sys = Scanner(B0, B1, Gmax, Smax, ADC_Δt, seq_Δt, GR_Δt, RF_Δt,
         RF_ring_down_T, RF_dead_time_T, ADC_dead_time_T)
@@ -28,19 +69,51 @@ julia> sys = Scanner()
 julia> sys.B0
 ```
 """
-@with_kw mutable struct Scanner
-    #Main
-    B0::Real=1.5
-    B1::Real=10e-6
-    Gmax::Real=60e-3
-    Smax::Real=500
-    #Sampling
-    ADC_Δt::Real=2e-6
-    seq_Δt::Real=1e-5
-    GR_Δt::Real=1e-5
-    RF_Δt::Real=1e-6
-    #Secondary
-    RF_ring_down_T::Real=20e-6
-    RF_dead_time_T::Real=100e-6
-    ADC_dead_time_T::Real=10e-6
+@with_kw mutable struct Scanner{T}
+    limits::HardwareLimits{T} = HardwareLimits{Float64}()
+    gradients::Gradients{T} = LinearXYZGradients{Float64}()
+    rf_coils::RFCoils{T} = UniformRFCoils{Float64}()
+end
+
+function Base.view(sys::Scanner, p)
+    return Scanner(limits=sys.limits, gradients=sys.gradients, rf_coils=@view(sys.rf_coils[p]))
+end
+
+function Base.view(rf_coils::RFCoilsSensDefinedAtPhantomPositions, p)
+    return RFCoilsSensDefinedAtPhantomPositions(@view rf_coils.coil_sens[p,:])
 end
+
+function Base.view(rf_coils::UniformRFCoils, p)
+    return rf_coils
+end
+
+function acquire_signal!(sig, rf_coils::UniformRFCoils, Mxy)
+    sig .= transpose(sum(Mxy; dims=1))
+    return nothing
+end
+
+
+function acquire_signal!(sig, rf_coils::RFCoilsSensDefinedAtPhantomPositions, Mxy)
+    for i in 1:size(rf_coils.coil_sens, 2)
+        sig[:, i] .= transpose(sum(rf_coils.coil_sens[:, i] .* Mxy; dims=1))
+    end
+    return nothing
+end
+
+function getproperty(sys::Scanner, key::Symbol)
+    if key in fieldnames(HardwareLimits)
+        return getfield(sys.limits, key)
+    else
+        return getfield(sys, key)
+    end
+end
+
+function get_n_coils(rf_coils::UniformRFCoils)
+    return 1
+end
+
+function get_n_coils(rf_coils::RFCoils)
+    return size(rf_coils.coil_sens, 2)
+end
+
+export ArbitraryRFCoils, RFCoilsSensDefinedAtPhantomPositions, UniformRFCoils, acquire_signal!, HardwareLimits, LinearXYZGradients, get_n_coils

KomaMRICore/src/simulation/Functors.jl

@@ -105,6 +105,7 @@ adapt_storage(T::Type{<:Real}, xs::MotionList) = MotionList(paramtype.(T, xs.mot
 # Phantom
 @functor Phantom
 @functor Motion
+@functor Scanner
 @functor Translate
 @functor Rotate
 @functor HeartBeat

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

@@ -48,6 +48,7 @@ that they can be re-used from block to block.
 function run_spin_precession!(
     p::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner{T},
     sig::AbstractArray{Complex{T}},
     M::Mag{T},
     sim_method::Bloch,

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

@@ -125,6 +125,7 @@ end
 function run_spin_precession!(
     p::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner{T},
     sig::AbstractArray{Complex{T}},
     M::Mag{T},
     sim_method::Bloch,
@@ -159,7 +160,7 @@ function run_spin_precession!(
             outflow_spin_reset!(pre.Mxy, seq_block.tp_ADC', p.motion; seq_t=seq.t)
         end
 
-        sig .= transpose(sum(pre.Mxy; dims=1))
+        acquire_signal!(sig, sys.rf_coils, pre.Mxy)
     end
 
     #Mxy precession and relaxation, and Mz relaxation

KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl

@@ -32,6 +32,7 @@ precession.
 function run_spin_precession!(
     p::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner{T},
     sig::AbstractArray{Complex{T}},
     M::Mag{T},
     sim_method::BlochDict,

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

@@ -24,6 +24,7 @@ precession.
 function run_spin_precession!(
     p::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner,
     sig::AbstractArray{Complex{T}},
     M::Mag{T},
     sim_method::SimulationMethod,
@@ -51,7 +52,7 @@ function run_spin_precession!(
     outflow_spin_reset!(Mxy, seq.t', p.motion)
     outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ)
     #Acquired signal
-    sig .= transpose(sum(Mxy[:, findall(seq.ADC)]; dims=1)) #<--- TODO: add coil sensitivities
+    acquire_signal!(sig, sys.rf_coils, Mxy[:, findall(seq.ADC)])
     return nothing
 end
 

KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl

@@ -10,7 +10,9 @@ include("Magnetization.jl")
 function sim_output_dim(
     obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::SimulationMethod
 ) where {T<:Real}
-    return (sum(seq.ADC.N), 1) #Nt x Ncoils, This should consider the coil info from sys
+    # Determine the number of coils
+    n_coils = get_n_coils(sys.rf_coils)
+    return (sum(seq.ADC.N), n_coils) # Nt x Ncoils
 end
 
 """Magnetization initialization for Bloch simulation method."""

KomaMRICore/src/simulation/SimulatorCore.jl

@@ -79,6 +79,7 @@ separating the spins of the phantom `obj` in `Nthreads`.
 function run_spin_precession_parallel!(
     obj::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner,
     sig::AbstractArray{Complex{T}},
     Xt::SpinStateRepresentation{T},
     sim_method::SimulationMethod,
@@ -91,7 +92,7 @@ function run_spin_precession_parallel!(
 
     ThreadsX.foreach(enumerate(parts)) do (i, p)
         run_spin_precession!(
-            @view(obj[p]), seq, @view(sig[dims..., i]), @view(Xt[p]), sim_method, backend, @view(prealloc[p])
+            @view(obj[p]), seq, @view(sys[p]), @view(sig[dims..., i]), @view(Xt[p]), sim_method, backend, @view(prealloc[p])
         )
     end
 
@@ -167,6 +168,7 @@ take advantage of CPU parallel processing.
 function run_sim_time_iter!(
     obj::Phantom,
     seq::DiscreteSequence,
+    sys::Scanner,
     sig::AbstractArray{Complex{T}},
     Xt::SpinStateRepresentation{T},
     sim_method::SimulationMethod,
@@ -199,7 +201,7 @@ function run_sim_time_iter!(
             rfs += 1
         else
             run_spin_precession_parallel!(
-                obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method, backend, prealloc_block(prealloc_result, block); Nthreads
+                obj, seq_block, sys, @view(sig[acq_samples, dims...]), Xt, sim_method, backend, prealloc_block(prealloc_result, block); Nthreads
             )
         end
         samples += Nadc
@@ -378,6 +380,7 @@ function simulate(
         seqd = seqd |> gpu #DiscreteSequence
         Xt = Xt |> gpu #SpinStateRepresentation
         sig = sig |> gpu #Signal
+        sys = sys |> gpu #Scanner
         precalc = precalc |> gpu #Info calculated prior to simulation
     end
 
@@ -389,6 +392,7 @@ function simulate(
     @time timed_tuple = @timed run_sim_time_iter!(
         obj,
         seqd,
+        sys,
         sig,
         Xt,
         sim_params["sim_method"],

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

@@ -0,0 +1,57 @@
+# # Experimental: Simulating with realistic coils
+
+using KomaMRI # hide
+obj = brain_phantom2D()
+coil_sens1 = exp.(-π * (((obj.x) .+ 0.1) .^ 2 / 0.02) .+ ((obj.y) .^ 2 / 0.02))
+coil_sens2 = exp.(-π * (((obj.x) .- 0.1) .^ 2 / 0.02) .+ ((obj.y) .^ 2 / 0.02))
+coil_sens3 = exp.(-π * ((obj.x) .^ 2 / 0.02) .+ (((obj.y) .+ 0.1) .^ 2 / 0.02))
+coil_sens4 = exp.(-π * ((obj.x) .^ 2 / 0.02) .+ (((obj.y) .- 0.1) .^ 2 / 0.02))
+coil_sens = hcat(coil_sens1, coil_sens2, coil_sens3, coil_sens4)
+sys = Scanner()
+sys.rf_coils = RFCoilsSensDefinedAtPhantomPositions(complex.(coil_sens))
+#sys.rf_coils = ArbitraryRFCoils(obj.x, obj.y, obj.z, complex.(coil_sens), complex.(coil_sens))
+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:
+
+sim_params = KomaMRICore.default_sim_params()
+sim_params["sim_method"] = Bloch()
+raw = simulate(obj, seq, sys; sim_params)
+
+acq = AcquisitionData(raw) # hide
+acq.traj[1].circular = false # hide
+Nx, Ny = raw.params["reconSize"][1:2] # hide
+# reco parameters
+params = Dict{Symbol,Any}()
+params[:reco] = "multiCoil"
+params[:reconSize] = (Nx, Ny)
+params[:iterations] = 40
+
+# Coil sensitivities interpolated to reconstruction points
+FOV = 230e-3
+xq = range(-FOV / 2, FOV / 2, Nx)
+yq = range(-FOV / 2, FOV / 2, Ny)
+coil_sens1 = exp.(-π * (((xq) .+ 0.1) .^ 2 / 0.02) .+ ((yq') .^ 2 / 0.02))
+coil_sens2 = exp.(-π * (((xq) .- 0.1) .^ 2 / 0.02) .+ ((yq') .^ 2 / 0.02))
+coil_sens3 = exp.(-π * ((xq) .^ 2 / 0.1) .+ (((yq') .+ 0.1) .^ 2 / 0.02))
+coil_sens4 = exp.(-π * ((xq) .^ 2 / 0.1) .+ (((yq') .- 0.1) .^ 2 / 0.02))
+coil_sens_recon = Float32.([coil_sens1[:] coil_sens2[:] coil_sens3[:] coil_sens4[:]])
+params[:senseMaps] = reshape(complex.(coil_sens_recon), Nx, Ny, 1, size(coil_sens_recon, 2))
+
+# do reconstruction
+Ireco = reconstruction(acq, params)
+recon_image = plot_image(abs.(Ireco.data[:, :]))
+reconParams = Dict{Symbol,Any}(:reco => "direct", :reconSize => (Nx, Ny)) # hide
+image = reconstruction(acq, reconParams) # hide
+slice_abs1 = abs.(image[:, :, 1, 1, 1, 1]) # hide
+slice_abs2 = abs.(image[:, :, 1, 1, 2, 1]) # hide
+slice_abs3 = abs.(image[:, :, 1, 1, 3, 1]) # hide
+slice_abs4 = abs.(image[:, :, 1, 1, 4, 1]) # hide
+p3 = plot_image(slice_abs1; height=400)
+p4 = plot_image(slice_abs2; height=400) # hide
+p5 = plot_image(slice_abs3; height=400)
+p6 = plot_image(slice_abs4; height=400) # hide
+[recon_image; p3 p4 p5 p6]