Add key time points in DiscreteSequence when evaluating periodic and flow-related motions

KomaMRIBase/Project.toml

@@ -1,6 +1,7 @@
 name = "KomaMRIBase"
 uuid = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
-version = "0.9.6"
+
+version = "0.9.7"
 authors = ["Carlos Castillo Passi <[email protected]>"]
 
 [workspace]

KomaMRIBase/src/KomaMRIBase.jl

@@ -53,7 +53,7 @@ export brain_phantom2D, brain_phantom3D, pelvis_phantom2D, heart_phantom
 export MotionList, NoMotion, Motion
 export translate, rotate, heartbeat, path, flowpath
 export Translate, TranslateX, TranslateY, TranslateZ
-export Rotate, RotateX, RotateY, RotateZ 
+export Rotate, RotateX, RotateY, RotateZ, CenterOfMass
 export HeartBeat, Path, FlowPath
 export TimeRange, Periodic, TimeCurve
 export SpinRange, AllSpins

KomaMRIBase/src/datatypes/Phantom.jl

@@ -72,6 +72,9 @@ function Base.:(≈)(obj1::Phantom, obj2::Phantom)
     if length(obj1) != length(obj2) return false end
     return reduce(&, [getfield(obj1, field)  ≈ getfield(obj2, field) for field in NON_STRING_PHANTOM_FIELDS])
 end
+# Comparison between two different motion types is always false:
+Base.:(==)(::Union{NoMotion, Motion, MotionList}, ::Union{NoMotion, Motion, MotionList}) = false
+Base.:(≈)(::Union{NoMotion, Motion, MotionList},  ::Union{NoMotion, Motion, MotionList}) = false
 
 """Separate object spins in a sub-group"""
 function Base.getindex(obj::Phantom, p)

KomaMRIBase/src/motion/Action.jl

@@ -1,8 +1,7 @@
 abstract type AbstractAction{T<:Real} end
 
-Base.:(≈)(::Nothing, ::Nothing) = true
-Base.:(≈)(::Nothing, _) = false; Base.:(≈)(_, ::Nothing) = false
-Base.:(≈)(a1::AbstractAction, a2::AbstractAction) = (typeof(a1) == typeof(a2)) & reduce(&, [getfield(a1, field)  ≈ getfield(a2, field) for field in fieldnames(typeof(a1))])
+Base.:(==)(a1::AbstractAction, a2::AbstractAction) = (typeof(a1) == typeof(a2)) & reduce(&, [getfield(a1, field) == getfield(a2, field) for field in fieldnames(typeof(a1))])
+Base.:(≈)(a1::AbstractAction,  a2::AbstractAction) = (typeof(a1) == typeof(a2)) & reduce(&, [getfield(a1, field)  ≈ getfield(a2, field) for field in fieldnames(typeof(a1))])
 is_composable(m::AbstractAction) = true
 
 # Simple actions

KomaMRIBase/src/motion/Motion.jl

@@ -80,7 +80,7 @@ end
 - `spins`: (`::AbstractSpinSpan`) spin indexes affected by the motion
 
 # Keywords
-- `center`: (`::NTuple{3,Real}` or `nothing`) optional center of rotation, given in global coordinates. If `nothing` (default), the rotation is performed around the phantom’s center of mass.
+- `center`: (`::NTuple{3,Real}` or `::CenterOfMass`) center of rotation, given in global coordinates. Default is center of mass.
 
 # Returns
 - `rt`: (`::Motion`) Motion struct with [`Rotate`](@ref) action
