New Motion

KomaMRIBase/src/KomaMRIBase.jl

@@ -31,7 +31,6 @@ include("datatypes/Phantom.jl")
 include("datatypes/simulation/DiscreteSequence.jl")
 include("timing/TimeStepCalculation.jl")
 include("timing/TrapezoidalIntegration.jl")
-include("timing/UnitTime.jl")
 
 # Main
 export γ # gyro-magnetic ratio [Hz/T]
@@ -47,12 +46,11 @@ export kfoldperm, trapz, cumtrapz
 # Phantom
 export brain_phantom2D, brain_phantom3D, pelvis_phantom2D, heart_phantom
 # Motion
-export MotionModel
-export NoMotion, SimpleMotion, ArbitraryMotion
-export SimpleMotionType
-export Translation, Rotation, HeartBeat
-export PeriodicTranslation, PeriodicRotation, PeriodicHeartBeat
-export get_spin_coords
+export AbstractMotion, Motion, MotionVector, NoMotion
+export SimpleMotion, ArbitraryMotion
+export Translation, Rotation, HeartBeat, Trajectory, FlowTrajectory
+export TimeScale, TimeRange, Periodic
+export sort_motions!, get_spin_coords
 # Secondary
 export get_kspace, rotx, roty, rotz
 # Additionals

KomaMRIBase/src/datatypes/Phantom.jl

