Trapezoidal rule for RFs

KomaMRIBase/src/datatypes/simulation/DiscreteSequence.jl

@@ -27,8 +27,6 @@ struct DiscreteSequence{T<:Real}
     t::AbstractVector{T}
     Δt::AbstractVector{T}
 end
-
-Base.length(seq::DiscreteSequence) = length(seq.Δt)
 Base.getindex(seq::DiscreteSequence, i::Integer) = begin
     DiscreteSequence(seq.Gx[i, :],
                      seq.Gy[i, :],
@@ -39,35 +37,34 @@ Base.getindex(seq::DiscreteSequence, i::Integer) = begin
                      seq.t[i, :],
                      seq.Δt[i, :])
 end
+# Use in precession trapezoidal rule
 Base.getindex(seq::DiscreteSequence, i::UnitRange) = begin
-    DiscreteSequence(seq.Gx[i.start:i.stop+1],
-                     seq.Gy[i.start:i.stop+1],
-                     seq.Gz[i.start:i.stop+1],
-                     seq.B1[i.start:i.stop+1],
-                     seq.Δf[i.start:i.stop+1],
+    DiscreteSequence(seq.Gx[i],
+                     seq.Gy[i],
+                     seq.Gz[i],
+                     seq.B1[i],
+                     seq.Δf[i],
                      seq.ADC[i],
-                     seq.t[i.start:i.stop+1],
-                     seq.Δt[i])
+                     seq.t[i],
+                     seq.Δt[i.start:i.stop-1])
 end
 Base.view(seq::DiscreteSequence, i::UnitRange) = begin
-    @views DiscreteSequence(seq.Gx[i.start:i.stop+1],
-                     seq.Gy[i.start:i.stop+1],
-                     seq.Gz[i.start:i.stop+1],
-                     seq.B1[i.start:i.stop+1],
-                     seq.Δf[i.start:i.stop+1],
+    @views DiscreteSequence(seq.Gx[i],
+                     seq.Gy[i],
+                     seq.Gz[i],
+                     seq.B1[i],
+                     seq.Δf[i],
                      seq.ADC[i],
-                     seq.t[i.start:i.stop+1],
-                     seq.Δt[i])
+                     seq.t[i],
+                     seq.Δt[i.start:i.stop-1])
 end
+Base.length(seq::DiscreteSequence)  = length(seq.Δt)
 Base.iterate(seq::DiscreteSequence) = (seq[1], 2)
 Base.iterate(seq::DiscreteSequence, i) = (i <= length(seq)) ? (seq[i], i+1) : nothing
 
