Optimize ArbitraryMotion (continuation)

KomaMRIBase/src/KomaMRIBase.jl

@@ -7,7 +7,7 @@ using Reexport
 #Datatypes
 using Parameters
 #Simulation
-using Interpolations
+@reexport using Interpolations
 #Reconstruction
 using MRIBase
 @reexport using MRIBase:
@@ -31,6 +31,7 @@ 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]
@@ -51,8 +52,8 @@ export NoMotion, SimpleMotion, ArbitraryMotion
 export SimpleMotionType
 export Translation, Rotation, HeartBeat
 export PeriodicTranslation, PeriodicRotation, PeriodicHeartBeat
-export get_spin_coords, sort_motions!
-export LinearInterpolator
+export get_spin_coords, initialize_motion!
+export Interpolator
 # Secondary
 export get_kspace, rotx, roty, rotz
 # Additionals

KomaMRIBase/src/datatypes/Phantom.jl

@@ -69,10 +69,7 @@ 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))
-    ],
+    [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
@@ -117,7 +114,7 @@ function get_dims(obj::Phantom)
     return dims
 end
 
-function sort_motions!(motion::MotionModel) 
+function initialize_motion!(motion::MotionModel) 
     return nothing
 end
 

KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl

@@ -2,13 +2,30 @@
 #       Interpolator{T,Degree,ETPType}, 
 #           Degree = Linear,Cubic.... 
 #           ETPType = Periodic, Flat...
-const LinearInterpolator = Interpolations.Extrapolation{
-    T,
-    1,
-    Interpolations.GriddedInterpolation{T,1,V,Gridded{Linear{Throw{OnGrid}}},Tuple{V}},
-    Gridded{Linear{Throw{OnGrid}}},
-    Periodic{Nothing},
-} where {T<:Real,V<:AbstractVector{T}}
+
+const Interpolator = Interpolations.GriddedInterpolation{
+    T,N,V,Itp,K
+} where {
+    T<:Real,
+    N,
+    V<:AbstractArray{T},
+    Itp<:Tuple{ Vararg{ Union{Interpolations.Gridded{Linear{Throw{OnGrid}}},Interpolations.NoInterp} } },
+    K<:Tuple{Vararg{AbstractVector{T}}},
+}
+
+function GriddedInterpolation(
+    nodes::NType,
+    A::AType
+) where {T<:Real, AType<:AbstractArray{T}, NType<:Tuple{Vararg{AbstractVector{T}}}}
+    Ns, _ = size(A)
+    if Ns > 1
+        ITPType = Tuple{NoInterp, Gridded{Linear{Throw{OnGrid}}}}
+        return Interpolations.GriddedInterpolation{T, 2, typeof(A), ITPType, typeof(nodes)}(nodes, A, (NoInterp(), Gridded(Linear())))
+    else
+        ITPType = Tuple{Gridded{Linear{Throw{OnGrid}}}}
+        return Interpolations.GriddedInterpolation{T, 1, typeof(A[:]), ITPType, typeof((nodes[2], ))}((nodes[2], ), A[:], (Gridded(Linear()), ))
+    end
+end
 
 """
     motion = ArbitraryMotion(period_durations, dx, dy, dz)
@@ -43,106 +60,87 @@ julia> motion = ArbitraryMotion(
         )
 ```
 """
-struct ArbitraryMotion{T<:Real,V<:AbstractVector{T}} <: MotionModel{T}
-    period_durations::Vector{T}
-    dx::Array{T,2}
-    dy::Array{T,2}
-    dz::Array{T,2}
-    ux::Vector{LinearInterpolator{T,V}}
-    uy::Vector{LinearInterpolator{T,V}}
-    uz::Vector{LinearInterpolator{T,V}}
-end
-
-function ArbitraryMotion(
-    period_durations::AbstractVector{T},
-    dx::AbstractArray{T,2},
-    dy::AbstractArray{T,2},
-    dz::AbstractArray{T,2},
-) where {T<:Real}
-    @warn "Note that ArbitraryMotion is under development so it is not optimized so far" maxlog = 1
-    Ns = size(dx)[1]
-    num_pieces = size(dx)[2] + 1
-    limits = times(period_durations, num_pieces)
-
-    #! format: off
-    Δ = zeros(Ns,length(limits),4)
-    Δ[:,:,1] = hcat(repeat(hcat(zeros(Ns,1),dx),1,length(period_durations)),zeros(Ns,1))
-    Δ[:,:,2] = hcat(repeat(hcat(zeros(Ns,1),dy),1,length(period_durations)),zeros(Ns,1))
-    Δ[:,:,3] = hcat(repeat(hcat(zeros(Ns,1),dz),1,length(period_durations)),zeros(Ns,1))
-
-    etpx = [extrapolate(interpolate((limits,), Δ[i,:,1], Gridded(Linear())), Periodic()) for i in 1:Ns]
-    etpy = [extrapolate(interpolate((limits,), Δ[i,:,2], Gridded(Linear())), Periodic()) for i in 1:Ns]
-    etpz = [extrapolate(interpolate((limits,), Δ[i,:,3], Gridded(Linear())), Periodic()) for i in 1:Ns]
-    #! format: on
-
-    return ArbitraryMotion(period_durations, dx, dy, dz, etpx, etpy, etpz)
+struct ArbitraryMotion{T} <: MotionModel{T}
+    t_start::T
+    t_end::T
+    dx::AbstractArray{T}
+    dy::AbstractArray{T}
+    dz::AbstractArray{T}
 end
 
 function Base.getindex(
     motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon}
 )
