Merge pull request #517 from JuliaHealth/easy-motion

Ease motion definition when only one motion is needed

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

@@ -178,11 +178,26 @@ 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

KomaMRIBase/src/motion/motionlist/MotionList.jl → KomaMRIBase/src/motion/MotionList.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 """

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