Incorrect gradient interpretation when waveforms do not end in zero

KomaMRIBase/src/datatypes/sequence/Grad.jl

@@ -74,6 +74,9 @@ mutable struct Grad
 	rise::Real
 	fall::Real
 	delay::Real
+    first
+    last
+    time_shape_id
     function Grad(A, T, rise, fall, delay)
 		all(T .< 0) || rise < 0 || fall < 0 || delay < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, fall, delay)
     end

KomaMRIBase/src/timing/KeyValuesCalculation.jl

@@ -32,6 +32,14 @@ get_theo_A(g::Grad; off_val=0) = begin
 	else
 		aux = [off_val; 0; A; 0]
 	end
+    # To remove abrut changes in the plots
+    if g.rise == 0
+        aux[2] = 0      # zero
+    end
+    if g.fall == 0
+        aux[end] = aux[end-1]   # keep the last value
+    end
+    # If no signal, then don't show samples
 	if sum(abs.(g.A)) == 0
 		aux *= off_val
 	end

KomaMRIFiles/src/Sequence/Pulseq.jl

@@ -434,9 +434,39 @@ function read_seq(filename)
     for i = 1:length(blockEvents)
         seq += get_block(obj, i)
     end
-    # Final details
     # Remove dummy seq block at the start, Issue #203
     seq = seq[2:end]
+
+    # Add first and last Pulseq points
+    for gi in 1:3
+        grad_prev_last = 0
+        for bi in 1:length(seq)
+            gr = seq.GR[gi, bi]
+            if sum(abs.(gr.A)) == 0     # this is for no-gradient case
+                grad_prev_last = 0
+            else
+                if length(gr.A) > 1   # this is for the uniformly-shaped or time-shaped case
+                    if gr.delay > 0
+                        grad_prev_last = 0
+                    end
+                    seq.GR[gi, bi].first = grad_prev_last
+                    if length(gr.T) > 1 || (isdefined(gr, :time_shape_id) && gr.time_shape_id != 0)    # this is for time-shaped case (I assume it is the extended trapezoid case)
+                        seq.GR[gi, bi].last = gr.A[end]   #I need to check this or [end-1]
+                    else
+                        odd_step1 = [seq.GR[gi, bi].first; 2 * gr.A]
+                        odd_step2 = odd_step1 .* (mod.(1:length(odd_step1), 2) * 2 .- 1)
+                        waveform_odd_rest = cumsum(odd_step2) .* (mod.(1:length(odd_step2), 2) * 2 .- 1)
+                        seq.GR[gi, bi].last = waveform_odd_rest[end]
+                    end
+                    grad_prev_last = seq.GR[gi, bi].last
+                else    # this is for the trapedoid case
+                    grad_prev_last = 0
+                end
+            end
+        end
+    end
+
+    # Final details
     # Hack for including extension and triggers
     seq.DEF["additional_text"] = read_Extension(extensionLibrary, triggerLibrary) #Temporary hack
     seq.DEF = KomaMRIBase.recursive_merge(obj["definitions"], seq.DEF)
@@ -458,6 +488,8 @@ function read_seq(filename)
     return seq
 end
 
+
+
 #To Sequence
 """
     grad = read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
@@ -492,16 +524,24 @@ function read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
         delay =         g[4]
         #Amplitude
         gA = amplitude * decompress_shape(shapeLibrary[amp_shape_id]...)
+        if length(gA) == 1
+            gA = gA[1]
+        end
         Nrf = length(gA) - 1
         #Creating timings
         if time_shape_id == 0 #no time waveform
             gT = Nrf * Δt_gr
             G = Grad(gA, gT, Δt_gr/2, Δt_gr/2, delay)
+            #G = Grad(gA, gT, 0, 0, delay)
         else
             gt = decompress_shape(shapeLibrary[time_shape_id]...)
             gT = (gt[2:end] .- gt[1:end-1]) * Δt_gr
+            if length(gT) == 1
+                gT = gT[1]
+            end
             G = Grad(gA,gT,0,0,delay)
         end
+        G.time_shape_id = time_shape_id
     end
     G
 end
@@ -521,6 +561,11 @@ Reads the RF. It is used internally by [`get_block`](@ref).
 - `rf`: (`1x1 ::Matrix{RF}`) RF struct
 """
 function read_RF(rfLibrary, shapeLibrary, Δt_rf, i)
