Merge branch 'master' into docs-develop

KomaMRIBase/src/KomaMRIBase.jl

@@ -26,8 +26,8 @@ include("timing/KeyValuesCalculation.jl")
 include("datatypes/Sequence.jl")
 include("datatypes/sequence/Delay.jl")
 # Motion
-include("motion/motionlist/MotionList.jl")
-include("motion/nomotion/NoMotion.jl")
+include("motion/MotionList.jl")
+include("motion/NoMotion.jl")
 # Phantom
 include("datatypes/Phantom.jl")
 # Simulator

KomaMRIBase/src/datatypes/Phantom.jl

@@ -17,7 +17,7 @@ a property value representing a spin. This struct serves as an input for the sim
 - `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
-- `motion`: (`::Union{NoMotion, MotionList{T<:Real}}`) motion
+- `motion`: (`::Union{NoMotion, Motion{T<:Real} MotionList{T<:Real}}`) motion
 
 # Returns
 - `obj`: (`::Phantom`) Phantom struct
@@ -47,7 +47,7 @@ julia> obj.ρ
     Dθ::AbstractVector{T}  = zeros(eltype(x), size(x))
     #Diff::Vector{DiffusionModel}  #Diffusion map
     #Motion
-    motion::Union{NoMotion, MotionList{T}} = NoMotion()
+    motion::Union{NoMotion, Motion{T}, MotionList{T}} = NoMotion()
 end
 
 const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String,         fieldnames(Phantom))

KomaMRIBase/src/datatypes/simulation/DiscreteSequence.jl

@@ -85,8 +85,8 @@ based on simulation parameters.
 # Returns
 - `seqd`: (`::DiscreteSequence`) DiscreteSequence struct
 """
-function discretize(seq::Sequence; sampling_params=default_sampling_params())
-    t, Δt      = get_variable_times(seq; Δt=sampling_params["Δt"], Δt_rf=sampling_params["Δt_rf"])
+function discretize(seq::Sequence; sampling_params=default_sampling_params(), motion=NoMotion())
+    t, Δt      = get_variable_times(seq; Δt=sampling_params["Δt"], Δt_rf=sampling_params["Δt_rf"], motion=motion)
     B1, Δf     = get_rfs(seq, t)
     Gx, Gy, Gz = get_grads(seq, t)
     tadc       = get_adc_sampling_times(seq)

KomaMRIBase/src/motion/motionlist/Motion.jl → KomaMRIBase/src/motion/Motion.jl

@@ -178,13 +178,30 @@ Base.:(≈)(m1::Motion, m2::Motion)  = (typeof(m1) == typeof(m2)) & reduce(&, [g
 """ Motion sub-group """
 function Base.getindex(m::Motion, p)
     idx, spin_range = m.spins[p]
-    return idx !== nothing ? Motion(m.action[idx], m.time, spin_range) : nothing
+    return idx !== nothing ? Motion(m.action[idx], m.time, spin_range) : NoMotion()
 end
 function Base.view(m::Motion, p)
     idx, spin_range = @view(m.spins[p])
-    return idx !== nothing ? Motion(@view(m.action[idx]), m.time, spin_range) : nothing
+    return idx !== nothing ? Motion(@view(m.action[idx]), m.time, spin_range) : NoMotion()
+end
+
+"""
+    x, y, z = get_spin_coords(motion, x, y, z, t)
+"""
+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)
+    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)
+    displacement_z!(@view(uz[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)
+    return x .+ ux, y .+ uy, z .+ uz
 end
 
 # Auxiliary functions
 times(m::Motion) = times(m.time)
-is_composable(m::Motion) = is_composable(m.action)
+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, ::AbstractTimeSpan) = nothing

KomaMRIBase/src/motion/motionlist/MotionList.jl → KomaMRIBase/src/motion/MotionList.jl

@@ -1,6 +1,6 @@
-include("Action.jl")
 include("SpinSpan.jl")
 include("TimeSpan.jl")
+include("Action.jl")
 include("Motion.jl")
 
 """
@@ -36,27 +36,24 @@ struct MotionList{T<:Real}
     motions::Vector{<:Motion{T}}
 end
 
+# NOTE: this constructor must be simplified once the Vector{<:Motion} approach is accomplished: 
+# https://github.com/JuliaHealth/KomaMRI.jl/issues/480
 """ Constructors """
-MotionList(motions::Motion...) = length([motions]) > 0 ? MotionList([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`"
-
-""" MotionList sub-group """
-function Base.getindex(mv::MotionList{T}, p) where {T<:Real}
-    motion_array_aux = Motion{T}[]
-    for m in mv.motions
-        m[p] !== nothing ? push!(motion_array_aux, m[p]) : nothing
-    end
-    return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion()
-end
-function Base.view(mv::MotionList{T}, p) where {T<:Real}
-    motion_array_aux = Motion{T}[]
-    for m in mv.motions
-        @view(m[p]) !== nothing ? push!(motion_array_aux, @view(m[p])) : nothing
-    end
-    return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion()
+function MotionList(motions::Motion...) 
+    if length(motions) == 0
+        return NoMotion()
+    elseif length(motions) == 1
+        return motions[1]
+    else
+        return MotionList([motions...])
+    end 
 end
 
-""" Addition of MotionLists """
-function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1::Int, Ns2::Int) where {T<:Real}
+# NOTE: these vcat methods must be simplified once the Vector{<:Motion} approach is accomplished: 
+# 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}
     mv_aux = Motion{T}[]
     for m in m1.motions
         m_aux = copy(m)
