Open cncastillo opened 1 year ago
I think this is related to the goal of adding automatic differentiation.
I already started doing this for the Grads in #321
using Parameters
abstract type MRISequenceEvent end
NumOrVect = Union{T, Vector{T}} where {T<:Number}
@with_kw_noshow mutable struct Grad{GradAmp<:NumOrVect, TimeDur<:NumOrVect} <: MRISequenceEvent
# Grad properties
A ::GradAmp = zero(eltype(GradAmp))
T ::TimeDur = zero(eltype(TimeDur))
rise ::eltype(TimeDur) = zero(eltype(T))
fall ::eltype(TimeDur) = rise
delay::eltype(TimeDur) = zero(eltype(T))
first::eltype(GradAmp) = zero(eltype(A))
last ::eltype(GradAmp) = first
# Checking validity of Grad
@assert length(A) >= length(T)
@assert all(T .>= zero(eltype(T)))
@assert rise >= zero(eltype(T))
@assert fall >= zero(eltype(T))
@assert delay >= zero(eltype(T))
end
# Gradient type aliases for dispatching
TrapezoidalGrad = Grad{G, T} where {G, T}
UniformlySampledGrad = Grad{Vector{G}, T} where {G, T}
TimeShapedGrad = Grad{Vector{G}, Vector{T}} where {G, T}
ValidGrad = Union{TrapezoidalGrad{G, T}, UniformlySampledGrad{G, T}, TimeShapedGrad{G, T}} where {G, T}
# MRISequenceEvent functions
# Event duration
dur(g::Grad) = g.delay + g.rise + sum(g.T) + g.fall
dur(G::Vector{Grad{G, T}}) where {G, T} = max(dur.(G)...)
# Event amplitude
ampl(g::TrapezoidalGrad{G,T}) where {G, T} = vcat(zero(G), g.first, g.A, g.A, g.last)
ampl(g::TimeShapedGrad{G,T}) where {G, T} = vcat(zero(G), g.first, g.A, g.last)
ampl(g::UniformlySampledGrad{G,T}) where {G, T} = vcat(zero(G), g.first, g.A, g.last)
# Event time
time(g::TrapezoidalGrad{G,T}) where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::TimeShapedGrad{G,T}) where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::UniformlySampledGrad{G,T}) where {G, T} = begin
NA = length(g.A) - 1
ΔT = repeat([g.T / NA], NA)
return cumsum(vcat(zero(T), g.delay, g.rise, ΔT, g.fall))
end
For Unitful support
using Unitful
using Unitful: 𝐈, 𝐌, 𝐓, 𝐋
using Unitful: T, m, mT, s, ms, Hz, MHz
Base.show(io::IO, g::ValidGrad{QA, QT}) where {QA <: Quantity, QT <: Quantity} = begin
printfield(g, f) = begin
val = getfield(g, f)
if eltype(val) <: QA && !all(abs.(val) .== zero(QA))
return "$f=$(ustrip(val)) $(unit(QA))"
elseif eltype(val) <: QT && !all(abs.(val) .== zero(QT))
return "$f=$(ustrip(val)) $(unit(QT))"
else
return ""
end
end
text = [printfield(g, f) for f in fieldnames(Grad)]
text = join(text[text.!=""], ", ")
print(io, "Grad($text)")
end
γ = 42.5774688MHz/T
# Defining gradient units
MagneticFieldGradient = Union{Quantity{T, 𝐈^-1*𝐌*𝐓^-2*𝐋^-1, U}, Level{L, S, Quantity{T, 𝐈^-1*𝐌*𝐓^-2*𝐋^-1, U}} where {L, S}} where {T, U}
# Promotion rules
Unitful.promote_to_derived()
Unitful.@derived_dimension BFieldGradient 𝐈^-1*𝐌*𝐓^-2*𝐋^-1 true
Unitful.@derived_dimension Frequency 𝐓^-1 true
Unitful.promote_unit(::S, ::T) where {S<:BFieldGradientFreeUnits, T<:BFieldGradientFreeUnits} = Unitful.T / Unitful.m
Unitful.promote_unit(::S, ::T) where {S<:FrequencyFreeUnits, T<:FrequencyFreeUnits} = Unitful.Hz
end
This would greatly improve the readability of the "Pulse Sequence" without needing to guess the units. Although we haven't tested it yet, it should not be difficult.