(WIP) Phantom time dependent properties

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

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
 

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

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)

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)

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

KomaMRIBase/src/motion/Interpolation.jl → 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