-    fields = []
-    for field in fieldnames(ArbitraryMotion)
-        if field in (:dx, :dy, :dz)
-            push!(fields, getfield(motion, field)[p, :])
-        elseif field in (:ux, :uy, :uz)
-            push!(fields, getfield(motion, field)[p])
-        else
-            push!(fields, getfield(motion, field))
-        end
-    end
-    return ArbitraryMotion(fields...)
+    return ArbitraryMotion(motion.t_start, motion.t_end, motion.dx[p,:], motion.dy[p,:], motion.dz[p,:])
 end
 
 Base.:(==)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(ArbitraryMotion)])
 Base.:(≈)(m1::ArbitraryMotion, m2::ArbitraryMotion)  = reduce(&, [getfield(m1, field) ≈ getfield(m2, field) for field in fieldnames(ArbitraryMotion)])
 
 function Base.vcat(m1::ArbitraryMotion, m2::ArbitraryMotion)
     fields = []
-    @assert m1.period_durations == m2.period_durations "period_durations of both ArbitraryMotions must be the same"
-    for field in
-        Iterators.filter(x -> !(x == :period_durations), fieldnames(ArbitraryMotion))
+    @assert (m1.t_start == m2.t_start) && (m1.t_end == m2.t_end) "starting and ending times must be the same"
+    for field in (:dx, :dy, :dz)
         push!(fields, [getfield(m1, field); getfield(m2, field)])
     end
-    return ArbitraryMotion(m1.period_durations, fields...)
+    return ArbitraryMotion(m1.t_start, m1.t_end, fields...)
 end
 
 """
     limits = times(obj.motion)
 """
 function times(motion::ArbitraryMotion)
-    period_durations = motion.period_durations
-    num_pieces = size(motion.dx)[2] + 1
-    return times(period_durations, num_pieces)
+    return collect(range(motion.t_start, motion.t_end, length=size(motion.dx)[2]))
 end
 
-function times(period_durations::AbstractVector, num_pieces::Int)
-    # Pre-allocating memory
-    limits = zeros(eltype(period_durations), num_pieces * length(period_durations) + 1)
-
-    idx = 1
-    for i in 1:length(period_durations)
-        segment_increment = period_durations[i] / num_pieces
-        cumulative_sum = limits[idx]  # Start from the last computed value in limits
-        for j in 1:num_pieces
-            cumulative_sum += segment_increment
-            limits[idx + 1] = cumulative_sum
-            idx += 1
-        end
+
+function get_itp_functions(
+    dx::AbstractArray{T},
+    dy::AbstractArray{T},
+    dz::AbstractArray{T}
+    ) where {T<:Real}
+    Ns, Nt = size(dx)
+
+    t = similar(dx, Nt)
+    t .= 0:(Nt-1)
+
+    id = similar(dx, Ns)
+    id .= 1:Ns
+
+    itpx = GriddedInterpolation((id, t ./ (Nt-1)), dx)
+    itpy = GriddedInterpolation((id, t ./ (Nt-1)), dy)
+    itpz = GriddedInterpolation((id, t ./ (Nt-1)), dz)
+
+    return itpx, itpy, itpz
+end
+
+function get_itp_results(
+    itpx::Interpolator{T}, 
+    itpy::Interpolator{T}, 
+    itpz::Interpolator{T}, 
+    t::AbstractArray{T}
+) where {T<:Real}
+    Ns = ndims(itpx.coefs) == 1 ? 1 : size(itpx.coefs,1)
+    if Ns > 1
+        id = similar(t, Ns)
+        id .= 1:Ns
+        # Grid
+        idx = 1*id .+ 0*t # spin id
+        t   = 0*id .+ 1*t # time instants
+        return itpx.(idx, t), itpy.(idx, t), itpz.(idx, t)
+    else
+        return itpx.(t), itpy.(t), itpz.(t)
     end