@@ -1,9 +1,9 @@
-abstract type MotionModel{T<:Real} end
-
-#Motion models:
-include("phantom/motion/SimpleMotion.jl")
-include("phantom/motion/ArbitraryMotion.jl")
-include("phantom/motion/NoMotion.jl")
+# TimeScale:
+include("../timing/TimeScale.jl")
+# Motion:
+abstract type AbstractMotion{T<:Real} end
+include("phantom/Motion.jl")
+include("phantom/NoMotion.jl")
 
 """
  obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion)
@@ -54,7 +54,7 @@
  Dθ::AbstractVector{T} = zeros(eltype(x), size(x))
  #Diff::Vector{DiffusionModel} #Diffusion map
  #Motion
- motion::MotionModel{T} = NoMotion{eltype(x)}()
+ motion::AbstractMotion{T} = NoMotion{eltype(x)}()
 end
 
 """Size and length of a phantom"""
@@ -67,18 +67,29 @@
 Base.getindex(x::Phantom, i::Integer) = x[i:i]
 
 """Compare two phantoms"""
-Base.:(==)(obj1::Phantom, obj2::Phantom) = reduce(
- &,
- [getfield(obj1, field) == getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))],
-)
-Base.:(≈)(obj1::Phantom, obj2::Phantom) = reduce(&, [getfield(obj1, field) ≈ getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))])
-Base.:(==)(m1::MotionModel, m2::MotionModel) = false
-Base.:(≈)(m1::MotionModel, m2::MotionModel) = false
+function Base.:(==)(obj1::Phantom, obj2::Phantom)
+ return reduce(
+ &,
+ [
+ getfield(obj1, field) == getfield(obj2, field) for
+ field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
+ ],
+ )
+end
+function Base.:(≈)(obj1::Phantom, obj2::Phantom)
+ return reduce(
+ &,
+ [
+ getfield(obj1, field) ≈ getfield(obj2, field) for
+ field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
+ ],
+ )
+end
 
 """Separate object spins in a sub-group"""
 Base.getindex(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin
  fields = []
- for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
+ for field in Iterators.filter(x -> x != :name, fieldnames(Phantom))
  push!(fields, (field, getfield(obj, field)[p]))
  end
  return Phantom(; name=obj.name, fields...)
@@ -87,21 +98,25 @@
 """Separate object spins in a sub-group (lightweigth)."""
 Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin
  fields = []
- for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
+ for field in Iterators.filter(x -> x != :name, fieldnames(Phantom))
  push!(fields, (field, @view(getfield(obj, field)[p])))
  end
  return Phantom(; name=obj.name, fields...)
 end
 
 """Addition of phantoms"""
 +(obj1::Phantom, obj2::Phantom) = begin
+ Nmaxchars = 50
+ name = first(obj1.name * "+" * obj2.name, Nmaxchars)
  fields = []
- for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
+ for field in Iterators.filter(x -> x != :name, fieldnames(Phantom))
  push!(fields, (field, [getfield(obj1, field); getfield(obj2, field)]))
  end
- Nmaxchars = 50
- name = first(obj1.name * "+" * obj2.name, Nmaxchars)
- return Phantom(; name=name, fields...)
+ return Phantom(; 
+ name=name, 
+ fields..., 
+ motion=vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))
+ )
 end
 
 """Scalar multiplication of a phantom"""
@@ -177,16 +192,18 @@
  Dλ1=Dλ1[ρ .!= 0],
  Dλ2=Dλ2[ρ .!= 0],
  Dθ=Dθ[ρ .!= 0],
- motion=SimpleMotion(
- PeriodicHeartBeat(;
- period=period,
- asymmetry=asymmetry,
+ motion=MotionVector(
+ HeartBeat(;
+ times=Periodic(; period=period, asymmetry=asymmetry),
  circumferential_strain=circumferential_strain,
  radial_strain=radial_strain,
  longitudinal_strain=0.0,
  ),
- PeriodicRotation(;
- period=period, asymmetry=asymmetry, yaw=rotation_angle, pitch=0.0, roll=0.0
+ Rotation(;
+ times=Periodic(; period=period, asymmetry=asymmetry),
+ yaw=rotation_angle,
+ pitch=0.0,
+ roll=0.0,
  ),
  ),
  )

KomaMRIBase/src/datatypes/phantom/Motion.jl

@@ -0,0 +1,120 @@
+abstract type Motion{T<:Real} end
+
+is_composable(m::Motion) = false
+
+struct MotionVector{T<:Real} <: AbstractMotion{T}
+ motions::Vector{<:Motion{T}}
+end
+
+MotionVector(motions...) = length([motions]) > 0 ? MotionVector([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`"
+
+include("motion/SimpleMotion.jl")
+include("motion/ArbitraryMotion.jl")
+
+Base.getindex(mv::MotionVector, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionVector(getindex.(mv.motions, Ref(p))) 
+Base.view(mv::MotionVector, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionVector(view.(mv.motions, Ref(p)))
+
+""" Addition of MotionVectors """
+function Base.vcat(m1::MotionVector{T}, m2::MotionVector{T}, Ns1::Int, Ns2::Int) where {T<:Real}
+ mv1 = m1.motions
+ mv1_aux = Motion{T}[]
+ for i in 1:length(mv1)
+ if typeof(mv1[i]) <: ArbitraryMotion
+ zeros1 = similar(mv1[i].dx, Ns2, size(mv1[i].dx, 2))
+ zeros1 .= zero(T)
+ push!(mv1_aux, typeof(mv1[i])(mv1[i].times, [[getfield(mv1[i], d); zeros1] for d in filter(x -> x != :times, fieldnames(typeof(mv1[i])))]...)) 
+ else
+ push!(mv1_aux, mv1[i])
+ end
+ end
+ mv2 = m2.motions
+ mv2_aux = Motion{T}[]
+ for i in 1:length(mv2)
+ if typeof(mv2[i]) <: ArbitraryMotion
+ zeros2 = similar(mv2[i].dx, Ns1, size(mv2[i].dx, 2))
+ zeros2 .= zero(T)
+ push!(mv2_aux, typeof(mv2[i])(mv2[i].times, [[zeros2; getfield(mv2[i], d)] for d in filter(x -> x != :times, fieldnames(typeof(mv2[i])))]...))
+ else
+ push!(mv2_aux, mv2[i])
+ end
+ end
+ return MotionVector([mv1_aux; mv2_aux])
+end
+
+""" Compare two motion vectors """
+function Base.:(==)(mv1::MotionVector{T}, mv2::MotionVector{T}) where {T<:Real}
+ sort_motions!(mv1)
+ sort_motions!(mv2)
+ return reduce(&, mv1.motions .== mv2.motions)
+end
+function Base.:(≈)(mv1::MotionVector{T}, mv2::MotionVector{T}) where {T<:Real} 
+ sort_motions!(mv1)
+ sort_motions!(mv2)
+ return reduce(&, mv1.motions .≈ mv2.motions)
+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_{\text{spins}}`` rows and `length(t)` columns.
+
+# Arguments
+- `motion`: (`::Vector{<:Motion{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`: (`::AbstractArray{T<:Real}`) horizontal array of time instants
+
+# Returns
+- `x, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time
+"""
+function get_spin_coords(
+ mv::MotionVector{T},
+ x::AbstractVector{T},
+ y::AbstractVector{T},
+ z::AbstractVector{T},
+ t::AbstractArray{T}
+) where {T<:Real}
+ # Buffers for positions:
+ xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t
+ # Buffers for displacements:
+ ux, uy, uz = similar(xt), similar(yt), similar(zt) 
+
+ # Composable motions: they need to be run sequentially
+ for m in Iterators.filter(is_composable, mv.motions)
+ displacement_x!(ux, m, xt, yt, zt, t)
+ displacement_y!(uy, m, xt, yt, zt, t)
+ displacement_z!(uz, m, xt, yt, zt, t)
+ xt .+= ux
+ yt .+= uy
+ zt .+= uz
+ end
+ # Additive motions: these motions can be run in parallel
+ for m in Iterators.filter(!is_composable, mv.motions)
+ displacement_x!(ux, m, x, y, z, t)
+ displacement_y!(uy, m, x, y, z, t)
+ displacement_z!(uz, m, x, y, z, t)
+ xt .+= ux
+ yt .+= uy
+ zt .+= uz
+ end
+ return xt, yt, zt
+end
+
+"""
+ times = times(motion)
+"""
+times(m::Motion) = times(m.times)
+times(mv::MotionVector{T}) where {T<:Real} = begin
+ nodes = reduce(vcat, [times(m) for m in mv.motions]; init=[zero(T)])
+ return unique(sort(nodes))
+end
+
+"""
+ sort_motions!
+"""
+function sort_motions!(mv::MotionVector{T}) where {T<:Real}
+ sort!(mv.motions; by=m -> times(m)[1])
+ return nothing
+end

KomaMRIBase/src/datatypes/phantom/NoMotion.jl

@@ -0,0 +1,31 @@
+struct NoMotion{T<:Real} <: AbstractMotion{T} end
+
+Base.getindex(mv::NoMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = mv 
+Base.view(mv::NoMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = mv
+
+""" Addition of NoMotions """
+Base.vcat(m1::NoMotion{T}, m2::AbstractMotion{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m2
+Base.vcat(m1::AbstractMotion{T}, m2::NoMotion{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m1
+
+Base.:(==)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true
+Base.:(≈)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true
+
+function get_spin_coords(
+ mv::NoMotion{T},
+ x::AbstractVector{T},
+ y::AbstractVector{T},
+ z::AbstractVector{T},
+ t::AbstractArray{T}
+) where {T<:Real}
+ return x, y, z
+end
+
+"""
+ times = times(motion)
+"""
+times(mv::NoMotion{T}) where {T<:Real} = [zero(T)]
+
+"""
+ sort_motions!
+"""
+sort_motions!(mv::NoMotion) = nothing