Add tests to MRIBuilder for pulseq and plotting

Project.toml

@@ -27,6 +27,9 @@ julia = "1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44"
+VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 
 [targets]
-test = ["Test"]
+test = ["Test", "CairoMakie", "Gtk", "VisualRegressionTests"]

test/example_pulseseq/01_from_FID_to_PRESS_v140/01_FID.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 0
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+RadiofrequencyRasterTime 1e-06 
+TotalDuration 5 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+1  62   1   0   0   0  0  0
+2 499938   0   0   0   0  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id time_shape_id delay freq phase
+# ..        Hz   ....     ....          ....    us   Hz   rad
+[RF]
+1          500 1 2 3 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 8192 31250 29730 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
+500
+
+
+[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 5caad2dc77321c47fe1cf56d9ef581b9

test/example_pulseseq/01_from_FID_to_PRESS_v140/02_FID_MFA.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 0
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+RadiofrequencyRasterTime 1e-06 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+ 1 112   1   0   0   0  0  0
+ 2 499888   0   0   0   0  1  0
+ 3 112   2   0   0   0  0  0
+ 4 499888   0   0   0   0  1  0
+ 5 112   3   0   0   0  0  0
+ 6 499888   0   0   0   0  1  0
+ 7 112   4   0   0   0  0  0
+ 8 499888   0   0   0   0  1  0
+ 9 112   5   0   0   0  0  0
+10 499888   0   0   0   0  1  0
+11 112   6   0   0   0  0  0
+12 499888   0   0   0   0  1  0
+13 112   7   0   0   0  0  0
+14 499888   0   0   0   0  1  0
+15 112   8   0   0   0  0  0
+16 499888   0   0   0   0  1  0
+17 112   9   0   0   0  0  0
+18 499888   0   0   0   0  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id time_shape_id delay freq phase
+# ..        Hz   ....     ....          ....    us   Hz   rad
+[RF]
+1      55.5556 1 2 3 100 0 0
+2      111.111 1 2 3 100 0 0
+3      166.667 1 2 3 100 0 0
+4      222.222 1 2 3 100 0 0
+5      277.778 1 2 3 100 0 0
+6      333.333 1 2 3 100 0 0
+7      388.889 1 2 3 100 0 0
+8      444.444 1 2 3 100 0 0
+9          500 1 2 3 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 8192 31250 29480 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
+1000
+
+
+[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 7e2b89c1ae6eca64fc825f6ccdc035ef

test/example_pulseseq/01_from_FID_to_PRESS_v140/03_SE_noSpoil.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 0
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+RadiofrequencyRasterTime 1e-06 
+TotalDuration 0.25 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+1 112   1   0   0   0  0  0
+2 5388   0   0   0   0  0  0
+3 112   2   0   0   0  0  0
+4 19388   0   0   0   0  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id time_shape_id delay freq phase
+# ..        Hz   ....     ....          ....    us   Hz   rad
+[RF]
+1          250 1 2 3 100 0 0
+2          500 1 2 3 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 8192 12500 3280 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
+1000
+
+
+[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 67f9e34da9b97a9f23af82d17863462c

test/example_pulseseq/01_from_FID_to_PRESS_v140/04_SE_BackgroundGrad.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 0
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+RadiofrequencyRasterTime 1e-06 
+TotalDuration 0.25006 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+1   3   0   1   0   0  0  0
+2 112   1   2   0   0  0  0
+3 5388   0   3   0   0  0  0
+4 112   2   2   0   0  0  0
+5 19388   0   4   0   0  1  0
+6   3   0   5   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          250 3 4 5 100 0 0
+2          500 3 4 5 100 0 0
+
+# Format of arbitrary gradients:
+#   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
+# id amplitude amp_shape_id time_shape_id delay
+# ..      Hz/m       ..         ..          us
+[GRADIENTS]
+1         5000 1 2 0
+2         5000 3 6 0
+3         5000 3 7 0
+4         5000 3 8 0
+5         5000 9 2 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 8192 12500 3280 0 0
+
+# Sequence Shapes
+[SHAPES]
+
+shape_id 1
+num_samples 2
+0
+1
+
+shape_id 2
+num_samples 2
+0
+3
+
+shape_id 3
+num_samples 2
+1
+1
+
+shape_id 4
+num_samples 2
+0
+0
+
+shape_id 5
+num_samples 2
+0
+1000
+
+shape_id 6
+num_samples 2
+0
+112
+
+shape_id 7
+num_samples 2
+0
+5388
+
+shape_id 8
+num_samples 2
+0
+19388
+
+shape_id 9
+num_samples 2
+1
+0
+
+
+[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 aac1d624e8bf86af3edd1f06085c1ae2

test/example_pulseseq/01_from_FID_to_PRESS_v140/05_SE_withSpoilers.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 4
+revision 0
+
+[DEFINITIONS]
+AdcRasterTime 1e-07 
+BlockDurationRaster 1e-05 
+GradientRasterTime 1e-05 
+RadiofrequencyRasterTime 1e-06 
+TotalDuration 0.25 
+
+# Format of blocks:
+# NUM DUR RF  GX  GY  GZ  ADC  EXT
+[BLOCKS]
+1 112   1   0   0   0  0  0
+2 5315   0   0   0   0  0  0
+3 185   2   0   0   1  0  0
+4 19388   0   0   0   1  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id time_shape_id delay freq phase
+# ..        Hz   ....     ....          ....    us   Hz   rad
+[RF]
+1          250 1 2 3 100 0 0
+2          500 1 2 3 830 0 0
+
+# Format of trapezoid gradients:
+# id amplitude rise flat fall delay
+# ..      Hz/m   us   us   us    us
+[TRAP]
+ 1  1.69492e+06 240  350 240   0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 8192 12500 3280 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
+1000
+
+
+[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 7a621205cecc36a2aadeeccf9d35399f

test/example_pulseseq/01_from_FID_to_PRESS_v140/README.md

+These example sequence files have been taken from the pulseq/tutorials github repository (https://github.com/pulseq/tutorials/tree/main/01_from_FID_to_PRESS)
+
+These are all v1.4.0 and are used to test the reading of such pulseq files.
\ No newline at end of file

test/example_pulseseq/fiddisp_v131.seq

+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 3
+revision 1
+
+[DEFINITIONS]
+Name fid 
+
+# Format of blocks:
+#  #  D RF  GX  GY  GZ ADC EXT
+[BLOCKS]
+1  0  1   0   0   0  0  0
+2  1  0   0   0   0  0  0
+3  0  0   0   0   0  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id delay freq phase
+# ..        Hz   ....     ....    us   Hz   rad
+[RF]
+1         2500 1 2 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 1024 312500 20 0 0
+
+# Format of delays:
+# id delay (us)
+[DELAYS]
+1 5000
+
+# Sequence Shapes
+[SHAPES]
+
+shape_id 1
+num_samples 120
+1
+0
+0
+97
+-1
+0
+0
+17
+
+shape_id 2
+num_samples 120
+0
+0
+118
+

test/runtests.jl

@@ -2,5 +2,6 @@ using MRIBuilder
 using Test
 
 @testset "MRIBuilder.jl" begin
-    # Write your tests here.
+    include("test_IO.jl")
+    include("test_plot.jl")
 end

test/test_IO.jl

+@testset "test_IO.jl" begin
+    directory = joinpath(pwd(), "example_pulseseq")
+    @testset "read v1.3.1 fiddisp.seq file" begin
+        seq = read_sequence(joinpath(directory, "fiddisp_v131.seq"))
+
+        @test length(seq) == 3
+
+        # starting with RF pulse
+        @test length(events(seq[1])) == 1
+        (index, pulse) = events(seq[1])[1]
+
+        @test start_time(seq, 1, index) == 0.1
+        @test end_time(seq, 1, index) == 0.22
+
+        for (time, ampl, phase_check) in [
+            (0.09, 0., isnan),
+            (0.11, 2.5, iszero),
+            (0.19, 2.5, iszero),
+            (0.21, 0., iszero),
+            (0.23, 0., isnan),
+        ]
+            @test phase_check(phase(seq, time))
+            @test amplitude(seq, time) ≈ ampl
+        end
+
+        # Single ADC event
+        @test length(readout_times(seq)) == 1024
+        @test readout_times(seq)[1] ≈ 0.22 + 5 + 0.02 + 0.5 * 0.3125
+        @test readout_times(seq)[end] ≈ 0.22 + 5 + 0.02 + 1023.5 * 0.3125
+
+        @test TR(seq) ≈ 0.22 + 5 + 0.02 + 1024 * 0.3125
+    end
+    @testset "read all v1.4.0 files in 01_from_FID_to_PRESS" begin
+        path = joinpath(directory, "01_from_FID_to_PRESS_v140")
+        for fn in readdir(path)
+            if !endswith(fn, ".seq")
+                continue
+            end
+            full_filename = joinpath(path, fn)
+            seq = read_sequence(full_filename)
+        end
+    end
+    @testset "check 01_from_FID_to_PRESS/01_FID.seq" begin
+        fn = joinpath(directory, "01_from_FID_to_PRESS_v140", "01_FID.seq")
+        seq = read_sequence(fn)
+        @test length(seq) == 2
+        (index, pulse) = events(seq[1])[1]
+        @test flip_angle(pulse) ≈ 90
+        @test start_time(seq, 1, index) ≈ 0.1  # determined by RF dead time
+        @test end_time(seq, 1, index) ≈ 0.6  # RF dead time + RF duration
+        @test phase(seq, 0.3) == 0
+        @test amplitude(seq, 0.3) ≈ 0.5
+
+        start_adc = (
+            0.6 +  # end of RF pulse
+            0.02 +  # RF ringdown time (added by block duration)
+            29.730  # Delay until ADC start (to get start at ADC at TE=30)
+        )
+        @test length(readout_times(seq)) == 8192
+        @test readout_times(seq)[1] ≈ start_adc + 0.5 * 0.03125
+        @test readout_times(seq)[end] ≈ start_adc + 8191.5 * 0.03125
+    end
+    @testset "check 01_from_FID_to_PRESS/06_PRESS_center.seq" begin
+        fn = joinpath(directory, "01_from_FID_to_PRESS_v140", "06_PRESS_center.seq")
+        seq = read_sequence(fn)
+        @test length(seq) == 7
+
+        (index, excitation) = events(seq[1])[1]
+        #@test flip_angle(excitation) ≈ 90
+        @test start_time(seq, 1, index) ≈ 0.1  # determined by RF dead time
+        @test end_time(seq, 1, index) ≈ 3.1  # RF dead time + RF duration
+        @test phase(seq, 0.3) ≈ 0 + rad2deg(0.5)
+        @test phase(seq, 1.6) ≈ 0
+        @test length(readout_times(seq)) == 4096
+
+        refocus_pulses = (
+            events(seq[3])[1][2],
+            events(seq[5])[1][2],
+        )
+        for p in refocus_pulses
+            @test duration(p) ≈ 3  # should have been 4 for refocus pulses, but there is an error in the matlab generation
+            @test phase(p, 0.) ≈ 90 rtol=1e-5
+            @test phase(p, 0.38) ≈ 90 + rad2deg(0.5) rtol=1e-5
+            @test phase(p, 1.5) ≈ 90 rtol=1e-5
+        end
+    end
+    if false
+        # JSON encoding has not been implemented yet
+        @testset "check that JSON encoding works for all v1.4.0 files in 01_from_FID_to_PRESS" begin
+            path = joinpath(directory, "01_from_FID_to_PRESS_v140")
+            for fn in readdir(path)
+                @testset "checking $fn" begin
+                    if !endswith(fn, ".seq")
+                        continue
+                    end
+                    full_filename = joinpath(path, fn)
+                    seq_orig = read_sequence(full_filename)
+
+                    
+                    io = IOBuffer()
+                    write_sequence(io, seq_orig)
+                    s = String(io.data)
+                    seq_json = read_sequence(s)
+                    @test seq_orig.TR == seq_json.TR
+                    @test length(seq_orig.gradients) == length(seq_json.gradients)
+                    for (g1, g2) in zip(seq_orig.gradients, seq_json.gradients)
+                        @test all(g1.origin .== g2.origin)
+                        @test all(g1.shape.times .== g2.shape.times)
+                        @test all(g1.shape.amplitudes .== g2.shape.amplitudes)
+                    end
+                    @test length(seq_orig.pulses) == length(seq_json.pulses)
+                    for (p1, p2) in zip(seq_orig.pulses, seq_json.pulses)
+                        @test p1.max_amplitude == p2.max_amplitude
+                        @test all(p1.amplitude.times .== p2.amplitude.times)
+                        @test all(p1.amplitude.amplitudes .== p2.amplitude.amplitudes)
+                        @test all(p1.phase.times .== p2.phase.times)
+                        @test all(p1.phase.amplitudes .== p2.phase.amplitudes)
+                    end
+                    @test iszero(length(seq_json.instants))
+                    @test length(seq_json.readout_times) > 0
+                    @test all(seq_json.readout_times .== seq_orig.readout_times)
+                end
+            end
+        end
+        @testset "check that JSON encoding works for some sequences with instant pulses/gradients" begin
+            for seq_orig in [
+                DWI(TE=80., bval=2., TR=100., scanner=Siemens_Connectom),
+                DWI(TE=80., bval=2., gradient_duration=0., TR=100, scanner=Siemens_Terra),
+                SpinEcho(TE=30., TR=100., scanner=Scanner(B0=1.5)),
+            ]
+                io = IOBuffer()
+                write_sequence(io, seq_orig)
+                s = String(io.data)
+                seq_json = read_sequence(s)
+                @test seq_json.TR == 100.
+                @test seq_json.scanner.B0 == seq_orig.scanner.B0
+                @test seq_json.scanner.gradient == seq_orig.scanner.gradient
+                @test seq_json.scanner.slew_rate == seq_orig.scanner.slew_rate
+                @test length(seq_orig.gradients) == length(seq_json.gradients)
+                for (g1, g2) in zip(seq_orig.gradients, seq_json.gradients)
+                    @test all(g1.origin .== g2.origin)
+                    @test all(g1.shape.times .== g2.shape.times)
+                    @test all(g1.shape.amplitudes .== g2.shape.amplitudes)
+                end
+                @test iszero(length(seq_json.pulses))
+                @test length(seq_json.instants) == length(seq_orig.instants)
+                for (i1, i2) in zip(seq_orig.instants, seq_json.instants)
+                    @test i1.time == i2.time
+                    if i1 isa mr.InstantRFPulse
+                        @test i2 isa mr.InstantRFPulse
+                        @test i1.flip_angle == i2.flip_angle
+                        @test i1.phase == i2.phase
+                    else
+                        @test i1 isa mr.InstantGradient
+                        @test i2 isa mr.InstantGradient
+                        @test all(i1.origin .== i2.origin)
+                        @test all(i1.qvec .== i2.qvec)
+                    end
+                end
+                @test length(seq_json.readout_times) > 0
+                @test all(seq_json.readout_times .== seq_orig.readout_times)
+            end
+        end
+    end
+end
\ No newline at end of file

test/test_plot.jl

+using CairoMakie
+using VisualRegressionTests
+using Gtk
+@testset "test_plot.jl" begin
+    dir = @__DIR__
+    isCI = get(ENV, "CI", "false") == "true"
+    @testset "pulseseq example PRESS sequence" begin
+        function plot_press(fname)
+            fn = joinpath(dir, "example_pulseseq", "01_from_FID_to_PRESS_v140", "06_PRESS_center.seq")
+            sequence = read_sequence(fn)
+            f = Figure()
+            Axis(f[1, 1])
+            plot!(sequence)
+            CairoMakie.save(fname, f)
+        end
+
+        @visualtest plot_press "$dir/plots/pulseq_press.png" !isCI
+    end
+
+end
\ No newline at end of file