Fix flow-related functions

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/Motion.jl

@@ -202,4 +202,6 @@ 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.jl

@@ -1,6 +1,6 @@
-include("Action.jl")
 include("SpinSpan.jl")
 include("TimeSpan.jl")
+include("Action.jl")
 include("Motion.jl")
 
 """
@@ -202,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.jl

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

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

KomaMRICore/src/simulation/Flow.jl

@@ -81,13 +81,7 @@ function replace_view(replace_by, idx)
    return replace_by
 end
 
-function get_mask(spin_reset, t::Real)
+function get_mask(spin_reset, t)
    itp  = KomaMRIBase.interpolate(spin_reset, KomaMRIBase.Gridded(KomaMRIBase.Constant{KomaMRIBase.Next}()), Val(size(spin_reset, 1)), t)
    return KomaMRIBase.resample(itp, t)
-end
-function get_mask(spin_reset, t::AbstractArray)
-   itp  = KomaMRIBase.interpolate(spin_reset, KomaMRIBase.Gridded(KomaMRIBase.Constant{KomaMRIBase.Next}()), Val(size(spin_reset, 1)), t)
-   mask = KomaMRIBase.resample(itp, t)
-   mask .= (cumsum(mask; dims=2) .== 1)
-   return mask
 end

KomaMRICore/src/simulation/SimulatorCore.jl

@@ -336,7 +336,7 @@ function simulate(
         """ maxlog=1
     end
     # Simulation init
-    seqd = discretize(seq; sampling_params=sim_params) # Sampling of Sequence waveforms
+    seqd = discretize(seq; sampling_params=sim_params, motion=obj.motion) # Sampling of Sequence waveforms
     parts, excitation_bool = get_sim_ranges(seqd; Nblocks=sim_params["Nblocks"]) # Generating simulation blocks
     t_sim_parts = [seqd.t[p[1]] for p in parts]
     append!(t_sim_parts, seqd.t[end])