Incorrect gradient interpretation when waveforms do not end in zero

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

@@ -498,6 +498,7 @@ function read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
         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
@@ -522,6 +523,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"]
@@ -534,8 +540,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))
@@ -551,6 +557,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
@@ -569,8 +576,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*γ*KomaMRIBase.get_theo_A(gr; off_val=0)[2:end]
+        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

KomaMRIFiles/test/test_files/pulseq_read_comparison/fid.seq

@@ -0,0 +1,120 @@
+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 1
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+Name fid 
+RadiofrequencyRasterTime 1e-06 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+ 1  23   1   0   0   0  0  0
+ 2 2000   0   0   0   0  0  0
+ 3 324   0   0   0   0  1  0
+ 4 100000   0   0   0   0  0  0
+ 5  23   1   0   0   0  0  0
+ 6 2000   0   0   0   0  0  0
+ 7 324   0   0   0   0  1  0
+ 8 100000   0   0   0   0  0  0
+ 9  23   1   0   0   0  0  0
+10 2000   0   0   0   0  0  0
+11 324   0   0   0   0  1  0
+12 100000   0   0   0   0  0  0
+13  23   1   0   0   0  0  0
+14 2000   0   0   0   0  0  0
+15 324   0   0   0   0  1  0
+16 100000   0   0   0   0  0  0
+17  23   1   0   0   0  0  0
+18 2000   0   0   0   0  0  0
+19 324   0   0   0   0  1  0
+20 100000   0   0   0   0  0  0
+21  23   1   0   0   0  0  0
+22 2000   0   0   0   0  0  0
+23 324   0   0   0   0  1  0
+24 100000   0   0   0   0  0  0
+25  23   1   0   0   0  0  0
+26 2000   0   0   0   0  0  0
+27 324   0   0   0   0  1  0
+28 100000   0   0   0   0  0  0
+29  23   1   0   0   0  0  0
+30 2000   0   0   0   0  0  0
+31 324   0   0   0   0  1  0
+32 100000   0   0   0   0  0  0
+33  23   1   0   0   0  0  0
+34 2000   0   0   0   0  0  0
+35 324   0   0   0   0  1  0
+36 100000   0   0   0   0  0  0
+37  23   1   0   0   0  0  0
+38 2000   0   0   0   0  0  0
+39 324   0   0   0   0  1  0
+40 100000   0   0   0   0  0  0
+41  23   1   0   0   0  0  0
+42 2000   0   0   0   0  0  0
+43 324   0   0   0   0  1  0
+44 100000   0   0   0   0  0  0
+45  23   1   0   0   0  0  0
+46 2000   0   0   0   0  0  0
+47 324   0   0   0   0  1  0
+48 100000   0   0   0   0  0  0
+49  23   1   0   0   0  0  0
+50 2000   0   0   0   0  0  0
+51 324   0   0   0   0  1  0
+52 100000   0   0   0   0  0  0
+53  23   1   0   0   0  0  0
+54 2000   0   0   0   0  0  0
+55 324   0   0   0   0  1  0
+56 100000   0   0   0   0  0  0
+57  23   1   0   0   0  0  0
+58 2000   0   0   0   0  0  0
+59 324   0   0   0   0  1  0
+60 100000   0   0   0   0  0  0
+61  23   1   0   0   0  0  0
+62 2000   0   0   0   0  0  0
+63 324   0   0   0   0  1  0
+64 100000   0   0   0   0  0  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id time_shape_id delay freq phase
+# ..        Hz   ....     ....          ....    us   Hz   rad
+[RF]
+1         2500 1 2 3 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 256 12500 20 0 0
+
+# Sequence Shapes
+[SHAPES]
+
+shape_id 1
+num_samples 2
+1
+1
+
+shape_id 2
+num_samples 2
+0
+0
+
+shape_id 3
+num_samples 2
+0
+100
+
+
+[SIGNATURE]
+# This is the hash of the Pulseq file, calculated right before the [SIGNATURE] section was added
+# It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
+# The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be sripped away for recalculating/verification)
+Type md5
+Hash cb837669c48fb79459d2b7d56453c696