diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl
index 0f3624e11..e0808bfa1 100644
--- a/KomaMRIBase/src/KomaMRIBase.jl
+++ b/KomaMRIBase/src/KomaMRIBase.jl
@@ -70,6 +70,9 @@ include("datatypes/sequence/TimingChecks.jl")
include("datatypes/sequence/HardwareChecks.jl")
include("datatypes/sequence/AddBlockMacro.jl")
include("datatypes/sequence/DelayDuration.jl")
+# Timing
+include("timing/TimeCurve.jl")
+include("timing/Interpolation.jl")
# Motion
include("motion/MotionList.jl")
include("motion/NoMotion.jl")
@@ -102,8 +105,12 @@ export Translate, TranslateX, TranslateY, TranslateZ
export Rotate, RotateX, RotateY, RotateZ, CenterOfMass
export HeartBeat, Path, FlowPath
export TimeRange, Periodic, TimeCurve
-export SpinRange, AllSpins
+export SpinRange, AllSpins, spin_indicator
export get_spin_coords
+export TimeDependentProperty, get_phantom_property
+export get_spin_property, get_spin_property_at_end
+export get_spin_properties, get_spin_properties_block_end
+
# Secondary
export get_kspace, rotx, roty, rotz
# Additionals
diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl
index b09e97b7c..8236cfa5e 100644
--- a/KomaMRIBase/src/datatypes/Phantom.jl
+++ b/KomaMRIBase/src/datatypes/Phantom.jl
@@ -1,3 +1,5 @@
+include("phantom/TimeDependentProperty.jl")
+
"""
obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion)
@@ -9,14 +11,14 @@ a property value representing a spin. This struct serves as an input for the sim
- `x`: (`::AbstractVector{T<:Real}`, `[m]`) spin x-position vector
- `y`: (`::AbstractVector{T<:Real}`, `[m]`) spin y-position vector
- `z`: (`::AbstractVector{T<:Real}`, `[m]`) spin z-position vector
-- `ρ`: (`::AbstractVector{T<:Real}`) spin proton density vector
-- `T1`: (`::AbstractVector{T<:Real}`, `[s]`) spin T1 parameter vector
-- `T2`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2 parameter vector
-- `T2s`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2s parameter vector
-- `Δw`: (`::AbstractVector{T<:Real}`, `[rad/s]`) spin off-resonance parameter vector
-- `Dλ1`: (`::AbstractVector{T<:Real}`) spin Dλ1 (diffusion) parameter vector
-- `Dλ2`: (`::AbstractVector{T<:Real}`) spin Dλ2 (diffusion) parameter vector
-- `Dθ`: (`::AbstractVector{T<:Real}`) spin Dθ (diffusion) parameter vector
+- `ρ`: (`::PhantomProperty{T}`) spin proton density
+- `T1`: (`::PhantomProperty{T}`, `[s]`) spin T1 parameter vector
+- `T2`: (`::PhantomProperty{T}`, `[s]`) spin T2 parameter vector
+- `T2s`: (`::PhantomProperty{T}`, `[s]`) spin T2s parameter vector
+- `Δw`: (`::PhantomProperty{T}`, `[rad/s]`) spin off-resonance parameter vector
+- `Dλ1`: (`::PhantomProperty{T}`) spin Dλ1 (diffusion) parameter vector
+- `Dλ2`: (`::PhantomProperty{T}`) spin Dλ2 (diffusion) parameter vector
+- `Dθ`: (`::PhantomProperty{T}`) spin Dθ (diffusion) parameter vector
- `motion`: (`::Union{NoMotion, Motion{T<:Real} MotionList{T<:Real}}`) motion
# Returns
@@ -29,33 +31,40 @@ julia> obj = Phantom(x=[0.0])
julia> obj.ρ
```
"""
+const PhantomProperty{T<:Real} = Union{AbstractVector{T}, TimeDependentProperty{T}}
+
@with_kw mutable struct Phantom{T<:Real}
name::String = "spins"
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))
- 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
+ ρ::PhantomProperty{T} = ones(eltype(x), size(x))
+ T1::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000
+ T2::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000
+ T2s::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000
#Off-resonance related
- Δw::AbstractVector{T} = zeros(eltype(x), size(x))
+ Δw::PhantomProperty{T} = zeros(eltype(x), size(x))
#χ::Vector{SusceptibilityModel}
#Diffusion
- Dλ1::AbstractVector{T} = zeros(eltype(x), size(x))
- Dλ2::AbstractVector{T} = zeros(eltype(x), size(x))
- Dθ::AbstractVector{T} = zeros(eltype(x), size(x))
+ Dλ1::PhantomProperty{T} = zeros(eltype(x), size(x))
+ Dλ2::PhantomProperty{T} = zeros(eltype(x), size(x))
+ Dθ::PhantomProperty{T} = zeros(eltype(x), size(x))
#Diff::Vector{DiffusionModel} #Diffusion map
#Motion
motion::Union{NoMotion, Motion{T}, MotionList{T}} = NoMotion()
end
-const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom))
-const VECTOR_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) <: AbstractVector, fieldnames(Phantom))
+const POSITION_PHANTOM_FIELDS = (:x, :y, :z)
+const PROPERTY_PHANTOM_FIELDS = (:ρ, :T1, :T2, :T2s, :Δw, :Dλ1, :Dλ2, :Dθ)
+
+const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom))
+
+get_phantom_property(obj::Phantom, key::Symbol) = get_property_value(getfield(obj, key))
+get_property_value(p) = p
"""Size and length of a phantom"""
-size(x::Phantom) = size(x.ρ)
-Base.length(x::Phantom) = length(x.ρ)
+size(x::Phantom) = size(x.x)
+Base.length(x::Phantom) = length(x.x)
# To enable to iterate and broadcast over the Phantom
Base.iterate(x::Phantom) = (x[1], 2)
Base.iterate(x::Phantom, i::Integer) = (i <= length(x)) ? (x[i], i + 1) : nothing
@@ -95,20 +104,28 @@ end
"""Addition of phantoms"""
+(obj1::Phantom, obj2::Phantom) = begin
name = first(obj1.name * "+" * obj2.name, 50) # The name is limited to 50 characters
- fields = []
- for field in VECTOR_PHANTOM_FIELDS
- push!(fields, (field, [getfield(obj1, field); getfield(obj2, field)]))
+ fields = Pair{Symbol,Any}[]
+ for field in POSITION_PHANTOM_FIELDS
+ push!(fields, field => [getfield(obj1, field); getfield(obj2, field)])
+ end
+ for field in PROPERTY_PHANTOM_FIELDS
+ push!(fields, field => vcat_property(getfield(obj1, field), getfield(obj2, field)))
end
- return Phantom(;
- name = name,
- fields...,
- motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)))
+ return Phantom(;
+ name=name,
+ fields...,
+ motion=vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)),
+ )
end
"""Scalar multiplication of a phantom"""
*(α::Real, obj::Phantom) = begin
obj1 = deepcopy(obj)
- obj1.ρ .*= α
+ if obj1.ρ isa AbstractVector
+ obj1.ρ .*= α
+ else
+ obj1.ρ = α * obj1.ρ
+ end
return obj1
end
diff --git a/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl
new file mode 100644
index 000000000..c9dbb84a4
--- /dev/null
+++ b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl
@@ -0,0 +1,74 @@
+struct TimeDependentProperty{T<:Real}
+ value::AbstractMatrix{T}
+ time::TimeCurve{T}
+end
+
+get_property_value(p::TimeDependentProperty) = p.value[:, 1]
+
+"""Spin-resolved phantom properties at simulation times `t` (vector or matrix columns)."""
+function get_spin_properties(p, t)
+ return (
+ ρ=get_spin_property(p.ρ, t),
+ T1=get_spin_property(p.T1, t),
+ T2=get_spin_property(p.T2, t),
+ T2s=get_spin_property(p.T2s, t),
+ Δw=get_spin_property(p.Δw, t),
+ Dλ1=get_spin_property(p.Dλ1, t),
+ Dλ2=get_spin_property(p.Dλ2, t),
+ Dθ=get_spin_property(p.Dθ, t),
+ )
+end
+
+"""Block-end spin properties used for GPU kernels and block-final relaxation."""
+function get_spin_properties_block_end(p, t)
+ props = get_spin_properties(p, t)
+ return (
+ ρ=get_spin_property_at_end(props.ρ),
+ T1=get_spin_property_at_end(props.T1),
+ T2=get_spin_property_at_end(props.T2),
+ T2s=get_spin_property_at_end(props.T2s),
+ Δw=get_spin_property_at_end(props.Δw),
+ Dλ1=get_spin_property_at_end(props.Dλ1),
+ Dλ2=get_spin_property_at_end(props.Dλ2),
+ Dθ=get_spin_property_at_end(props.Dθ),
+ )
+end
+
+get_spin_property(p::AbstractVector, t) = p
+function get_spin_property(p::TimeDependentProperty, t)
+ Ns = size(p.value, 1)
+ t_unit = unit_time(t, p.time)
+ return resample(interpolate(p.value, Gridded(Linear()), Val(Ns), t_unit), t_unit)
+end
+
+"""
+ get_spin_property_at_end(p::AbstractMatrix)
+
+Returns spin-wise property values at the last queried time in `p`, with shape
+`(Nspins, Ntimes)`.
+"""
+get_spin_property_at_end(p::AbstractMatrix) = vec(selectdim(p, 2, size(p, 2)))
+get_spin_property_at_end(p::AbstractVector) = p
+
+Base.:(==)(p1::TimeDependentProperty, p2::TimeDependentProperty) = p1.value == p2.value && p1.time == p2.time
+Base.:(≈)(p1::TimeDependentProperty, p2::TimeDependentProperty) = p1.value ≈ p2.value && p1.time ≈ p2.time
+Base.:(==)(::TimeDependentProperty, ::AbstractVector) = false
+Base.:(==)(::AbstractVector, ::TimeDependentProperty) = false
+Base.:(≈)(::TimeDependentProperty, ::AbstractVector) = false
+Base.:(≈)(::AbstractVector, ::TimeDependentProperty) = false
+
+Base.getindex(p::TimeDependentProperty, i) = TimeDependentProperty(p.value[i, :], p.time)
+Base.view(p::TimeDependentProperty, i) = TimeDependentProperty(@view(p.value[i, :]), p.time)
+
+times(::AbstractVector) = Float64[]
+times(p::TimeDependentProperty) = times(p.time)
+
+*(α::Real, p::TimeDependentProperty) = TimeDependentProperty(α .* p.value, p.time)
+
+vcat_property(p1::AbstractVector, p2::AbstractVector) = [p1; p2]
+vcat_property(p1::TimeDependentProperty, p2::AbstractVector) = TimeDependentProperty([p1.value; repeat(p2, 1, size(p1.value, 2))], p1.time)
+vcat_property(p1::AbstractVector, p2::TimeDependentProperty) = TimeDependentProperty([repeat(p1, 1, size(p2.value, 2)); p2.value], p2.time)
+vcat_property(p1::TimeDependentProperty, p2::TimeDependentProperty) = begin
+ @assert p1.time == p2.time "Time ranges must be the same"
+ return TimeDependentProperty([p1.value; p2.value], p1.time)
+end
\ No newline at end of file
diff --git a/KomaMRIBase/src/motion/MotionList.jl b/KomaMRIBase/src/motion/MotionList.jl
index 99d6ecbad..941bdad5f 100644
--- a/KomaMRIBase/src/motion/MotionList.jl
+++ b/KomaMRIBase/src/motion/MotionList.jl
@@ -1,6 +1,4 @@
-include("Interpolation.jl")
include("SpinSpan.jl")
-include("TimeCurve.jl")
include("Action.jl")
include("Motion.jl")
@@ -54,7 +52,7 @@ end
# https://github.com/JuliaHealth/KomaMRI.jl/issues/480
""" Addition of MotionLists """
# MotionList + MotionList
-function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real}
+function vcat_motion(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real}
mv_aux = Motion{T}[]
for m in m1.motions
m_aux = deepcopy(m)
@@ -70,7 +68,7 @@ function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Rea
return MotionList(mv_aux...)
end
# Motion + Motion
-function Base.vcat(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real}
+function vcat_motion(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real}
mv_aux = Motion{T}[]
m_aux = deepcopy(m1)
m_aux.spins = expand(m_aux.spins, Ns1)
@@ -82,8 +80,8 @@ function Base.vcat(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real}
return MotionList(mv_aux...)
end
# Motion + MotionList
-Base.vcat(m1::MotionList{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} = vcat(m2, m1, Ns2, Ns1)
-function Base.vcat(m1::Motion{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real}
+vcat_motion(m1::MotionList{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} = vcat_motion(m2, m1, Ns2, Ns1)
+function vcat_motion(m1::Motion{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real}
mv_aux = Motion{T}[]
m_aux = deepcopy(m1)
m_aux.spins = expand(m_aux.spins, Ns1)
diff --git a/KomaMRIBase/src/motion/NoMotion.jl b/KomaMRIBase/src/motion/NoMotion.jl
index 3c5561c43..9a639bdbd 100644
--- a/KomaMRIBase/src/motion/NoMotion.jl
+++ b/KomaMRIBase/src/motion/NoMotion.jl
@@ -18,10 +18,10 @@ Base.view(mv::NoMotion, p) = mv
""" Addition of NoMotions """
# NoMotion + NoMotion
-Base.vcat(m1::NoMotion, m2::NoMotion, Ns1, Ns2) = m1
+vcat_motion(m1::NoMotion, m2::NoMotion, Ns1, Ns2) = m1
# NoMotion + MotionList
-Base.vcat(m1::MotionList, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1)
-function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T}
+vcat_motion(m1::MotionList, m2::NoMotion, Ns1, Ns2) = vcat_motion(m2, m1, 0, Ns1)
+function vcat_motion(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T}
mv_aux = Motion{T}[]
for m in m2.motions
m_aux = deepcopy(m)
@@ -32,8 +32,8 @@ function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T}
return MotionList(mv_aux)
end
# NoMotion + Motion
-Base.vcat(m1::Motion, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1)
-function Base.vcat(m1::NoMotion, m2::Motion{T}, Ns1, Ns2) where {T}
+vcat_motion(m1::Motion, m2::NoMotion, Ns1, Ns2) = vcat_motion(m2, m1, 0, Ns1)
+function vcat_motion(m1::NoMotion, m2::Motion{T}, Ns1, Ns2) where {T}
m_aux = deepcopy(m2)
m_aux.spins = expand(m_aux.spins, Ns2)
m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
diff --git a/KomaMRIBase/src/motion/SpinSpan.jl b/KomaMRIBase/src/motion/SpinSpan.jl
index fa6060aa3..cd0d3f003 100644
--- a/KomaMRIBase/src/motion/SpinSpan.jl
+++ b/KomaMRIBase/src/motion/SpinSpan.jl
@@ -72,3 +72,18 @@ expand(sr::SpinRange, Ns::Int) = sr
intersect_idx(a, b) = findall(x -> x in a, b)
intersect_idx(a, b::BitVector) = findall(x -> x in a, findall(x->x==true, b))
intersect_idx(a::BitVector, b) = findall(x -> x in findall(x->x==true, a), b)
+
+"""Per-spin weights for broadcast masking (GPU-friendly; no gather on `replace_by`)."""
+function spin_indicator(spins::AllSpins, w::AbstractVector{T}) where {T}
+ fill!(w, one(T))
+ return w
+end
+
+function spin_indicator(spins::SpinRange, w::AbstractVector{T}) where {T}
+ fill!(w, zero(T))
+ r = spins.range
+ @inbounds for i in r
+ i <= length(w) && (w[i] = one(T))
+ end
+ return w
+end
diff --git a/KomaMRIBase/src/motion/Interpolation.jl b/KomaMRIBase/src/timing/Interpolation.jl
similarity index 93%
rename from KomaMRIBase/src/motion/Interpolation.jl
rename to KomaMRIBase/src/timing/Interpolation.jl
index 09b699867..1227ca68d 100644
--- a/KomaMRIBase/src/motion/Interpolation.jl
+++ b/KomaMRIBase/src/timing/Interpolation.jl
@@ -54,6 +54,10 @@ function resample(itp::Interpolator2D, t)
Ns = size(itp.coefs, 1)
id = _similar(t, Ns)
copyto!(id, collect(range(oneunit(eltype(t)), eltype(t)(Ns), Ns)))
+ if t isa AbstractVector
+ # Produce an outer-product evaluation: Ns spins × length(t) times.
+ return itp.(reshape(id, :, 1), reshape(t, 1, :))
+ end
return itp.(id, t)
end
diff --git a/KomaMRIBase/src/motion/TimeCurve.jl b/KomaMRIBase/src/timing/TimeCurve.jl
similarity index 100%
rename from KomaMRIBase/src/motion/TimeCurve.jl
rename to KomaMRIBase/src/timing/TimeCurve.jl
diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl
index 17ac2c5fe..e3a7c9067 100644
--- a/KomaMRIBase/test/runtests.jl
+++ b/KomaMRIBase/test/runtests.jl
@@ -1527,7 +1527,7 @@ end
name, [x; x[rng]], [y; y[rng]], [z; z[rng]],
[ρ; ρ[rng]], [T1; T1[rng]], [T2; T2[rng]], [T2s; T2s[rng]],
[Δw; Δw[rng]], [Dλ1; Dλ1[rng]], [Dλ2; Dλ2[rng]], [Dθ; Dθ[rng]],
- vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
)
# NOTE: these vcat methods must be simplified once the Vector{<:Motion} approach is accomplished:
# https://github.com/JuliaHealth/KomaMRI.jl/issues/480
@@ -1535,42 +1535,61 @@ end
@test obj1 + obj2 == oba
# NoMotion + MotionList
obj2.motion = MotionList(tr, rt)
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# MotionList + NoMotion
obj1.motion = MotionList(tr, rt)
obj2.motion = NoMotion()
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# NoMotion + Motion
obj1.motion = NoMotion()
obj2.motion = tr
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# Motion + NoMotion
obj1.motion = tr
obj2.motion = NoMotion()
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# MotionList + MotionList
obj1.motion = MotionList(tr, rt)
obj2.motion = MotionList(tr, rt)
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# Motion + Motion
obj1.motion = tr
obj2.motion = rt
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# Motion + MotionList
obj1.motion = tr
obj2.motion = MotionList(tr, rt)
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
@test obj1 + obj2 == oba
# MotionList + Motion
obj1.motion = MotionList(tr, rt)
obj2.motion = tr
- oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ @test obj1 + obj2 == oba
+
+ tc = TimeRange(0.0, 1.0)
+ pd = [1.0 0.0; 0.5 1.0]
+ ρt = TimeDependentProperty(pd, tc)
+ # TimeDependentProperty + AbstractVector
+ obj1.ρ = ρt
+ obj2.ρ = ρ
+ oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ)
+ @test obj1 + obj2 == oba
+ # AbstractVector + TimeDependentProperty
+ obj1.ρ = ρ
+ obj2.ρ = ρt
+ oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ)
+ @test obj1 + obj2 == oba
+ # TimeDependentProperty + TimeDependentProperty
+ obj1.ρ = ρt
+ obj2.ρ = ρt
+ oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ)
@test obj1 + obj2 == oba
end
@testset "Scalar multiplication" begin
@@ -1579,6 +1598,24 @@ end
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θ, motion=MotionList(tr, rt))
@test c * obj1 == obc
end
+ @testset "Time-dependent ρ" begin
+ tc = TimeRange(0.0, 1.0)
+ pd = [1.0 0.0; 0.5 1.0]
+ ρt = TimeDependentProperty(pd, tc)
+ obj = Phantom(ρ=ρt, x=x[1:2], y=y[1:2], z=z[1:2])
+ @test KomaMRIBase.get_spin_property(obj.ρ, [0.5]) ≈ [0.5, 0.75]
+ ρt2 = KomaMRIBase.get_spin_property(obj.ρ, [0.0, 1.0])
+ @test ρt2 ≈ [1.0 0.0; 0.5 1.0]
+ @test KomaMRIBase.get_spin_property_at_end(ρt2) ≈ [0.0, 1.0]
+ end
+ @testset "Time-dependent T1" begin
+ tc = TimeRange(0.0, 1.0)
+ T1t = TimeDependentProperty([0.5 1.0; 1.5 2.0], tc)
+ obj = Phantom(T1=T1t, x=x[1:2], y=y[1:2], z=z[1:2])
+ props = KomaMRIBase.get_spin_properties(obj, [0.5])
+ @test props.T1 ≈ [0.75, 1.75]
+ @test KomaMRIBase.get_spin_properties_block_end(obj, [0.0, 1.0]).T1 ≈ [1.0, 2.0]
+ end
@testset "Brain Phantom 2D" begin
ph = brain_phantom2D()
@test ph.name == "brain2D_axial"
diff --git a/KomaMRICore/src/simulation/Flow.jl b/KomaMRICore/src/simulation/Flow.jl
index 99a782ba8..495cb8dd5 100644
--- a/KomaMRICore/src/simulation/Flow.jl
+++ b/KomaMRICore/src/simulation/Flow.jl
@@ -3,17 +3,29 @@ function outflow_spin_reset!(args...; kwargs...)
end
function outflow_spin_reset!(spin_state_matrix, t, ml::MotionList; replace_by=0, seq_t=0, add_t0=false)
- for m in ml.motions
+ for m in ml.motions
outflow_spin_reset!(spin_state_matrix, t, m; replace_by=replace_by, seq_t=seq_t, add_t0=add_t0)
end
return nothing
end
-function outflow_spin_reset!(spin_state_matrix, t, m::Motion; replace_by=0, seq_t=0, add_t0=false)
+function outflow_spin_reset!(spin_state_matrix, t, m::Motion; replace_by=0, seq_t=0, add_t0=false)
outflow_spin_reset!(spin_state_matrix, t, m.action, m.time, m.spins; replace_by=replace_by, seq_t=seq_t, add_t0=add_t0)
return nothing
end
+function _flow_spin_weights(spin_span, prototype::AbstractVector{T}) where {T}
+ w = similar(prototype)
+ return spin_indicator(spin_span, w)
+end
+
+_flow_reset_mask(mask::AbstractVector) = @view mask[lastindex(mask):lastindex(mask)]
+_flow_reset_mask(mask::AbstractMatrix) = selectdim(mask, 2, size(mask, 2))
+
+flow_replace_mag(replace_by::AbstractVector) = replace_by
+flow_replace_mag(replace_by::AbstractMatrix) = selectdim(replace_by, 2, size(replace_by, 2))
+flow_replace_mag(replace_by::Number) = replace_by
+
function outflow_spin_reset!(
spin_state_matrix::AbstractArray,
t,
@@ -24,17 +36,14 @@ function outflow_spin_reset!(
seq_t=0,
add_t0=false,
)
- # Initialize time: add t0 and normalize
ts = KomaMRIBase.unit_time(init_time(t, seq_t, add_t0), time_curve)
- # Get spin state range affected by the spin span
- idx = KomaMRIBase.get_indexing_range(spin_span)
- spin_state_matrix = @view(spin_state_matrix[idx, :])
- replace_by = replace_view(replace_by, idx)
- # Obtain mask
+ prototype = ndims(spin_state_matrix) == 1 ? spin_state_matrix : selectdim(spin_state_matrix, 2, 1)
+ w = _flow_spin_weights(spin_span, prototype)
mask = get_mask(action.spin_reset, ts)
- # Modify spin state: reset and replace by initial value
- spin_state_matrix .*= (1 .- mask)
- spin_state_matrix .+= replace_by .* mask
+ indicator = reshape(w, :, 1)
+ weighted_mask = mask .* indicator
+ spin_state_matrix .*= (1 .- weighted_mask)
+ spin_state_matrix .+= replace_by .* weighted_mask
return nothing
end
@@ -48,20 +57,14 @@ function outflow_spin_reset!(
seq_t=0,
add_t0=false,
)
- # Initialize time: add t0 and normalize
ts = KomaMRIBase.unit_time(init_time(t, seq_t, add_t0), time_curve)
- # Get spin state range affected by the spin span
- idx = KomaMRIBase.get_indexing_range(spin_span)
- M = @view(M[idx])
- replace_by = replace_view(replace_by, idx)
- # Obtain mask
- mask = get_mask(action.spin_reset, ts)
- mask = @view(mask[:, end])
- # Modify spin state: reset and replace by initial value
- M.xy .*= (1 .- mask)
- M.z .*= (1 .- mask)
- M.xy .+= 0 .* mask
- M.z .+= replace_by .* mask
+ w = _flow_spin_weights(spin_span, M.z)
+ mask = _flow_reset_mask(get_mask(action.spin_reset, ts))
+ ρ = flow_replace_mag(replace_by)
+ M.xy .*= (1 .- mask .* w)
+ M.z .*= (1 .- mask .* w)
+ M.xy .+= zero(eltype(M.z)) .* mask .* w
+ M.z .+= ρ .* mask .* w
return nothing
end
@@ -73,13 +76,6 @@ function init_time(t, seq_t, add_t0)
return t
end
-function replace_view(replace_by::AbstractArray, idx)
- return @view(replace_by[idx])
-end
-function replace_view(replace_by, idx)
- return replace_by
-end
-
function get_mask(spin_reset, t)
itp = KomaMRIBase.interpolate(spin_reset, KomaMRIBase.Gridded(KomaMRIBase.Constant{KomaMRIBase.Next}()), Val(size(spin_reset, 1)), t)
return KomaMRIBase.resample(itp, t)
diff --git a/KomaMRICore/src/simulation/Functors.jl b/KomaMRICore/src/simulation/Functors.jl
index 52dcc4fb2..4aac9de7d 100644
--- a/KomaMRICore/src/simulation/Functors.jl
+++ b/KomaMRICore/src/simulation/Functors.jl
@@ -112,6 +112,7 @@ adapt_storage(T::Type{<:Real}, xs::MotionList) = MotionList(paramtype.(T, xs.mot
@functor Path
@functor FlowPath
@functor TimeCurve
+@functor TimeDependentProperty
# Spinor
@functor Spinor
# DiscreteSequence
diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl
index bf3f6004b..f70009e30 100644
--- a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl
+++ b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl
@@ -33,7 +33,7 @@ function prealloc(sim_method::Bloch, backend::KA.CPU, obj::Phantom{T}, M::Mag{T}
similar(M.xy),
similar(M.xy)
),
- obj.Δw ./ T(2π .* γ)
+ zeros(T, size(obj.x)),
)
end
@@ -65,21 +65,28 @@ function run_spin_precession!(
fill!(ϕ, zero(T))
block_time = zero(T)
sample = 1
+ props_end = get_spin_properties_block_end(p, seq.t')
+ ρ_end = props_end.ρ
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[1])
+ Δw = get_spin_property(p.Δw, seq.t[1])
+ @. ΔBz = Δw / T(2π * γ)
@. Bz_old = x * seq.Gx[1] + y * seq.Gy[1] + z * seq.Gz[1] + ΔBz
#Simulation
for i in eachindex(seq.Δt)
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i + 1])
#Effective Field
+ Δw = get_spin_property(p.Δw, seq.t[i + 1])
+ @. ΔBz = Δw / T(2π * γ)
@. Bz_new = x * seq.Gx[i + 1] + y * seq.Gy[i + 1] + z * seq.Gz[i + 1] + ΔBz
#Rotation
@. ϕ += (Bz_old + Bz_new) * T(-π * γ) * seq.Δt[i]
block_time += seq.Δt[i]
#Acquired Signal
if seq.ADC[i + 1]
+ T2 = get_spin_property(p.T2, seq.t[i + 1])
#Update signal
- @. Mxy = exp(-block_time / p.T2) * M.xy * cis(ϕ)
+ @. Mxy = exp(-block_time / T2) * M.xy * cis(ϕ)
#Reset Spin-State (Magnetization). Only for FlowPath
outflow_spin_reset!(Mxy, seq.t[i + 1], p.motion)
#Acquired signal
@@ -90,10 +97,11 @@ function run_spin_precession!(
Bz_old .= Bz_new
end
#Final Spin-State
- @. M.xy = M.xy * exp(-block_time / p.T2) * cis(ϕ)
- @. M.z = M.z * exp(-block_time / p.T1) + p.ρ * (T(1) - exp(-block_time / p.T1))
+ T2_end = props_end.T2
+ @. M.xy = M.xy * exp(-block_time / T2_end) * cis(ϕ)
+ @. M.z = M.z * exp(-block_time / props_end.T1) + ρ_end * (T(1) - exp(-block_time / props_end.T1))
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end)
return nothing
end
@@ -133,8 +141,12 @@ function run_spin_excitation!(
for i in eachindex(seq.Δt)
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i])
+ props = get_spin_properties(p, seq.t[i, :])
+ ρ = props.ρ
+ ρ_end = get_spin_property_at_end(ρ)
#Effective field
- @. Bz = (seq.Gx[i] * x + seq.Gy[i] * y + seq.Gz[i] * z) + ΔBz - seq.Δf[i] / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z)
+ @. ΔBz = props.Δw / T(2π * γ) - seq.Δf[i] / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z)
+ @. Bz = seq.Gx[i] * x + seq.Gy[i] * y + seq.Gz[i] * z + ΔBz
@. B = sqrt(abs(seq.B1[i])^2 + abs(Bz)^2)
@. B[B == 0] = eps(T)
#Spinor Rotation
@@ -143,10 +155,10 @@ function run_spin_excitation!(
@. β = -Complex{T}(im) * (seq.B1[i] / B) * sin(φ_half)
mul!(Spinor(α, β), M, Maux_xy, Maux_z)
#Relaxation
- @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2)
- @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + p.ρ * (T(1) - exp(-seq.Δt[i] / p.T1))
+ @. M.xy = M.xy * exp(-seq.Δt[i] / props.T2)
+ @. M.z = M.z * exp(-seq.Δt[i] / props.T1) + ρ * (T(1) - exp(-seq.Δt[i] / props.T1))
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end)
#Acquire signal
if seq.ADC[i + 1] # ADC at the end of the time step
sig[sample] = sum(M.xy)
diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl
index dba40b09f..139f3443c 100644
--- a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl
+++ b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl
@@ -21,7 +21,7 @@ function prealloc(
return BlochGPUPrealloc(
KA.zeros(backend, Complex{T}, (cld(size(obj.x, 1), groupsize), max_block_length)),
KA.zeros(backend, Complex{T}, 1, max_block_length),
- obj.Δw ./ T(2π .* γ)
+ obj.Δw ./ T(2π .* γ),
)
end
@@ -37,13 +37,15 @@ function run_spin_precession!(
) where {T<:Real, SM<:BlochLikeSimMethods}
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
+ props_end = get_spin_properties_block_end(p, seq.t')
+ @. pre.ΔBz = props_end.Δw / T(2π * γ)
has_adc = length(sig) > 0
#Precession
precession_kernel!(backend, groupsize)(
pre.sig_output,
M.xy, M.z,
- x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)),
+ x, y, z, pre.ΔBz, props_end.T1, props_end.T2, props_end.ρ, UInt32(length(M.xy)),
seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.ADC, UInt32(length(seq.t)),
Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc),
sim_method,
@@ -57,7 +59,7 @@ function run_spin_precession!(
end
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t', p.motion; replace_by=props_end.ρ)
return nothing
end
@@ -74,13 +76,15 @@ function run_spin_excitation!(
) where {T<:Real, SM<:BlochLikeSimMethods}
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
+ props_end = get_spin_properties_block_end(p, seq.t')
+ @. pre.ΔBz = props_end.Δw / T(2π * γ)
has_adc = length(sig) > 0
#Excitation
excitation_kernel!(backend, groupsize)(
pre.sig_output,
M.xy, M.z,
- x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)),
+ x, y, z, pre.ΔBz, props_end.T1, props_end.T2, props_end.ρ, UInt32(length(M.xy)),
seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)),
Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc),
sim_method,
@@ -94,7 +98,7 @@ function run_spin_excitation!(
end
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ) # TODO: reset state inside kernel
+ outflow_spin_reset!(M, seq.t', p.motion; replace_by=props_end.ρ) # TODO: reset state inside kernel
return nothing
end
diff --git a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl
index 7db8d6f0f..52f7b9923 100644
--- a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl
+++ b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl
@@ -46,8 +46,12 @@ function run_spin_precession!(
) where {T<:Real}
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
+ t_acq = seq.t[2:end]'
+ props = get_spin_properties(p, t_acq)
+ ρ = props.ρ
+ ρ_end = get_spin_property_at_end(ρ)
#Effective field
- Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ)
+ Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ props.Δw ./ T(2π .* γ)
#Rotation
if is_ADC_on(seq)
ϕ = T(-2π .* γ) .* cumtrapz(seq.Δt', Bz)
@@ -57,7 +61,7 @@ function run_spin_precession!(
#Mxy precession and relaxation, and Mz relaxation
tp = cumsum(seq.Δt) # t' = t - t0
dur = sum(seq.Δt) # Total length, used for signal relaxation
- Mxy = M.xy .* exp.(-tp' ./ p.T2) .* cis.(ϕ) #This assumes Δw and T2 are constant in time
+ Mxy = M.xy .* exp.(-tp' ./ props.T2) .* cis.(ϕ)
M.xy .= Mxy[:, end]
#Reset Spin-State (Magnetization). Only for FlowPath
outflow_spin_reset!(Mxy, seq.t[2:end]', p.motion)
@@ -65,16 +69,17 @@ function run_spin_precession!(
sig[:, :, 1] .= @views transpose(Mxy[:, findall(seq.ADC[2:end])])
if sim_method.save_Mz
- Mz = M.z .* exp.(-tp' ./ p.T1) .+ p.ρ .* (1 .- exp.(-tp' ./ p.T1)) #Calculate intermediate points
+ Mz = M.z .* exp.(-tp' ./ props.T1) .+ ρ .* (1 .- exp.(-tp' ./ props.T1)) #Calculate intermediate points
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(Mz, seq.t[2:end]', p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(Mz, seq.t[2:end]', p.motion; replace_by=ρ_end)
sig[:, :, 2] .= @views transpose(Mz[:, findall(seq.ADC[2:end])]) #Save state to signal
M.z .= Mz[:, end]
else
- M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1)) #Jump to the last point
+ T1_end = get_spin_property_at_end(props.T1)
+ M.z .= M.z .* exp.(-dur ./ T1_end) .+ ρ_end .* (1 .- exp.(-dur ./ T1_end)) #Jump to the last point
end
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=ρ_end)
return nothing
end
diff --git a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl
index 4e65121ae..97acd8e4c 100644
--- a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl
+++ b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl
@@ -45,7 +45,7 @@ function prealloc(sim_method::BlochMagnus, backend::KA.CPU, obj::Phantom{T}, M::
similar(M.xy),
similar(M.xy)
),
- obj.Δw ./ T(2π .* γ)
+ zeros(T, size(obj.x)),
)
end
@@ -99,7 +99,11 @@ function run_spin_excitation!(
for i in eachindex(seq.Δt)
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i + 1])
+ props = get_spin_properties(p, seq.t[i, :])
+ ρ = props.ρ
+ ρ_end = get_spin_property_at_end(ρ)
#Effective field
+ @. ΔBz = props.Δw / T(2π * γ)
@. ωxy_new = seq.B1[i + 1] * B_to_ω
@. ωz_new = (seq.Gx[i + 1] * x + seq.Gy[i + 1] * y + seq.Gz[i + 1] * z + ΔBz) * B_to_ω + seq.Δf[i + 1] * T(2π)
effective_rotation_vector!(θxy, θz, ωxy_old, ωz_old, ωxy_new, ωz_new, seq.Δt[i], sim_method)
@@ -110,10 +114,10 @@ function run_spin_excitation!(
@. β = -Complex{T}(im) * (θxy / θ) * sin(θ / T(2))
mul!(Spinor(α, β), M, Maux_xy, Maux_z)
#Relaxation
- @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2)
- @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + p.ρ * (T(1) - exp(-seq.Δt[i] / p.T1))
+ @. M.xy = M.xy * exp(-seq.Δt[i] / props.T2)
+ @. M.z = M.z * exp(-seq.Δt[i] / props.T1) + ρ * (T(1) - exp(-seq.Δt[i] / props.T1))
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end)
#Acquire signal
if seq.ADC[i + 1] # ADC at the end of the time step
sig[sample] = sum(M.xy)
diff --git a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl
index 4c63d8b44..9fc2899ec 100644
--- a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl
+++ b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl
@@ -33,8 +33,10 @@ function run_spin_precession!(
) where {T<:Real}
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
+ props = get_spin_properties(p, seq.t')
+ ρ_end = get_spin_property_at_end(props.ρ)
#Effective field
- Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ)
+ Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ props.Δw ./ T(2π .* γ)
#Rotation
if is_ADC_on(seq)
ϕ = T(-2π .* γ) .* cumtrapz(seq.Δt', Bz)
@@ -44,12 +46,13 @@ function run_spin_precession!(
#Mxy precession and relaxation, and Mz relaxation
tp = cumsum(seq.Δt) # t' = t - t0
dur = sum(seq.Δt) # Total length, used for signal relaxation
- Mxy = M.xy .* exp.(-tp' ./ p.T2) .* cis.(ϕ) #This assumes Δw and T2 are constant in time
+ Mxy = M.xy .* exp.(-tp' ./ props.T2) .* cis.(ϕ)
+ Mz = M.z .* exp.(-dur ./ get_spin_property_at_end(props.T1)) .+ ρ_end .* (1 .- exp.(-dur ./ get_spin_property_at_end(props.T1)))
M.xy .= Mxy[:, end]
- M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1))
+ M.z .= Mz[:, end]
#Reset Spin-State (Magnetization). Only for FlowPath
outflow_spin_reset!(Mxy, seq.t[2:end]', p.motion)
- outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=ρ_end)
#Acquired signal
sig .= @views transpose(sum(Mxy[:, findall(seq.ADC[2:end])]; dims=1)) #<--- TODO: add coil sensitivities
return nothing
@@ -94,8 +97,11 @@ function run_spin_excitation!(
)
#Motion
x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, s.t)
+ props = get_spin_properties(p, s.t)
+ ρ = props.ρ
+ ρ_end = get_spin_property_at_end(ρ)
#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 = props.Δ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)
B .+= (B .== 0) .* eps(T)
@@ -103,10 +109,10 @@ function run_spin_excitation!(
φ = T(-2π .* γ) .* (B .* s.Δt)
mul!(Q(φ, s.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))
+ @. M.xy = M.xy * exp(-s.Δt / props.T2)
+ @. M.z = M.z * exp(-s.Δt / props.T1) + ρ * (1 - exp(-s.Δt / props.T1))
#Reset Spin-State (Magnetization). Only for FlowPath
- outflow_spin_reset!(M, s.tnew, p.motion; replace_by=p.ρ)
+ outflow_spin_reset!(M, s.tnew, p.motion; replace_by=ρ_end)
#Acquire signal
# TODO: use sim_method and sys to modify sig
if s.ADC # ADC at the end of the time step
diff --git a/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl b/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl
index 23b93c350..cdb4453cb 100644
--- a/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl
+++ b/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl
@@ -28,7 +28,7 @@ function initialize_spins_state(
) where {T<:Real}
Nspins = length(obj)
Mxy = zeros(T, Nspins)
- Mz = obj.ρ
+ Mz = get_phantom_property(obj, :ρ)
Xt = Mag{T}(Mxy, Mz)
return Xt, obj
end
diff --git a/KomaMRIFiles/src/Phantom/Phantom.jl b/KomaMRIFiles/src/Phantom/Phantom.jl
index 3589b7f8d..ba3c7eeee 100644
--- a/KomaMRIFiles/src/Phantom/Phantom.jl
+++ b/KomaMRIFiles/src/Phantom/Phantom.jl
@@ -19,12 +19,14 @@ function read_phantom(filename::String)
name = read_attribute(fid, "Name")
push!(phantom_fields, (:name, name))
# Position and contrast
- for key in ["position", "contrast"]
- group = fid[key]
- for label in HDF5.keys(group)
- values = read(group[label])
- push!(phantom_fields, (Symbol(label), values))
- end
+ pos = fid["position"]
+ for label in HDF5.keys(pos)
+ push!(phantom_fields, (Symbol(label), read(pos[label])))
+ end
+ contrast = fid["contrast"]
+ T = eltype(phantom_fields[end][2])
+ for label in HDF5.keys(contrast)
+ import_contrast_property!(phantom_fields, contrast, label, T)
end
# Motion
if "motion" in keys(fid)
@@ -88,6 +90,48 @@ function import_motion_subfield!(motion_subfields::Array, subfield_value::Union{
push!(motion_subfields, subfield_value)
return nothing
end
+function read_timecurve_group(group::HDF5.Group)
+ t = read(group, "t")
+ t_unit = read(group, "t_unit")
+ periodic = haskey(HDF5.attributes(group), "periodic") ? read_attribute(group, "periodic") : false
+ periods = haskey(HDF5.attributes(group), "periods") ? read_attribute(group, "periods") : 1.0
+ return TimeCurve(; t, t_unit, periodic, periods)
+end
+
+function write_timecurve_group!(group::HDF5.Group, tc::TimeCurve)
+ HDF5.attributes(group)["type"] = "TimeCurve"
+ group["t"] = tc.t
+ group["t_unit"] = tc.t_unit
+ HDF5.attributes(group)["periodic"] = tc.periodic
+ HDF5.attributes(group)["periods"] = tc.periods
+end
+
+function export_contrast_property!(contrast::HDF5.Group, label::Symbol, prop)
+ if prop isa TimeDependentProperty
+ grp = create_group(contrast, string(label))
+ HDF5.attributes(grp)["type"] = "TimeDependentProperty"
+ grp["value"] = prop.value
+ time_grp = create_group(grp, "time")
+ write_timecurve_group!(time_grp, prop.time)
+ else
+ contrast[string(label)] = prop
+ end
+ return nothing
+end
+
+function import_contrast_property!(phantom_fields::Array, contrast::HDF5.Group, label::String, ::Type{T}) where {T<:Real}
+ node = contrast[label]
+ if node isa HDF5.Dataset
+ push!(phantom_fields, (Symbol(label), read(node)))
+ else
+ @assert read_attribute(node, "type") == "TimeDependentProperty"
+ value = read(node["value"])
+ time = read_timecurve_group(node["time"])
+ push!(phantom_fields, (Symbol(label), TimeDependentProperty(value, time)))
+ end
+ return nothing
+end
+
function import_motion_subfield!(motion_subfields::Array, subfield_value::String, key::String, T::Type{<:Real})
if subfield_value in ["true", "false"]
return push!(motion_subfields, subfield_value == "true" ? true : false)
@@ -132,7 +176,7 @@ function write_phantom(
# Contrast (Rho, T1, T2, T2s Deltaw)
contrast = create_group(fid, "contrast")
for x in store_contrasts
- contrast[String(x)] = getfield(obj, x)
+ export_contrast_property!(contrast, x, getfield(obj, x))
end
# Motion
if !(obj.motion isa NoMotion) & store_motion
diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl
index 0a1a935c4..24d5fb5c9 100644
--- a/KomaMRIFiles/test/runtests.jl
+++ b/KomaMRIFiles/test/runtests.jl
@@ -41,6 +41,18 @@ end
obj2 = read_phantom(filename)
@test obj1 == obj2
end
+ @testset "TimeDependentProperty" begin
+ pth = @__DIR__
+ filename = pth * "/test_files/phantom/brain_timedep_w.phantom"
+ obj1 = brain_phantom2D()
+ tc = TimeRange(0.0, 1.0)
+ Ns = length(obj1)
+ obj1.ρ = TimeDependentProperty([obj1.ρ obj1.ρ], tc)
+ obj1.T1 = TimeDependentProperty([obj1.T1 2 .* obj1.T1], tc)
+ write_phantom(obj1, filename)
+ obj2 = read_phantom(filename)
+ @test obj1 == obj2
+ end
@testset "ArbitraryAction" begin
pth = @__DIR__
filename = pth * "/test_files/phantom/brain_arbitrarymotion_w.phantom"
diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl
index c696674cf..f6b22e04f 100644
--- a/KomaMRIPlots/src/ui/DisplayFunctions.jl
+++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl
@@ -1127,6 +1127,68 @@ function plot_kspace(seq::Sequence; width=nothing, height=nothing, darkmode=fals
return plot_koma(p, l; config)
end
+motion_time_nodes(::NoMotion) = Float64[]
+motion_time_nodes(motion) = times(motion)
+
+function append_time_dependent_property_nodes!(nodes, p::TimeDependentProperty)
+ append!(nodes, times(p))
+ return nothing
+end
+append_time_dependent_property_nodes!(nodes, ::Any) = nothing
+
+function property_time_nodes(obj::Phantom)
+ nodes = Float64[]
+ for field in fieldnames(Phantom)
+ append_time_dependent_property_nodes!(nodes, getfield(obj, field))
+ end
+ return nodes
+end
+
+function phantom_time_nodes(obj::Phantom)
+ return unique(sort([motion_time_nodes(obj.motion); property_time_nodes(obj)]))
+end
+
+phantom_time_node_count(obj::Phantom) = length(phantom_time_nodes(obj))
+
+property_extrema(p::AbstractVector) = extrema(p)
+property_extrema(p::TimeDependentProperty) = extrema(p.value)
+
+function spin_coords_for_plot(mv::NoMotion, x, y, z, t)
+ nt = length(t)
+ nt == 1 && return x, y, z
+ return repeat(x, 1, nt), repeat(y, 1, nt), repeat(z, 1, nt)
+end
+spin_coords_for_plot(motion, x, y, z, t) = get_spin_coords(motion, x, y, z, t')
+
+spin_property_frames(p::AbstractVector, t) =
+ repeat(reshape(KomaMRIBase.get_spin_property(p, t'), :, 1), 1, length(t))
+spin_property_frames(p::TimeDependentProperty, t) = KomaMRIBase.get_spin_property(p, t')
+
+function process_phantom_times(time_nodes, time_samples, ::Type{T}) where {T}
+ if isempty(time_nodes)
+ return [zero(T)]
+ elseif length(time_nodes) > 1
+ time_min, time_max = minimum(time_nodes), maximum(time_nodes)
+ t = collect(range(time_min, time_max, length=time_samples))
+ merged_times = sort(unique([t; time_nodes]))
+ if length(merged_times) > time_samples
+ to_keep = time_nodes
+ remaining = setdiff(merged_times, to_keep)
+ num_needed = time_samples - length(to_keep)
+ extra = remaining[round.(Int, range(1, length(remaining), length=num_needed))]
+ merged_times = sort([to_keep; extra])
+ elseif length(merged_times) < time_samples
+ extra_points = setdiff(t, merged_times)
+ extra_needed = time_samples - length(merged_times)
+ additional = sort(extra_points)[1:extra_needed]
+ merged_times = sort([merged_times; additional])
+ end
+ return merged_times
+ else
+ return time_nodes
+ end
+end
+
"""
p = plot_phantom_map(obj::Phantom, key::Symbol; kwargs...)
@@ -1171,37 +1233,9 @@ function plot_phantom_map(
view_2d=sum(KomaMRIBase.get_dims(obj)) < 3,
colorbar=true,
max_spins=20_000,
- time_samples= obj.motion isa NoMotion ? 0 : length(times(obj.motion)),
+ time_samples=phantom_time_node_count(obj),
kwargs...,
)
- function process_times(::NoMotion)
- return [zero(eltype(obj.x))]
- end
-
- function process_times(motion)
- time_nodes = times(motion)
- if length(time_nodes) > 1
- time_min, time_max = minimum(time_nodes), maximum(time_nodes)
- t = collect(range(time_min, time_max, length=time_samples))
- merged_times = sort(unique([t; time_nodes]))
- if length(merged_times) > time_samples
- to_keep = time_nodes
- remaining = setdiff(merged_times, to_keep)
- num_needed = time_samples - length(to_keep)
- extra = remaining[round.(Int, range(1, length(remaining), length=num_needed))]
- merged_times = sort([to_keep; extra])
- elseif length(merged_times) < time_samples
- extra_points = setdiff(t, merged_times)
- extra_needed = time_samples - length(merged_times)
- additional = sort(extra_points)[1:extra_needed]
- merged_times = sort([merged_times; additional])
- end
- return merged_times
- else
- return time_nodes
- end
- end
-
function decimate_uniform_phantom(obj, num_points::Int)
ss = Int(ceil(length(obj) / num_points))
return obj[1:ss:end]
@@ -1213,8 +1247,8 @@ function plot_phantom_map(
end
path = @__DIR__
- cmin_key = minimum(getproperty(obj, key))
- cmax_key = maximum(getproperty(obj, key))
+ prop = getproperty(obj, key)
+ cmin_key, cmax_key = property_extrema(prop)
if key == :T1 || key == :T2 || key == :T2s
cmin_key = 0
factor = 1e3
@@ -1253,8 +1287,9 @@ function plot_phantom_map(
cmin_key = get(kwargs, :zmin, factor * cmin_key)
cmax_key = get(kwargs, :zmax, factor * cmax_key)
- t = process_times(obj.motion)
- x, y, z = get_spin_coords(obj.motion, obj.x, obj.y, obj.z, t')
+ t = process_phantom_times(phantom_time_nodes(obj), time_samples, eltype(obj.x))
+ x, y, z = spin_coords_for_plot(obj.motion, obj.x, obj.y, obj.z, t)
+ values = spin_property_frames(prop, t)
x0 = -maximum(abs.([x y z])) * 1e2
xf = maximum(abs.([x y z])) * 1e2
@@ -1304,7 +1339,7 @@ function plot_phantom_map(
x=dims[1] ? (x[:,i])*1e2 : (y[:,i])*1e2,
y=dims[1] & dims[2] ? (y[:,i])*1e2 : (z[:,i])*1e2,
mode="markers",
- marker=attr(color=getproperty(obj,key)*factor,
+ marker=attr(color=values[:, i]*factor,
showscale=colorbar,
colorscale=colormap,
colorbar=attr(ticksuffix=unit, title=string(key)),
@@ -1314,7 +1349,7 @@ function plot_phantom_map(
),
visible=i==1,
showlegend=false,
- text=round.(getproperty(obj,key)*factor,digits=4),
+ text=round.(values[:, i]*factor,digits=4),
hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
$(string(key)): %{text}$unit"))
end
else # 3D
@@ -1357,7 +1392,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=values[:, i]*factor,
showscale=colorbar,
colorscale=colormap,
colorbar=attr(ticksuffix=unit, title=string(key)),
@@ -1367,7 +1402,7 @@ function plot_phantom_map(
),
visible=i==1,
showlegend=false,
- text=round.(getproperty(obj,key)*factor,digits=4),
+ text=round.(values[:, i]*factor,digits=4),
hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
z: %{z:.1f} cm
$(string(key)): %{text}$unit"))
end
end