+
+    if isempty(rfLibrary) || i==0
+        return reshape([RF(0,0)], 1, 1)
+    end
+
     #Unpacking
     #(1)amplitude (2)mag_id (3)phase_id (4)time_shape_id (5)delay (6)freq (7)phase
     r = rfLibrary[i]["data"]
@@ -533,8 +578,8 @@ function read_RF(rfLibrary, shapeLibrary, Δt_rf, i)
     phase =         r[7]
     #Amplitude and phase waveforms
     if amplitude != 0 && mag_id != 0
-        rfA = decompress_shape(shapeLibrary[mag_id]...)[1:end-1]
-        rfϕ = decompress_shape(shapeLibrary[phase_id]...)[1:end-1]
+        rfA = decompress_shape(shapeLibrary[mag_id]...)#[1:end-1]   # Temporal change for reading all samples from .seq
+        rfϕ = decompress_shape(shapeLibrary[phase_id]...)#[1:end-1]     # Temporal change for reading all samples from .seq
         @assert all(rfϕ.>=0) "[RF id $i] Phase waveform rfϕ must have non-negative samples (1.>=rfϕ.>=0). "
         Nrf = shapeLibrary[mag_id][1] - 1
         rfAϕ = amplitude .* rfA .* exp.(1im*(2π*rfϕ .+ phase))
@@ -550,6 +595,7 @@ function read_RF(rfLibrary, shapeLibrary, Δt_rf, i)
     else
         rft = decompress_shape(shapeLibrary[time_shape_id]...)
         rfT = (rft[2:end] .- rft[1:end-1]) * Δt_rf
+        rfAϕ = rfAϕ[1:length(rfT)]  # Temporal change for RF-Pulses when reading all samples from .seq
     end
     R = reshape([RF(rfAϕ,rfT,freq,delay)],1,1)#[RF(rfAϕ,rfT,freq,delay);;]
     R
