fix knot collapse at large absolute time

Author: Arthur Allilaire

KomaMRIBase/src/datatypes/Sequence.jl

@@ -703,8 +703,8 @@ function get_samples(seq::Sequence, range; events=[:rf, :gr, :adc], freq_in_phas
     fill_if_empty(x) = isempty(x.t) && length(range) == length(seq) ? merge(x, (t=[0.0; dur(seq)], A=zeros(eltype(x.A), 2))) : x
     # RF
     if :rf in events
-        t_rf = reduce(vcat, [T0[i] .+ times(seq.RF[1,i], :A)   for i in range])
-        t_Δf = reduce(vcat, [T0[i] .+ times(seq.RF[1,i], :Δf)  for i in range])
+        t_rf = reduce(vcat, [_reseparate_closing_knot!(T0[i] .+ times(seq.RF[1,i], :A))  for i in range])
+        t_Δf = reduce(vcat, [_reseparate_closing_knot!(T0[i] .+ times(seq.RF[1,i], :Δf)) for i in range])
         A_rf = reduce(vcat, [ampls(seq.RF[1,i]; freq_in_phase) for i in range])
         A_Δf = reduce(vcat, [freqs(seq.RF[1,i])                for i in range])
         A_ψ  = reduce(vcat, [rf_frame_phase(seq.RF[1,i])       for i in range])
@@ -719,9 +719,9 @@ function get_samples(seq::Sequence, range; events=[:rf, :gr, :adc], freq_in_phas
         gx = [_gradient_interpolation_samples(seq.GR[1,i]) for i in range]
         gy = [_gradient_interpolation_samples(seq.GR[2,i]) for i in range]
         gz = [_gradient_interpolation_samples(seq.GR[3,i]) for i in range]
-        t_gx = reduce(vcat, [T0[i] .+ gx[j].t for (j, i) in enumerate(range)])
-        t_gy = reduce(vcat, [T0[i] .+ gy[j].t for (j, i) in enumerate(range)])
-        t_gz = reduce(vcat, [T0[i] .+ gz[j].t for (j, i) in enumerate(range)])
+        t_gx = reduce(vcat, [_strictly_increasing_knots!(_reseparate_closing_knot!(T0[i] .+ gx[j].t)) for (j, i) in enumerate(range)])
+        t_gy = reduce(vcat, [_strictly_increasing_knots!(_reseparate_closing_knot!(T0[i] .+ gy[j].t)) for (j, i) in enumerate(range)])
+        t_gz = reduce(vcat, [_strictly_increasing_knots!(_reseparate_closing_knot!(T0[i] .+ gz[j].t)) for (j, i) in enumerate(range)])
         A_gx = reduce(vcat, [g.A for g in gx])
         A_gy = reduce(vcat, [g.A for g in gy])
         A_gz = reduce(vcat, [g.A for g in gz])

KomaMRIBase/src/datatypes/sequence/Grad.jl

@@ -253,7 +253,10 @@ _amplitude_roundoff_tol(A) = 1000 * eps(float(maximum(abs, A)))
 
 function _strictly_increasing_knots!(t)
     for i in 2:length(t)
-        t[i] <= t[i - 1] && (t[i] = t[i - 1] + MIN_RISE_TIME)
+        if t[i] <= t[i - 1]
+            gap = max(MIN_RISE_TIME, eps(t[i - 1]))
+            t[i] = t[i - 1] + gap
+        end
     end
     return t
 end

KomaMRIBase/src/timing/KeyValuesCalculation.jl

@@ -105,6 +105,13 @@ function _separate_closing_knot!(t)
     return t
 end
 
+# Idempotent re-separation of knots that could collapse after rebasing to absolute time
+function _reseparate_closing_knot!(t)
+    length(t) >= 2 || return t
+    t[end - 1] = min(t[end - 1], t[end] - max(MIN_RISE_TIME, eps(t[end])))
+    return t
+end
+
 _gradient_times(gr::TrapezoidalGrad) = cumsum([gr.delay; gr.rise; gr.T; gr.fall])
 
 times(gr::TrapezoidalGrad) =

KomaMRIBase/src/timing/TimeStepCalculation.jl

@@ -114,7 +114,7 @@ function get_variable_times(seq; Δt=1e-3, Δt_rf=1e-5, motion=NoMotion())
 			if is_Gy_on(s) append!(active_gradients, s.GR.y) end
 			if is_Gz_on(s) append!(active_gradients, s.GR.z) end
 			for y = active_gradients
-				ts = times(y) .+ t0
+				ts = _reseparate_closing_knot!(times(y) .+ t0)
 				taux = points_from_key_times([next_time(ts[1]); ts; prev_time(ts[end])]; dt=Δt)
                 append!(t_block, taux)
 			end

KomaMRIBase/test/runtests.jl

@@ -813,6 +813,25 @@ using TestItems, TestItemRunner
         @test all(diff(seqd.t) .> 0)
     end
 
+    @testset "closing knot survives FP rebasing at large t0" begin
+        for t0 in (270.0, 1500.0), event in (RF(1e-5, 1e-3), Grad(1e-3, 1e-3))
+            ts = KomaMRIBase._reseparate_closing_knot!(t0 .+ KomaMRIBase.times(event))
+            @test ts[end-1] < ts[end]
+        end
+    end
+
+    @testset "get_samples and get_variable_times keep knots separated at large t0" begin
+        T, offset = 1e-3, 300.0
+        rf_seq = Sequence(); rf_seq += Delay(offset); @addblock rf_seq += RF(1e-5, T)
+        gr_seq = Sequence(); gr_seq += Delay(offset); @addblock gr_seq += Grad(1e-3, T)
+        ref_seq = Sequence(); ref_seq += Delay(0.5); @addblock ref_seq += Grad(1e-3, T)
+        rf_t = KomaMRIBase.get_samples(rf_seq).rf.t
+        @test rf_t[end-1] < rf_t[end]
+        @test all(diff(KomaMRIBase.get_samples(gr_seq).gx.t) .> 0)
+        @test length(KomaMRIBase.get_variable_times(gr_seq)[1]) ==
+              length(KomaMRIBase.get_variable_times(ref_seq)[1])
+    end
+
      @testset "SequenceFunctions" begin
         seq = PulseDesigner.EPI_example()
         t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1)