-is_GR_on(seq::DiscreteSequence) =  sum(abs.([seq.Gx[1:end-1]; seq.Gy[1:end-1]; seq.Gz[1:end-1]])) != 0
-is_RF_on(seq::DiscreteSequence) =  sum(abs.(seq.B1[1:end-1])) != 0
-is_ADC_on(seq::DiscreteSequence) = sum(abs.(seq.ADC[1:end-1])) != 0
-is_GR_off(seq::DiscreteSequence) =  !is_GR_on(seq)
-is_RF_off(seq::DiscreteSequence) =  !is_RF_on(seq)
-is_ADC_off(seq::DiscreteSequence) = !is_ADC_on(seq)
+is_GR_on(seq::DiscreteSequence)   = sum(abs.([seq.Gx; seq.Gy; seq.Gz])) != 0
+is_RF_on(seq::DiscreteSequence)   = sum(abs.(seq.B1)) != 0
+is_ADC_on(seq::DiscreteSequence)  = sum(abs.(seq.ADC)) != 0
 
 """
     seqd = discretize(seq::Sequence; sampling_params=default_sampling_params())

KomaMRICore/src/datatypes/VectorSU2.jl

@@ -0,0 +1,31 @@
+mutable struct VectorSU2{T}
+    xy::AbstractVector{Complex{T}}
+    z::AbstractVector{T}
+end
+
+# Required indexing operations
+Base.view(V::VectorSU2, i::UnitRange) = @views VectorSU2(V.xy[i], V.z[i])
+
+# Maths
+Base.abs(V::VectorSU2) = sqrt.(abs.(V.xy).^2 + (V.z).^2)
+Base.:/(V::VectorSU2, α) = VectorSU2(V.xy/α, V.z/α)
+Base.:+(V1::VectorSU2, V2::VectorSU2) = VectorSU2(V1.xy .+ V2.xy, V1.z .+ V2.z)
+
+# Definition of Spinor from VectorSU2
+function calculateRot!(pre::PreallocResult{T}, Δt::T, sim_method::Bloch) where {T}
+    if sim_method.rf_ode_order == 2
+        # 1st order
+        Beff = (pre.B_old + pre.B_new) / 2
+        absBeff = abs(Beff)
+    elseif sim_method.rf_ode_order == 1
+        # 0th order
+        Beff = pre.B_old
+        absBeff = abs(Beff)
+        absBeff[absBeff .== zero(T)] .= eps(T) # not sure why this is needed, in RF block this should be non-zero
+    end
+    # Rot
+    @. pre.φ = T(-2π * γ) * absBeff * Δt 
+    @. pre.Rot.α = cos(pre.φ/2) - Complex{T}(im) * (Beff.z / absBeff) * sin(pre.φ/2)
+    @. pre.Rot.β = -Complex{T}(im) * (Beff.xy / absBeff) * sin(pre.φ/2)
+    return nothing
+end

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

@@ -1,39 +1,27 @@
 """Stores preallocated structs for use in Bloch CPU run_spin_precession function."""
-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)
-    ϕ::AbstractVector{T}                    # Vector{T}(Nspins x 1)
-    φ::AbstractVector{T}                    # Vector{T}(Nspins x 1)
-    Rot::Spinor{T}                          # Spinor{T}
+mutable struct BlochCPUPrealloc{T} <: PreallocResult{T}
+    B_old::VectorSU2{T}         # Vector{T}(Nspins x 1)
+    B_new::VectorSU2{T}         # Vector{T}(Nspins x 1)
+    φ::AbstractVector{T}        # Vector{T}(Nspins x 1)
+    Rot::Spinor{T}              # Spinor{T}
 end
 
 Base.view(p::BlochCPUPrealloc, i::UnitRange) = begin
-    @views BlochCPUPrealloc(
-        p.M[i],
-        p.Bz_old[i],
-        p.Bz_new[i],
-        p.ϕ[i],
-        p.φ[i],
-        p.Rot[i]
+    BlochCPUPrealloc(
+        @view(p.B_old[i]),
+        @view(p.B_new[i]),
+        @view(p.φ[i]),
+        @view(p.Rot[i])
     )
 end
 
 """Preallocates arrays for use in run_spin_precession."""
 function prealloc(sim_method::Bloch, backend::KA.CPU, obj::Phantom{T}, M::Mag{T}) where {T<:Real}
     return BlochCPUPrealloc(
-        Mag(
-            similar(M.xy),
-            similar(M.z)
-        ),
-        similar(obj.x),
-        similar(obj.x),
-        similar(obj.x),
-        similar(obj.x),
-        Spinor(
-            similar(M.xy),
-            similar(M.xy)
-        )
+        VectorSU2(zero(M.xy), zero(M.z)), #TODO: change this to not use M.xy
+        VectorSU2(zero(M.xy), zero(M.z)),
+        zero(obj.x),
+        Spinor(zero(M.xy), zero(M.xy))
     )
 end
 
@@ -52,52 +40,39 @@ function run_spin_precession!(
     M::Mag{T},
     sim_method::Bloch,
     backend::KA.CPU,
-    prealloc::BlochCPUPrealloc
+    pre::BlochCPUPrealloc
 ) where {T<:Real}
-    #Simulation
-    #Motion
-    x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[1,:]')
-
-    #Initialize arrays
-    Bz_old = prealloc.Bz_old
-    Bz_new = prealloc.Bz_new
-    ϕ = prealloc.ϕ
-    Mxy = prealloc.M.xy
-    fill!(ϕ, zero(T))
-    @. Bz_old = x[:,1] * seq.Gx[1] + y[:,1] * seq.Gy[1] + z[:,1] * seq.Gz[1] + p.Δw / T(2π * γ)
-
-    # Fill sig[1] if needed
+    # Initialize
     ADC_idx = 1
-    if (seq.ADC[1])
-        sig[1] = sum(M.xy)
-        ADC_idx += 1
-    end
-
-    t_seq = zero(T) # Time
-    for seq_idx=2:length(seq.t)
-        x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[seq_idx,:]')
-        t_seq += seq.Δt[seq_idx-1]
-
-        #Effective Field
-        @. Bz_new = x * seq.Gx[seq_idx] + y * seq.Gy[seq_idx] + z * seq.Gz[seq_idx] + p.Δw / T(2π * γ)
-
-        #Rotation
-        @. ϕ += (Bz_old + Bz_new) * T(-π * γ) * seq.Δt[seq_idx-1]
-
-        #Acquired Signal
-        if seq_idx <= length(seq.ADC) && seq.ADC[seq_idx]
-            @. Mxy = exp(-t_seq / p.T2) * M.xy * cis(ϕ)
-            sig[ADC_idx] = sum(Mxy) 
+    Δt = zero(T) # Time
+    pre.φ .= zero(T)
+    # Simulation    
+    for i in 1:length(seq.Δt)
+        # Bz_eff(tn+1)
+        x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i+1,:]')
+        @. pre.B_new.z = x * seq.Gx[i+1] + y * seq.Gy[i+1] + z * seq.Gz[i+1] + p.Δw / T(2π * γ)
+        # Rotation angle
+        @. pre.φ += T(-2π * γ) * (pre.B_old.z + pre.B_new.z) * seq.Δt[i] / 2
+        Δt += seq.Δt[i]
+        # Acquire signal
+        if seq.ADC[i+1] # TODO: || sim_method.sample_all
+            # Decay + Rot
+            @. M.xy *= exp(-Δt / p.T2) * cis(pre.φ)
+            # Sample
+            sig[ADC_idx] = sum(M.xy) 
             ADC_idx += 1
+            # Reset, decay and rot operators
+            Δt = zero(T)
+            pre.φ .= zero(T)
         end
-
-        Bz_old, Bz_new = Bz_new, Bz_old
+        # Update
+        pre.B_old.z .= pre.B_new.z
     end
-
-    #Final Spin-State
-    @. M.xy = M.xy * exp(-t_seq / p.T2) * cis(ϕ)
-    @. M.z = M.z * exp(-t_seq / p.T1) + p.ρ * (T(1) - exp(-t_seq / p.T1))
-
+    pre.B_old.xy .= zero(Complex{T})
+    # Final Spin-State
+    t_total = sum(seq.Δt)
+    @. M.xy *= exp(-Δt / p.T2) * cis(pre.φ)
+    @. M.z = M.z * exp(-t_total / p.T1) + p.ρ * (1 - exp(-t_total / p.T1))
     return nothing
 end
 
@@ -114,38 +89,31 @@ function run_spin_excitation!(
     M::Mag{T},
     sim_method::Bloch,
     backend::KA.CPU,
-    prealloc::BlochCPUPrealloc
+    pre::BlochCPUPrealloc
 ) where {T<:Real}
-    ΔBz = prealloc.Bz_old
-    Bz = prealloc.Bz_new
-    B = prealloc.ϕ
-    φ = prealloc.φ
-    α = prealloc.Rot.α
-    β = prealloc.Rot.β
-    Maux_xy = prealloc.M.xy
-    Maux_z = prealloc.M.z
-
-    #Can be calculated outside of loop
-    @. ΔBz = p.Δw / T(2π * γ)
-
-    #Simulation
-    for s in seq #This iterates over seq, "s = seq[i,:]"
-        #Motion
-        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)
-        @. 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(φ)
-        mul!(Spinor(α, β), M, Maux_xy, Maux_z)
-        #Relaxation
-        @. M.xy = M.xy * exp(-s.Δt / p.T2)
-        @. M.z = M.z * exp(-s.Δt / p.T1) + p.ρ * (T(1) - exp(-s.Δt / p.T1))
+    # Init
+    sample = 1
+    # Simulation
+    for i in 1:length(seq.Δt)
+        # Motion
+        x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i+1, :])
+        # Effective field
+        @. pre.B_new.xy = seq.B1[i+1]
+        @. pre.B_new.z  = (seq.Gx[i+1] * x + seq.Gy[i+1] * y + seq.Gz[i+1] * z) + p.Δw / T(2π * γ) - seq.Δf[i+1] / T(γ)
+        # Rotation
+        calculateRot!(pre, seq.Δt[i], sim_method)
+        mul!(pre.Rot, M)
+        # Relaxation
+        @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2)
+        @. M.z  = M.z * exp(-seq.Δt[i] / p.T1) + p.ρ * (1 - exp(-seq.Δt[i] / p.T1))
+        # Sample
+        if seq.ADC[i+1]
+            sig[sample, :] .= sum(M.xy)
+            sample += 1
+        end
+        # Update
+        pre.B_old.xy .= pre.B_new.xy
+        pre.B_old.z  .= pre.B_new.z
     end
-    #Acquired signal
-    #sig .= -1.4im #<-- This was to test if an ADC point was inside an RF block
     return nothing
 end

KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl

@@ -1,8 +1,10 @@
 #Bloch is the default simulation method
-struct Bloch <: SimulationMethod end
+struct Bloch <: SimulationMethod
+    rf_ode_order::Int
+end
+Bloch() = Bloch(1)
 
 export Bloch
-
 include("Magnetization.jl")
 
 @functor Mag #Gives gpu acceleration capabilities, see GPUFunctions.jl
@@ -34,6 +36,7 @@ Base.view(p::DefaultPreAlloc, i::UnitRange) = p
 """Default preallocation function."""
 prealloc(sim_method::SimulationMethod, backend::KA.Backend, obj::Phantom{T}, M::Mag{T}) where {T<:Real} = DefaultPreAlloc{T}()
 
+include("../../datatypes/VectorSU2.jl")
 include("KernelFunctions.jl")
 include("BlochSimple/BlochSimple.jl")
 include("Bloch/BlochCPU.jl")

KomaMRICore/src/simulation/SimulatorCore.jl

@@ -187,7 +187,7 @@ function run_sim_time_iter!(
         seq_block = @view seq[p]
         # Params
         # excitation_bool = is_RF_on(seq_block) #&& is_ADC_off(seq_block) #PATCH: the ADC part should not be necessary, but sometimes 1 sample is identified as RF in an ADC block
-        Nadc = sum(seq_block.ADC)
+        Nadc = sum(seq_block.ADC[2:end]) # The first sample is considered in the previous block
         acq_samples = samples:(samples + Nadc - 1)
         dims = [Colon() for i in 1:(ndims(sig) - 1)] # :,:,:,... Ndim times
         # Simulation wrappers
@@ -230,7 +230,7 @@ function get_sim_ranges(seqd::DiscreteSequence; Nblocks)
     start_idx_rf_block = 0
     start_idx_gr_block = 0
     #Split 1:N into Nblocks like kfoldperm
-    N = length(seqd.Δt)
+    N = length(seqd.t)
     k = min(N, Nblocks)
     n, r = divrem(N, k) #N >= k, N < k
     breaks = collect(1:n:(N + 1))
@@ -239,43 +239,43 @@ function get_sim_ranges(seqd::DiscreteSequence; Nblocks)
     end
     breaks = breaks[2:(end - 1)] #Remove borders,
     #Iterate over B1 values to decide the simulation UnitRanges
-    for i in eachindex(seqd.Δt)
+    for i in eachindex(seqd.t)
         if abs(seqd.B1[i]) > 1e-9 #TODO: This is needed as the function ⏢ in get_rfs is not very accurate
             if start_idx_rf_block == 0 #End RF block
-                start_idx_rf_block = i
+                start_idx_rf_block = i - 1
             end
             if start_idx_gr_block > 0 #End of GR block
-                push!(ranges, start_idx_gr_block:(i - 1))
+                push!(ranges, start_idx_gr_block:i-1)
                 push!(ranges_bool, false)
                 start_idx_gr_block = 0
             end
         else
             if start_idx_gr_block == 0 #Start GR block
-                start_idx_gr_block = i
+                start_idx_gr_block = i - 1
             end
             if start_idx_rf_block > 0 #End of RF block
-                push!(ranges, start_idx_rf_block:(i - 1))
+                push!(ranges, start_idx_rf_block:i-1)
                 push!(ranges_bool, true)
                 start_idx_rf_block = 0
             end
         end
         #More subdivisions
-        if i in breaks
-            if start_idx_rf_block > 0 #End of RF block
-                if length(start_idx_rf_block:(i - 1)) > 1
-                    push!(ranges, start_idx_rf_block:(i - 1))
-                    push!(ranges_bool, true)
-                    start_idx_rf_block = i
-                end
-            end
-            if start_idx_gr_block > 0 #End of RF block
-                if length(start_idx_gr_block:(i - 1)) > 1
-                    push!(ranges, start_idx_gr_block:(i - 1))
-                    push!(ranges_bool, false)
-                    start_idx_gr_block = i
-                end
-            end
-        end
+        # if i in breaks
+        #     if start_idx_rf_block > 0 #End of RF block
+        #         if length(start_idx_rf_block:(i - 1)) > 1
+        #             push!(ranges, start_idx_rf_block:(i - 1))
+        #             push!(ranges_bool, true)
+        #             start_idx_rf_block = i - 1
+        #         end
+        #     end
+        #     if start_idx_gr_block > 0 #End of RF block
+        #         if length(start_idx_gr_block:(i - 1)) > 1
+        #             push!(ranges, start_idx_gr_block:(i - 1))
+        #             push!(ranges_bool, false)
+        #             start_idx_gr_block = i
+        #         end
+        #     end
+        # end
     end
     #Finishing the UnitRange's
     if start_idx_rf_block > 0