@@ -568,8 +614,14 @@ Reads the ADC. It is used internally by [`get_block`](@ref).
 - `adc`: (`1x1 ::Vector{ADC}`) ADC struct
 """
 function read_ADC(adcLibrary, i)
+
+    if isempty(adcLibrary) || i==0
+        return [ADC(0, 0)]
+    end
+
     #Unpacking
     #(1)num (2)dwell (3)delay (4)freq (5)phase
+    # It is needed a review for this, the dwell definition may cause problems
     if !isempty(adcLibrary) # Is this the best? maybe defining i=0 is better, it works with RFs(?)
         a = adcLibrary[i]["data"]
     else

KomaMRIFiles/test/Project.toml

@@ -1,4 +1,7 @@
 [deps]
+KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"

KomaMRIFiles/test/runtests.jl

@@ -49,3 +49,192 @@ using TestItems, TestItemRunner
         @test obj.name == "brain_mrilab.mat"
     end
 end
+
+@testitem "Pulseq_Read_Comparison" tags=[:files] begin
+    using KomaMRIBase, MAT, PrettyTables
+    TOLERANCE = 1e-6
+
+    # Auxiliar functions
+    function get_theo_t_aux(rf::RF)
+        Namp, Ntim = length(rf.A), length(rf.T)
+        if Namp == 1 && Ntim == 1
+            return KomaMRIBase.get_theo_t(rf; max_rf_samples=Inf)[2:end]
+        elseif Namp > 1 && Ntim == 1
+            amps = KomaMRIBase.get_theo_t(rf; max_rf_samples=Inf)[2:end]
+            return [amps[1]; amps[2:end-1][[i for i in 1:length(amps)-1 if i % 2 == 0]]; amps[end]]
+        end
+        return []#KomaMRIBase.get_theo_t(rf; max_rf_samples=Inf)[2:end]
+    end
+    function get_theo_A_aux(rf::RF)
+        Namp, Ntim = length(rf.A), length(rf.T)
+        if Namp == 1 && Ntim == 1
+            return γ*KomaMRIBase.get_theo_A(rf; off_val=NaN, max_rf_samples=Inf)[2:end]
+        elseif Namp > 1 && Ntim == 1
+            amps = γ*KomaMRIBase.get_theo_A(rf; off_val=NaN, max_rf_samples=Inf)[2:end]
+            return [amps[1]; amps[2:end-1][[i for i in 1:length(amps)-1 if i % 2 == 0]]; amps[end]]
+        end
+        return []#γ*KomaMRIBase.get_theo_A(rf; off_val=NaN, max_rf_samples=Inf)[2:end]
+    end
+    function get_theo_t_aux(gr::Grad)
+        Namp, Ntim = length(gr.A), length(gr.T)
+        if Namp == 1 && Ntim == 1
+            return KomaMRIBase.get_theo_t(gr)
+        elseif Namp > 1 && Ntim == 1
+            return KomaMRIBase.get_theo_t(gr)[2:end]
+        elseif Namp > 1 && Ntim > 1
+            return KomaMRIBase.get_theo_t(gr)[2:end]
+        end
+        return []
+    end
+    function get_theo_A_aux(gr::Grad)
+        Namp, Ntim = length(gr.A), length(gr.T)
+        if Namp == 1 && Ntim == 1
+            return 1e-3*γ*KomaMRIBase.get_theo_A(gr; off_val=0)
+        elseif  Namp > 1 && Ntim == 1
+            return 1e-3*γ*[gr.first; KomaMRIBase.get_theo_A(gr; off_val=0)[3:end-1]; gr.last]
+        elseif Namp > 1 && Ntim > 1
+            return 1e-3*γ*KomaMRIBase.get_theo_A(gr; off_val=0)[2:end]
+        end
+        return []
+    end
+    function get_theo_t_aux(adc::ADC)
+        return (adc.N == 1) ? [adc.T/2] .+ adc.delay : [range(0, adc.T; length=adc.N)...] .+ adc.delay
+    end
+
+    function koma_read(filename)
+        path = @__DIR__
+        seq = read_seq("$(path)/test_files/pulseq_read_comparison/$(filename).seq")
+        N = length(seq)
+        T0 = KomaMRIBase.get_block_start_times(seq)
+        return Dict(
+            "rf" => Dict(
+                    "blocks"     => [i for i in 1:N if KomaMRIBase.is_RF_on(seq[i])],
+                    "times"      => [get_theo_t_aux(seq[i].RF[1]) .+ T0[i] for i in 1:N if KomaMRIBase.is_RF_on(seq[i])],
+                    "amplitudes" => [get_theo_A_aux(seq[i].RF[1]) for i in 1:N if KomaMRIBase.is_RF_on(seq[i])]
+                ),
+            "gx" => Dict(
+                    "blocks"     => [i for i in 1:N if KomaMRIBase.is_Gx_on(seq[i])],
+                    "times"      => [get_theo_t_aux(seq.GR[1,i]) .+ T0[i] for i in 1:N if KomaMRIBase.is_Gx_on(seq[i])],
+                    "amplitudes" => [get_theo_A_aux(seq.GR[1,i]) for i in 1:N if KomaMRIBase.is_Gx_on(seq[i])]
+                ),
+            "gy" => Dict(
+                    "blocks"     => [i for i in 1:N if KomaMRIBase.is_Gy_on(seq[i])],
+                    "times"      => [get_theo_t_aux(seq.GR[2,i]) .+ T0[i] for i in 1:N if KomaMRIBase.is_Gy_on(seq[i])],
+                    "amplitudes" => [get_theo_A_aux(seq.GR[2,i]) for i in 1:N if KomaMRIBase.is_Gy_on(seq[i])]
+                ),
+            "gz" => Dict(
+                    "blocks"     => [i for i in 1:N if KomaMRIBase.is_Gz_on(seq[i])],
+                    "times"      => [get_theo_t_aux(seq.GR[3,i]) .+ T0[i] for i in 1:N if KomaMRIBase.is_Gz_on(seq[i])],
+                    "amplitudes" => [get_theo_A_aux(seq.GR[3,i]) for i in 1:N if KomaMRIBase.is_Gz_on(seq[i])]
+                ),
+            "adc" => Dict(
+                    "blocks"     => [i for i in 1:N if KomaMRIBase.is_ADC_on(seq[i])],
+                    "times"      => [get_theo_t_aux(seq.ADC[i]) .+ T0[i] for i in 1:N if KomaMRIBase.is_ADC_on(seq[i])],
+                )
+        )
+    end
+
+    function pulseq_read(filename)
+        path = @__DIR__
+        seq_pulseq = matread("$(path)/test_files/pulseq_read_comparison/$(filename).mat")
+        return Dict(
+            "rf" => Dict(
+                    "blocks"     => Int.(vec(seq_pulseq["rfBlocks"])),
+                    "times"      => vec([vec(blockRfTimes) for blockRfTimes in seq_pulseq["rfTimes"]]),
+                    "amplitudes" => vec([vec(blockRfAmplitudes) for blockRfAmplitudes in seq_pulseq["rfAmplitudes"]])
+                ),
+            "gx" => Dict(
+                    "blocks"     => Int.(vec(seq_pulseq["gxBlocks"])),
+                    "times"      => vec([vec(blockGxTimes) for blockGxTimes in seq_pulseq["gxTimes"]]),
+                    "amplitudes" => vec([vec(blockGxAmplitudes) for blockGxAmplitudes in seq_pulseq["gxAmplitudes"]])
+                ),
+            "gy" => Dict(
+                    "blocks"     => Int.(vec(seq_pulseq["gyBlocks"])),
+                    "times"      => vec([vec(blockGyTimes) for blockGyTimes in seq_pulseq["gyTimes"]]),
+                    "amplitudes" => vec([vec(blockGyAmplitudes) for blockGyAmplitudes in seq_pulseq["gyAmplitudes"]])
+                ),
+            "gz" => Dict(
+                    "blocks"     => Int.(vec(seq_pulseq["gzBlocks"])),
+                    "times"      => vec([vec(blockGzTimes) for blockGzTimes in seq_pulseq["gzTimes"]]),
+                    "amplitudes" => vec([vec(blockGzAmplitudes) for blockGzAmplitudes in seq_pulseq["gzAmplitudes"]])
+                ),
+            "adc" => Dict(
+                    "blocks"     => Int.(vec(seq_pulseq["adcBlocks"])),
+                    "times"      => vec([vec(blockAdcTimes) for blockAdcTimes in seq_pulseq["adcTimes"]]),
+                )
+        )
+    end
+
+    function read_comparison(filename)
+        # Define some names of the sequence dictionaries
+        event_names = ["rf", "gx", "gy", "gz", "adc"]
+        sample_names = ["times", "amplitudes"]
+        # Read data for Koma and Pulseq
+        seq_koma = koma_read(filename)      # it reads the .seq
+        seq_pulseq = pulseq_read(filename)  # it reads the .mat (ground truth)
+        # Checks for the blocks IDs of each event
+        for event in event_names
+            blocks_koma = seq_koma[event]["blocks"]
+            blocks_pulseq = seq_pulseq[event]["blocks"]
+            same_number_of_blocks = length(blocks_koma) == length(blocks_pulseq)
+            @test same_number_of_blocks
+            if same_number_of_blocks
+                same_blocks = all(blocks_koma .== blocks_pulseq)
+                @test same_blocks
+                if !same_blocks
+                    println("Mismatch for $(event) blocks")
+                    data = [blocks_koma blocks_pulseq]
+                    pretty_table(data; header=["koma", "pulseq"])
+                end
+            else
+                println("Length mismatch for $(event):")
+                println("Koma $(length(seq_koma[event])) != Pulseq $(length(seq_pulseq[event]))")
+            end
+        end
+        # Checks for the timings of each event
+        for event in event_names
+            for sample in sample_names
+                if !(event == "adc" && sample == "amplitudes")
+                    for i in 1:length(seq_pulseq[event]["blocks"])
+                        samples_koma = seq_koma[event][sample][i]
+                        samples_pulseq = seq_pulseq[event][sample][i]
+                        same_points = all(.≈(samples_koma, samples_pulseq, atol=TOLERANCE))
+                        @test same_points
+                        if !same_points
+                            println("Mismatch for $(event) $(sample). Block $(seq_pulseq[event]["blocks"][i])")
+                            data = [samples_koma samples_pulseq]
+                            pretty_table(data; header=["koma", "pulseq"])
+                            # This is just an example about how to debug locally when there are problems
+                            #if event == "gx" && sample == "amplitudes" && seq_pulseq[event]["blocks"][i] == 3
+                            #    pretty_table(data[end-10:end,:]; header=["koma", "pulseq"])
+                            #end
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    @testset "Trapezoidal_Grad" begin
+        read_comparison("gr-trapezoidal")
+    end
+    @testset "Uniformly_Shaped_Grad" begin
+        read_comparison("gr-uniformly-shaped")
+    end
+    @testset "Time_Shaped_Grad" begin
+        read_comparison("gr-time-shaped")
+    end
+    @testset "Pulse_RF" begin
+        read_comparison("rf-pulse")
+    end
+    @testset "Uniformly_Shaped_RF" begin
+        read_comparison("rf-uniformly-shaped")
+    end
+    @testset "FID" begin
+        read_comparison("fid")
+    end
+    @testset "Spiral" begin
+        read_comparison("spiral")
+    end
+
+end