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
@@ -289,7 +289,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
@@ -401,7 +401,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
@@ -483,7 +483,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,59 @@ 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 Base.setproperty!(sys::Scanner, key::Symbol, value)
+    if key in fieldnames(HardwareLimits)
+        setfield!(sys.limits, key, value)  # Convert value to the correct type
+    else
+        setfield!(sys, key, value)
+    end
+end
+
+function get_n_coils(rf_coils::RFCoils)
+    return 1
+end
+
+function get_n_coils(rf_coils::RFCoilsSensDefinedAtPhantomPositions)
+    return size(rf_coils.coil_sens, 2)
+end
+
+export ArbitraryRFCoils, RFCoilsSensDefinedAtPhantomPositions, UniformRFCoils, acquire_signal!, HardwareLimits, LinearXYZGradients, get_n_coils

KomaMRIBase/test/runtests.jl

@@ -345,7 +345,7 @@ end
     end
     @testset "Spiral" begin
         sys = Scanner()
-        sys.Smax = 150    # [mT/m/ms]
+        sys.Smax = 150.0    # [mT/m/ms]
         sys.Gmax = 500e-3 # [T/m]
         sys.GR_Δt = 4e-6  # [s]
         FOV = 0.2       # [m]
@@ -734,10 +734,23 @@ end
 end
 
 @testitem "Scanner" tags=[:base] begin
-    B0, B1, Gmax, Smax = 1.5, 10e-6, 60e-3, 500
+    B0, B1, Gmax, Smax = 1.5, 10e-6, 60e-3, 500.0
     ADC_Δt, seq_Δt, GR_Δt, RF_Δt = 2e-6, 1e-5, 1e-5, 1e-6
     RF_ring_down_T, RF_dead_time_T, ADC_dead_time_T = 20e-6, 100e-6, 10e-6
-    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)
+    limits = HardwareLimits(
+        B0 = B0, 
+        B1 = B1, 
+        Gmax = Gmax, 
+        Smax = Smax, 
+        ADC_Δt = ADC_Δt, 
+        seq_Δt = seq_Δt, 
+        GR_Δt = GR_Δt, 
+        RF_Δt = RF_Δt, 
+        RF_ring_down_T = RF_ring_down_T, 
+        RF_dead_time_T = RF_dead_time_T, 
+        ADC_dead_time_T = ADC_dead_time_T
+    )
+    sys = Scanner(limits = limits)
     @test sys.B0 ≈ B0 && sys.B1 ≈ B1 && sys.Gmax ≈ Gmax && sys.Smax ≈ Smax
 end
 

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,
     sig::AbstractArray{Complex{T}},
     M::Mag{T},
     sim_method::Bloch,
@@ -66,11 +67,10 @@ function run_spin_precession!(
     ΔBz = prealloc.ΔBz
     fill!(ϕ, zero(T))
     @. Bz_old = x[:,1] * seq.Gx[1] + y[:,1] * seq.Gy[1] + z[:,1] * seq.Gz[1] + ΔBz
-
     # Fill sig[1] if needed
     ADC_idx = 1
     if (seq.ADC[1])
-        sig[1] = sum(M.xy)
+        sig[1] = sum(sys.rf_coils.coil_sens .* M.xy)
         ADC_idx += 1
     end
 
@@ -91,8 +91,7 @@ function run_spin_precession!(
 
             #Reset Spin-State (Magnetization). Only for FlowPath
             outflow_spin_reset!(Mxy, seq.t[seq_idx], p.motion)
-
-            sig[ADC_idx] = sum(Mxy) 
+            acquire_signal!(transpose(@view(sig[ADC_idx,:])), sys.rf_coils, Mxy)
             ADC_idx += 1
         end
 

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

@@ -125,6 +125,7 @@ end
 function run_spin_precession!(
     p::Phantom{T},
     seq::DiscreteSequence{T},
+    sys::Scanner,
     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,
     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"],

KomaMRICore/test/runtests.jl

@@ -370,7 +370,7 @@ end
 
     # This is a sequence with a sinc RF 30° excitation pulse
     sys = Scanner()
-    sys.Smax = 50
+    sys.Smax = 50.0
     B1 = 4.92e-6
     Trf = 3.2e-3
     zmax = 2e-2

benchmarks/setup.jl

@@ -7,7 +7,7 @@ end
 
 function setup_MRILab_benchmark()
     sys = Scanner()
-    sys.Smax = 150    # [mT/m/ms]
+    sys.Smax = 150.0    # [mT/m/ms]
     sys.Gmax = 500e-3 # [T/m]
     sys.GR_Δt = 4e-6  # [s]
     FOV = 0.2         # [m]

examples/3.tutorials/lit-02-SmallTipApproximation.jl

@@ -4,7 +4,7 @@
 
 using KomaMRI # hide
 sys = Scanner() # hide
-sys.Smax = 50 # hide
+sys.Smax = 50.0 # hide
 
 # In this example, we will showcase a common approximation in MRI, the small tip angle approximation.
 # For this, we will simulate a slice profile for spins with positions ``z\in[-2,\,2]\,\mathrm{cm}``