@@ -69,7 +66,50 @@ function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1::Int, Ns2::Int) whe
         m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
         push!(mv_aux, m_aux)
     end
-    return MotionList(mv_aux)
+    return MotionList(mv_aux...)
+end
+# Motion + Motion
+function Base.vcat(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real}
+    mv_aux = Motion{T}[]
+    m_aux = copy(m1)
+    m_aux.spins = expand(m_aux.spins, Ns1)
+    push!(mv_aux, m_aux)
+    m_aux = copy(m2)
+    m_aux.spins = expand(m_aux.spins, Ns2)
+    m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
+    push!(mv_aux, m_aux)
+    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}
+    mv_aux = Motion{T}[]
+    m_aux = copy(m1)
+    m_aux.spins = expand(m_aux.spins, Ns1)
+    push!(mv_aux, m_aux)
+    for m in m2.motions
+        m_aux = copy(m)
+        m_aux.spins = expand(m_aux.spins, Ns2)
+        m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
+        push!(mv_aux, m_aux)
+    end
+    return MotionList(mv_aux...)
+end
+
+""" MotionList sub-group """
+function Base.getindex(mv::MotionList{T}, p) where {T<:Real}
+    motion_array_aux = Motion{T}[]
+    for m in mv.motions
+        m[p] isa NoMotion ? nothing : push!(motion_array_aux, m[p])
+    end
+    return MotionList(motion_array_aux...)
+end
+function Base.view(mv::MotionList{T}, p) where {T<:Real}
+    motion_array_aux = Motion{T}[]
+    for m in mv.motions
+        @view(m[p]) isa NoMotion ? nothing : push!(motion_array_aux, @view(m[p]))
+    end
+    return MotionList(motion_array_aux...)
 end
 
 """ Compare two MotionLists """
@@ -162,3 +202,9 @@ function sort_motions!(m::MotionList)
     sort!(m.motions; by=m -> times(m)[1])
     return nothing
 end
+
+function add_jump_times!(t, ml::MotionList)
+    for m in ml.motions
+        add_jump_times!(t, m)
+    end
+end

KomaMRIBase/src/motion/nomotion/NoMotion.jl → KomaMRIBase/src/motion/NoMotion.jl

@@ -17,7 +17,9 @@ Base.getindex(mv::NoMotion, p) = mv
 Base.view(mv::NoMotion, p)     = mv
 
 """ Addition of NoMotions """
+# NoMotion + NoMotion
 Base.vcat(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}
     mv_aux = Motion{T}[]
@@ -29,6 +31,14 @@ function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T}
     end
     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}
+    m_aux = copy(m2)
+    m_aux.spins = expand(m_aux.spins, Ns2)
+    m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
+    return m_aux
+end
 
 """ Compare two NoMotions """
 Base.:(==)(m1::NoMotion, m2::NoMotion) = true
@@ -39,3 +49,4 @@ function get_spin_coords(
 ) where {T<:Real}
     return x, y, z
 end
+add_jump_times!(t, ::NoMotion) = nothing

KomaMRIBase/src/motion/motionlist/SpinSpan.jl → KomaMRIBase/src/motion/SpinSpan.jl

@@ -69,3 +69,5 @@ Base.length(sr::SpinRange) = length(sr.range)
 get_indexing_range(spins::SpinRange) = spins.range
 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)

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

@@ -34,4 +34,9 @@ julia> flowpath = FlowPath(
     dy::AbstractArray{T}
     dz::AbstractArray{T}
     spin_reset::AbstractArray{Bool}
+end
+
+function add_jump_times!(t, a::FlowPath, time_span::AbstractTimeSpan) 
+    jump_times = (times(time_span)[end] - times(time_span)[1])/(size(a.spin_reset)[2]-1) * (getindex.(findall(a.spin_reset .== 1), 2) .- 1) .- 1e-6
+    append!(t, jump_times)
 end

KomaMRIBase/src/timing/TimeStepCalculation.jl

@@ -81,13 +81,15 @@ This function returns non-uniform time points that are relevant in the sequence
     samples for RF excitation (by nominal we mean that the time separation should be at most
     `Δt_rf` when the samples are regarded by [`KomaMRI.is_RF_on`](@ref), otherwise the time
     points are not necessary and the separation will be bigger)
+- `motion`: (`::Union{NoMotion, Motion, MotionList}`, `=NoMotion()`) phantom motion,
+	from which it may be necessary to extract key time points relevant for the simulation
 
 # Returns
 - `t`: (`::Vector{Float64}`, `[s]`) time array with non-uniform time values
 - `Δt`: (`::Vector{Float64}`, `[s]`) delta time array with the separation between two
     adjacent time points of the `t` time array
 """
-function get_variable_times(seq; Δt=1e-3, Δt_rf=1e-5)
+function get_variable_times(seq; Δt=1e-3, Δt_rf=1e-5, motion=NoMotion())
 	t = Float64[]
 	ϵ = MIN_RISE_TIME # Small Float64
 	T0 = get_block_start_times(seq)
@@ -120,6 +122,7 @@ function get_variable_times(seq; Δt=1e-3, Δt_rf=1e-5)
         end
         append!(t, t_block)
 	end
+	add_jump_times!(t, motion)
 	# Removing repeated points
 	sort!(unique!(t))
 	# Fixes a problem with ADC at the start and end of the seq