@@ -90,7 +90,7 @@ end
 julia> rt = rotate(15.0, 0.0, 20.0, TimeRange(0.0, 1.0), SpinRange(1:10))
 ```
 """
-function rotate(pitch, roll, yaw, time=TimeRange(t_start=zero(eltype(pitch)), t_end=eps(eltype(pitch))), spins=AllSpins(); center=nothing)
+function rotate(pitch, roll, yaw, time=TimeRange(t_start=zero(eltype(pitch)), t_end=eps(eltype(pitch))), spins=AllSpins(); center=CenterOfMass())
     return Motion(Rotate(pitch, roll, yaw, center), time, spins)
 end
 
@@ -190,12 +190,25 @@ end
 
 """
     x, y, z = get_spin_coords(motion, x, y, z, t)
+
+Calculates the position of each spin at a set of arbitrary time instants, i.e. the time steps of the simulation. 
+For each dimension (x, y, z), the output matrix has ``N_{\t{spins}}`` rows and `length(t)` columns.
+
+# Arguments
+- `motion`: (`::Union{NoMotion, Motion{T<:Real} MotionList{T<:Real}}`) phantom motion
+- `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
+- `t`: horizontal array of time instants
+
+# Returns
+- `x, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time
 """
 function get_spin_coords(
     m::Motion{T}, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, t
 ) where {T<:Real}
     ux, uy, uz = x .* (0*t), y .* (0*t), z .* (0*t) # Buffers for displacements
-    t_unit = unit_time(t, m.time.t, m.time.t_unit, m.time.periodic, m.time.periods)
+    t_unit = unit_time(t, m.time)
     idx = get_indexing_range(m.spins)
     displacement_x!(@view(ux[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)
     displacement_y!(@view(uy[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)
@@ -204,7 +217,58 @@ function get_spin_coords(
 end
 
 # Auxiliary functions
-times(m::Motion) = times(m.time.t, m.time.periods)
+times(m::Motion) = times(m.time)
 is_composable(m::Motion) = is_composable(m.action)
-add_jump_times!(t, m::Motion) = add_jump_times!(t, m.action, m.time)
-add_jump_times!(t, ::AbstractAction, ::TimeCurve) = nothing
+
+"""
+    add_key_time_points!(t, motion)
+"""
+function add_key_time_points!(t, m::Motion)
+    add_key_time_points!(t, m.action, m.time.t_start, m.time.t_end, m.time.periods, m.time.periodic)
+    return nothing
+end
+function add_key_time_points!(t, a, t_start::T, t_end::T, periods, periodic) where T
+    isempty(t) && return
+    aux = T[] 
+    period = sum((t_end - t_start) .* periods)
+    t_max = maximum(t)
+    add_period_times!(aux, t_start, t_end, periods)
+    add_reset_times!(aux, a, t_start, t_end, periods)
+    extend_periodic!(aux, t_max, period, Val(periodic))
+    append!(t, aux[aux .<= t_max])
+    return nothing
+end
+
+"""
+    extend_periodic!(aux, t_max, period, periodic)
+"""
+function extend_periodic!(aux, t_max, period, periodic::Val{false})
+    return nothing
+end
+function extend_periodic!(aux, t_max, period, periodic::Val{true})
+    n_periods = floor(Int, t_max / period)
+    if n_periods > 0
+        initial_size = length(aux)
+        sizehint!(aux, initial_size * (n_periods + 1))
+        for n in 1:n_periods
+            append!(aux, aux[1:initial_size] .+ n*period)
+        end
+    end
+    return nothing
+end
+
+"""
+    add_period_times!(t, t_start, t_end, periods)
+"""
+function add_period_times!(t, t_start, t_end, periods)
+    period_times = times([t_start, t_end], t_start, t_end, periods)
+    append!(t, period_times .+ MIN_RISE_TIME .* ((-1) .^ ((1:length(period_times)) .+ 1)))
+    return nothing
+end
+
+"""
+    add_reset_times!(t, action, t_start, t_end, periods)
+""" 
+function add_reset_times!(t, ::AbstractAction, t_start, t_end, periods)
+    return nothing 
+end

KomaMRIBase/src/motion/MotionList.jl

@@ -130,22 +130,6 @@ end
 """ MotionList length """
 Base.length(m::MotionList) = length(m.motions)
 
-"""
-    x, y, z = get_spin_coords(motionset, x, y, z, t)
-
-Calculates the position of each spin at a set of arbitrary time instants, i.e. the time steps of the simulation. 
-For each dimension (x, y, z), the output matrix has ``N_{\t{spins}}`` rows and `length(t)` columns.
-
-# Arguments
-- `motion`: (`::Union{NoMotion, MotionList{T<:Real}}`) phantom motion
-- `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
-- `t`: horizontal array of time instants
-
-# Returns
-- `x, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time
-"""
 function get_spin_coords(
     ml::MotionList{T}, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, t
 ) where {T<:Real}
@@ -157,7 +141,7 @@ function get_spin_coords(
     ux, uy, uz = xt .* zero(T), yt .* zero(T), zt .* zero(T)
     # Composable motions: they need to be run sequentially. Note that they depend on xt, yt, and zt
     for m in Iterators.filter(is_composable, ml.motions)
-        t_unit = unit_time(t, m.time.t, m.time.t_unit, m.time.periodic, m.time.periods)
+        t_unit = unit_time(t, m.time)
         idx = get_indexing_range(m.spins)
         displacement_x!(@view(ux[idx, :]), m.action, @view(xt[idx, :]), @view(yt[idx, :]), @view(zt[idx, :]), t_unit)
         displacement_y!(@view(uy[idx, :]), m.action, @view(xt[idx, :]), @view(yt[idx, :]), @view(zt[idx, :]), t_unit)
@@ -167,7 +151,7 @@ function get_spin_coords(
     end
     # Additive motions: these motions can be run in parallel
     for m in Iterators.filter(!is_composable, ml.motions)
-        t_unit = unit_time(t, m.time.t, m.time.t_unit, m.time.periodic, m.time.periods)
+        t_unit = unit_time(t, m.time)
         idx = get_indexing_range(m.spins)
         displacement_x!(@view(ux[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)
         displacement_y!(@view(uy[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)
@@ -182,19 +166,19 @@ end
     times = times(motion)
 """
 function times(ml::MotionList)
-    nodes = reduce(vcat, [times(m) for m in ml.motions])
+    nodes = reduce(vcat, [times(m.time) for m in ml.motions])
     return unique(sort(nodes))
 end
 
 """
     sort_motions!(motion)
 
 Sorts motions in a list according to their starting time. It modifies the original list.
-If `motionset::NoMotion`, this function does nothing.
-If `motionset::MotionList`, this function sorts its motions.
+If `motion::Union{NoMotion, Motion}`, this function does nothing.
+If `motion::MotionList`, this function sorts its motions.
 
 # Arguments
-- `motion`: (`::Union{NoMotion, MotionList{T<:Real}}`) phantom motion
+- `motion`: (`::Union{NoMotion, Motion{T<:Real} MotionList{T<:Real}}`) phantom motion
 
 # Returns
 - `nothing`
@@ -204,8 +188,8 @@ function sort_motions!(m::MotionList)
     return nothing
 end
 
-function add_jump_times!(t, ml::MotionList)
+function add_key_time_points!(t, ml::MotionList)
     for m in ml.motions
-        add_jump_times!(t, m)
+        add_key_time_points!(t, m)
     end
 end

KomaMRIBase/src/motion/NoMotion.jl

@@ -49,4 +49,4 @@ function get_spin_coords(
 ) where {T<:Real}
     return x, y, z
 end
-add_jump_times!(t, ::NoMotion) = nothing
+add_key_time_points!(t, ::NoMotion) = nothing

KomaMRIBase/src/motion/TimeCurve.jl

@@ -129,20 +129,19 @@ Base.:(==)(t1::TimeCurve, t2::TimeCurve) = reduce(&, [getfield(t1, field) == get
 Base.:(≈)(t1::TimeCurve, t2::TimeCurve)  = reduce(&, [getfield(t1, field)  ≈ getfield(t2, field) for field in fieldnames(typeof(t1))])
 
 """ times & unit_time """
-# Although the implementation of these two functions when `per` is a vector is valid 
-# for all cases, it performs unnecessary and costly operations when `per` is a scalar. 
+times(tc::TimeCurve) = times(tc.t, tc.t_start, tc.t_end, tc.periods)
+# Although the implementation of these two functions when `periods` is a vector is valid 
+# for all cases, it performs unnecessary and costly operations when `periods` is a scalar.
 # Therefore, it has been decided to use method dispatch between these two cases.
-function times(t, per::Real)
-    return per .* t
+function times(t, t_start, t_end, periods::Real)
+    return t_start .+ periods .* (t .- t_start)
 end
-function times(t, per::AbstractVector)
-    scale   = repeat(per, inner=[length(t)])
-    offsets = repeat(cumsum(vcat(0, per[1:end-1]*t[end])), inner=[length(t)])
-    return (repeat(t, length(per)) .* scale) .+ offsets
-end
-function unit_time(tq, t, t_unit, periodic, per::Real)
-    return interpolate_times(t .* per, t_unit, periodic, tq)
-end
-function unit_time(tq, t, t_unit, periodic, per::AbstractVector)
-    return interpolate_times(times(t, per), repeat(t_unit, length(per)), periodic, tq)
+function times(t, t_start, t_end, periods::AbstractVector)
+    scale   = repeat(periods, inner=[length(t)])
+    offsets = repeat(cumsum(vcat(0, periods[1:end-1]*(t_end - t_start))), inner=[length(t)])
+    return t_start .+ ((repeat(t, length(periods)) .- t_start).* scale) .+ offsets
 end
+
+function unit_time(tq, tc::TimeCurve)
+    return interpolate_times(times(tc), repeat(tc.t_unit, length(tc.periods)), tc.periodic, tq)
+end

KomaMRIBase/src/motion/actions/ArbitraryAction.jl

@@ -7,9 +7,6 @@ function Base.view(action::ArbitraryAction, p)
     return typeof(action)([@view(getfield(action, d)[p,:]) for d in fieldnames(typeof(action))]...)
 end
 
-Base.:(==)(m1::ArbitraryAction, m2::ArbitraryAction) = (typeof(m1) == typeof(m2)) & reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(typeof(m1))])
-Base.:(≈)(m1::ArbitraryAction, m2::ArbitraryAction)  = (typeof(m1) == typeof(m2)) & reduce(&, [getfield(m1, field)  ≈ getfield(m2, field) for field in fieldnames(typeof(m1))])
-
 function displacement_x!(ux, action::ArbitraryAction, x, y, z, t)
     itp = interpolate(action.dx, Gridded(Linear()), Val(size(action.dx,1)), t)
     ux .= resample(itp, t)

KomaMRIBase/src/motion/actions/arbitraryactions/FlowPath.jl

@@ -38,7 +38,7 @@ end
 
 FlowPath(dx::AbstractArray{T}, dy::AbstractArray{T}, dz::AbstractArray{T}, spin_reset::BitMatrix) where T<:Real = FlowPath(dx, dy, dz, collect(spin_reset))
 
-function add_jump_times!(t, a::FlowPath, tc::TimeCurve) 
-    jump_times = (tc.t_end - tc.t_start)/(size(a.spin_reset)[2]-1) * (getindex.(findall(a.spin_reset .== 1), 2) .- 1) .- 1e-6
-    append!(t, jump_times)
+function add_reset_times!(t, a::FlowPath, t_start, t_end, periods)
+    aux = t_start .+ (t_end - t_start)/(size(a.spin_reset)[2]-1) * (getindex.(findall(a.spin_reset .== 1), 2) .- 1)
+    append!(t, times(aux, t_start, t_end, periods) .- MIN_RISE_TIME)
 end

KomaMRIBase/src/motion/actions/simpleactions/Rotate.jl

@@ -1,5 +1,10 @@
+struct CenterOfMass end
+Base.:(≈)(::CenterOfMass, ::CenterOfMass) = true
+Base.:(≈)(a::CenterOfMass, b) = false
+Base.:(≈)(a, b::CenterOfMass) = false
+
 @doc raw"""
-    r = Rotate(pitch, roll, yaw, center=nothing)
+    r = Rotate(pitch, roll, yaw, center=CenterOfMass())
 
 Rotate struct. It produces a rotation in the three axes: 
 x (pitch), y (roll), and z (yaw).
@@ -41,11 +46,11 @@ R &= R_z(\alpha) R_y(\beta) R_x(\gamma) \\
 - `pitch`: (`::Real`, `[º]`) rotation in x
 - `roll`: (`::Real`, `[º]`) rotation in y 
 - `yaw`: (`::Real`, `[º]`) rotation in z
-- `center`: (`::NTuple{3,Real}` or `nothing`) optional center of rotation, given in global coordinates. If `nothing` (default), the rotation is performed around the phantom’s center of mass.
+- `center`: (`::NTuple{3,Real}` or `::CenterOfMass`) optional center of rotation, given in global coordinates. Default is center of mass.
 
 # Notes
 - Rotations are applied around the point specified in `center`. If omitted, the rotation is centered at the phantom’s center of mass.
-- If `center` is not null, the rotation center is interpreted as a fixed point in space (absolute/global coordinates).
+- If `center` is not `::CenterOfMass`, the rotation center is interpreted as a fixed point in space (absolute/global coordinates).
 - This design ensures that consecutive or inverse rotations behave consistently and predictably, since the rotation center does not change with object transformations.
 
 # Returns
@@ -63,25 +68,24 @@ julia> r = Rotate(pitch=0.0, roll=45.0, yaw=0.0, center=(5e-3,0.0,0.0))
     pitch      :: T
     roll       :: T
     yaw        :: T
-    center     :: Union{Nothing,NTuple{3,T}} = nothing
+    center     :: Union{CenterOfMass,NTuple{3,T}} = CenterOfMass()
 end
 
-RotateX(pitch::T) where {T<:Real} = Rotate(pitch, zero(T), zero(T))
-RotateY(roll::T) where {T<:Real}  = Rotate(zero(T), roll, zero(T))
-RotateZ(yaw::T) where {T<:Real}   = Rotate(zero(T), zero(T), yaw)
+RotateX(pitch::T) where {T<:Real} = Rotate(pitch=pitch,   roll=zero(T), yaw=zero(T))
+RotateY(roll::T)  where {T<:Real} = Rotate(pitch=zero(T), roll=roll,    yaw=zero(T))
+RotateZ(yaw::T)   where {T<:Real} = Rotate(pitch=zero(T), roll=zero(T), yaw=yaw)
+
+get_center(center::CenterOfMass, x, y, z) = (sum(x) / length(x), sum(y) / length(y), sum(z) / length(z))
+get_center(center::NTuple, x, y, z)       = center
 
-function displacement_x!(ux, action::Rotate, x, y, z, t)   
+function displacement_x!(ux, action::Rotate, x, y, z, t)
     # Not using sind and cosd functions until bug with oneAPI is solved: 
     # https://github.com/JuliaGPU/oneAPI.jl/issues/65
     α = t .* (action.yaw*π/180)
     β = t .* (action.roll*π/180)
     γ = t .* (action.pitch*π/180)
-    cx = isnothing(action.center) ? sum(x) / length(x) : action.center[1]
-    cy = isnothing(action.center) ? sum(y) / length(y) : action.center[2]
-    cz = isnothing(action.center) ? sum(z) / length(z) : action.center[3]
-    x0 = x .- cx
-    y0 = y .- cy
-    z0 = z .- cz
+    cx, cy, cz = get_center(action.center, x, y, z)
+    x0, y0, z0 = x .- cx, y .- cy, z .- cz
     ux .= cos.(α) .* cos.(β) .* x0 +
          (cos.(α) .* sin.(β) .* sin.(γ) .- sin.(α) .* cos.(γ)) .* y0 +
          (cos.(α) .* sin.(β) .* cos.(γ) .+ sin.(α) .* sin.(γ)) .* z0 .+ cx .- x
@@ -92,12 +96,8 @@ function displacement_y!(uy, action::Rotate, x, y, z, t)
     α = t .* (action.yaw*π/180)
     β = t .* (action.roll*π/180)
     γ = t .* (action.pitch*π/180)
-    cx = isnothing(action.center) ? sum(x) / length(x) : action.center[1]
-    cy = isnothing(action.center) ? sum(y) / length(y) : action.center[2]
-    cz = isnothing(action.center) ? sum(z) / length(z) : action.center[3]
-    x0 = x .- cx
-    y0 = y .- cy
-    z0 = z .- cz
+    cx, cy, cz = get_center(action.center, x, y, z)
+    x0, y0, z0 = x .- cx, y .- cy, z .- cz
     uy .= sin.(α) .* cos.(β) .* x0 +
          (sin.(α) .* sin.(β) .* sin.(γ) .+ cos.(α) .* cos.(γ)) .* y0 +
          (sin.(α) .* sin.(β) .* cos.(γ) .- cos.(α) .* sin.(γ)) .* z0 .+ cy .- y
@@ -108,12 +108,8 @@ function displacement_z!(uz, action::Rotate, x, y, z, t)
     α = t .* (action.yaw*π/180)
     β = t .* (action.roll*π/180)
     γ = t .* (action.pitch*π/180)
-    cx = isnothing(action.center) ? sum(x) / length(x) : action.center[1]
-    cy = isnothing(action.center) ? sum(y) / length(y) : action.center[2]
-    cz = isnothing(action.center) ? sum(z) / length(z) : action.center[3]
-    x0 = x .- cx
-    y0 = y .- cy
-    z0 = z .- cz
+    cx, cy, cz = get_center(action.center, x, y, z)
+    x0, y0, z0 = x .- cx, y .- cy, z .- cz
     uz .=  -sin.(β) .* x0 + 
             cos.(β) .* sin.(γ) .* y0 +
             cos.(β) .* cos.(γ) .* z0 .+ cz .- z

KomaMRIBase/src/timing/TimeStepCalculation.jl

@@ -122,7 +122,7 @@ function get_variable_times(seq; Δt=1e-3, Δt_rf=1e-5, motion=NoMotion())
         end
         append!(t, t_block)
 	end
-	add_jump_times!(t, motion)
+	add_key_time_points!(t, motion)
 	# Removing repeated points
 	sort!(unique!(t))
 	# Fixes a problem with ADC at the start and end of the seq