Fix ADC Pulseq raster timing (#756)

Author: cncastillo

KomaMRIBase/Project.toml

@@ -1,7 +1,7 @@
 name = "KomaMRIBase"
 uuid = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
 
-version = "0.10.0"
+version = "0.10.1"
 authors = ["Carlos Castillo Passi <cncastillo@uc.cl>"]
 
 [workspace]

KomaMRIBase/src/datatypes/sequence/TimingChecks.jl

@@ -131,11 +131,11 @@ function _check_adc_dwell(dwell, raster, block_id)
     error("Block $block_id ADC dwell ($(dwell) s) is not aligned to raster $(raster) s.")
 end
 
-function _check_adc_timing(adc, raster, block_id)
+function _check_adc_timing(adc, adc_raster, rf_raster, block_id)
     dwell = adc.N == 1 ? adc.T : adc.T / (adc.N - 1)
-    _check_adc_dwell(dwell, raster, block_id)
+    _check_adc_dwell(dwell, adc_raster, block_id)
     adc.delay + PULSEQ_TIME_TOL >= dwell / 2 || error("Block $block_id ADC delay ($(adc.delay) s) is smaller than dwell/2 ($(dwell / 2) s).")
-    _check_raster_multiple(adc.delay - dwell / 2, raster, block_id, "ADC delay")
+    _check_raster_multiple(adc.delay - dwell / 2, rf_raster, block_id, "ADC delay")
     return nothing
 end
 
@@ -169,7 +169,7 @@ function check_timing(seq::Sequence, raster::NamedTuple)
             _check_grad_timing(gr, raster.GradientRasterTime, i, name)
         end
         adc = seq.ADC[i]
-        is_ADC_on(adc) && _check_adc_timing(adc, raster.AdcRasterTime, i)
+        is_ADC_on(adc) && _check_adc_timing(adc, raster.AdcRasterTime, raster.RadiofrequencyRasterTime, i)
     end
     return nothing
 end

KomaMRIBase/test/runtests.jl

@@ -941,8 +941,14 @@ using TestItems, TestItemRunner
         @addblock compact_adc += (ADC(2, 100e-9, 50e-9), Duration(10e-6))
         @test isnothing(check_timing(compact_adc))
         adc_raster_seq = Sequence(Scanner(RF_Δt=10e-6, ADC_Δt=100e-9))
-        @addblock adc_raster_seq += (ADC(111, 4.4e-6 * 110, 151e-6), Duration(650e-6))
+        @addblock adc_raster_seq += (ADC(111, 4.4e-6 * 110, 152.2e-6), Duration(650e-6))
         @test isnothing(check_timing(adc_raster_seq))
+        adc_bad_dwell = Sequence(Scanner(ADC_Δt=100e-9))
+        @addblock adc_bad_dwell += (ADC(2, 150e-9, 75e-9), Duration(10e-6))
+        @test_throws ErrorException check_timing(adc_bad_dwell)
+        adc_bad_delay = Sequence(Scanner(RF_Δt=10e-6, ADC_Δt=100e-9))
+        @addblock adc_bad_delay += (ADC(2, 100e-9, 150e-9), Duration(10e-6))
+        @test_throws ErrorException check_timing(adc_bad_delay)
         seq_bad = Sequence()
         @addblock seq_bad += RF(1e-6, 10.5e-6)
         @test_throws ErrorException check_timing(seq_bad)

KomaMRIFiles/Project.toml

@@ -1,7 +1,7 @@
 name = "KomaMRIFiles"
 uuid = "fcf631a6-1c7e-4e88-9e64-b8888386d9dc"
 
-version = "0.10.1"
+version = "0.10.2"
 authors = ["Carlos Castillo Passi <cncastillo@uc.cl>"]
 
 [workspace]

KomaMRIFiles/src/Sequence/pulseq/WritePulseq.jl

@@ -991,7 +991,7 @@ function quantize_adc!(adc, raster, block_id, warn_count)
     qt(value, field; step, step_name) =
         quantize_time_value(value, step, warn_ctx, field; step_name)
     dwell = adc.N == 1 ? adc.T : adc.T / (adc.N - 1)
-    dwell = qt(dwell, :dwell_time; step=1e-9, step_name=:nanoseconds)
+    dwell = qt(dwell, :dwell_time; step=raster.AdcRasterTime, step_name=:AdcRasterTime)
     adc.T = adc.N == 1 ? dwell : (adc.N - 1) * dwell
     pulseq_delay = adc.delay - dwell / 2
     if pulseq_delay < 0
@@ -1001,7 +1001,7 @@ This means the Koma delay must be >= dwell/2. Clamping it to this minimum value.
 See https://pulseq.github.io/pulseq_shapes_and_times.pdf#page=10"
         pulseq_delay = 0
     end
-    pulseq_delay = qt(pulseq_delay, :pulseq_delay; step=1e-6, step_name=:microseconds)
+    pulseq_delay = qt(pulseq_delay, :pulseq_delay; step=raster.RadiofrequencyRasterTime, step_name=:RadiofrequencyRasterTime)
     adc.delay = pulseq_delay + dwell / 2
     return adc.delay + adc.T + dwell / 2
 end

KomaMRIFiles/test/runtests.jl

@@ -261,20 +261,17 @@ end
         generated_prefix = "koma-generated-"
         pth = joinpath(@__DIR__, "test_files/pulseq/read_comparison/")
         versions = ["v1.2", "v1.3", "v1.4", "v1.5"]
-        constraint_sys = Scanner(B1=1.0, Gmax=100.0, Smax=1e12, ADC_Δt=100e-9)
         mktempdir() do tmpdir
             for v in versions
                 pulseq_files = filter(endswith(".seq"), readdir(pth * v)) .|> x -> splitext(x)[1]
                 for pulseq_file in pulseq_files
                     label = "$v/$pulseq_file"
-                    seq = @suppress read_seq("$pth$v/$pulseq_file.seq")
+                    seq = read_seq("$pth$v/$pulseq_file.seq"; verbose=false)
                     filename = joinpath(tmpdir, "$(generated_prefix)$(v)-$(pulseq_file).seq")
-                    if occursin("JEMRIS", label) || occursin("gammaSTAR", label)
-                        @suppress write_seq(seq, filename; check_timing=false)
-                    else
-                        @suppress write_seq(seq, filename; sys=constraint_sys)
-                    end
-                    seq2 = @suppress read_seq(filename)
+                    # Legacy JEMRIS files omit AdcRasterTime but use 1 ns dwell precision.
+                    occursin("JEMRIS", label) && (seq.DEF["AdcRasterTime"] = 1e-9)
+                    write_seq(seq, filename; verbose=false)
+                    seq2 = read_seq(filename; verbose=false)
                     @testset "$label" begin
                         @test seq2 ≈ seq
                     end