-    return limits
 end
 
-# TODO: Calculate interpolation functions "on the fly"
 function get_spin_coords(
     motion::ArbitraryMotion{T},
     x::AbstractVector{T},
     y::AbstractVector{T},
     z::AbstractVector{T},
-    t::AbstractArray{T},
+    t::AbstractArray{T}
 ) where {T<:Real}
-    xt = x .+ reduce(vcat, [etp.(t) for etp in motion.ux])
-    yt = y .+ reduce(vcat, [etp.(t) for etp in motion.uy])
-    zt = z .+ reduce(vcat, [etp.(t) for etp in motion.uz])
-    return xt, yt, zt
-end
+    itp = get_itp_functions(motion.dx, motion.dy, motion.dz)
+    ux, uy, uz = get_itp_results(itp..., unit_time(t, motion.t_start, motion.t_end))
+    return x .+ ux, y .+ uy, z .+ uz
+end

KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl

@@ -18,7 +18,7 @@ function get_spin_coords(
     x::AbstractVector{T},
     y::AbstractVector{T},
     z::AbstractVector{T},
-    t::AbstractArray{T},
+    t::AbstractArray{T}
 ) where {T<:Real}
     return x, y, z
 end

KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl

@@ -71,7 +71,7 @@
     x::AbstractVector{T},
     y::AbstractVector{T},
     z::AbstractVector{T},
-    t::AbstractArray{T},
+    t::AbstractArray{T}
 ) where {T<:Real}
     xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t
     # Composable motions: they need to be run sequentially
@@ -96,8 +96,8 @@
     return nodes
 end
 
-function sort_motions!(motion::SimpleMotion)
-    return sort!(motion.types; by=m -> times(m)[1])
+function initialize_motion!(motion::SimpleMotion)
+    sort!(motion.types; by=m -> times(m)[1])
 end
 
 # --------- Simple Motion Types: -------------
@@ -106,28 +106,8 @@
 include("simplemotion/Rotation.jl")
 include("simplemotion/HeartBeat.jl")
 
-function unit_time(t::AbstractArray{T}, t_start::T, t_end::T) where {T<:Real}
-    if t_start == t_end
-        return (t .>= t_start) .* one(T) # The problem with this is that it returns a BitVector (type stability issues)
-    else
-        return min.(max.((t .- t_start) ./ (t_end - t_start), zero(T)), one(T))
-    end
-end
-
 # Periodic types: defined by the period, the temporal symmetry and a displacement (amplitude)
 include("simplemotion/PeriodicTranslation.jl")
 include("simplemotion/PeriodicRotation.jl")
 include("simplemotion/PeriodicHeartBeat.jl")
 
-function unit_time_triangular(t, period, asymmetry)
-    t_rise = period * asymmetry
-    t_fall = period * (1 - asymmetry)
-    t_relative = mod.(t, period)
-    t_unit =
-        ifelse.(
-            t_relative .< t_rise,
-            t_relative ./ t_rise,
-            1 .- (t_relative .- t_rise) ./ t_fall,
-        )
-    return t_unit
-end

KomaMRIBase/src/timing/UnitTime.jl

@@ -0,0 +1,30 @@
+"""
+    t_unit = unit_time(t, t_start, t_end)
+
+For non-periodic motions
+"""
+function unit_time(t::AbstractArray{T}, t_start::T, t_end::T) where {T<:Real}
+    if t_start == t_end
+        return (t .>= t_start) .* one(T)
+    else
+        return min.(max.((t .- t_start) ./ (t_end - t_start), zero(T)), one(T))
+    end
+end
+
+"""
+    t_unit = unit_time_triangular(t, period, asymmetry)
+
+For periodic motions
+"""
+function unit_time_triangular(t::AbstractArray{T}, period::T, asymmetry::T) where {T<:Real}
+    t_rise = period * asymmetry
+    t_fall = period * (1 - asymmetry)
+    t_relative = mod.(t, period)
+    t_unit =
+        ifelse.(
+            t_relative .< t_rise,
+            t_relative ./ t_rise,
+            1 .- (t_relative .- t_rise) ./ t_fall,
+        )
+    return t_unit
+end