src/sequence_io/pulseq.jl

@@ -66,7 +66,7 @@ function BuildingBlock(pulseq::PulseqBlock; version, BlockDurationRaster, Radiof
     )
 
     if min_duration > stated_duration
-        if version == v"1.3.1"
+        if version == v"1.3.1" || isapprox(min_duration, stated_duration, rtol=1e-3)
             stated_duration = min_duration
         else
             error("Minimum duration to play all RF/gradient/ADC events exceeds stated duration.")
@@ -177,12 +177,13 @@ function PulseqBlock(block::BaseBuildingBlock; BlockDurationRaster, AdcRasterTim
     times = [t for (t, _) in waveform(block)]
     for dim in 1:3
         amplitudes = [g[dim] for (_, g) in waveform(block)]
-        if iszero(maximum(abs.(amplitudes); init=0.))
+        max_ampl = maximum(abs.(amplitudes); init=0.)
+        if iszero(max_ampl)
             push!(grads, nothing)
         else
             push!(grads, PulseqGradient(
-                maximum(amplitudes) * 1e3,
-                PulseqShape(amplitudes ./ maximum(amplitudes)),
+                max_ampl * 1e3,
+                PulseqShape(amplitudes ./ max_ampl),
                 PulseqShape(times .* 1e-3),
                 0.,
             ))

src/sequence_io/pulseq_io/parsers/shapes.jl

 struct CompressedPulseqShape
     num :: Int
-    samples :: Vector{Float64}
+    samples :: Vector{Number}
 end
 
 
@@ -57,14 +57,48 @@ function uncompress(compressed::CompressedPulseqShape)
     return PulseqShape(amplitudes)
 end
 
+function compress(shape::PulseqShape)
+    amplitudes = shape.samples
+    derivative = amplitudes[2:end] - amplitudes[1:end-1]
+
+    current_value = amplitudes[1]
+    compressed = Number[amplitudes[1]]
+
+    nrepeats = -1
+    for sample in derivative
+        if sample ≈ current_value
+            nrepeats += 1
+            if iszero(nrepeats)
+                push!(compressed, current_value)
+            end
+        else
+            if nrepeats != -1
+                push!(compressed, nrepeats)
+                nrepeats = -1
+            end
+            current_value = sample
+            push!(compressed, current_value)
+        end
+    end
+    if nrepeats >= 0
+        push!(compressed, nrepeats)
+    end
+    if length(compressed) >= length(amplitudes)
+        return CompressedPulseqShape(length(amplitudes), amplitudes)
+    else
+        return CompressedPulseqShape(length(amplitudes), compressed)
+    end
+end
+
 function gen_section(comp:: PulseqComponents, ::Val{:shapes})
     res = PulseqSection{:shapes}(String[])
     for index in sort([keys(comp.shapes)...])
-        shape = comp.shapes[index]
+        shape = compress(comp.shapes[index])
+        @assert all(isfinite.(shape.samples))
         append!(res.content, [
             "",
             "shape_id $index",
-            "num_samples $(length(shape.samples))"
+            "num_samples $(shape.num)"
         ])
         for sample in shape.samples
             push!(res.content, string(sample))

test/test_IO.jl

@@ -111,6 +111,7 @@
         @testset "check that format $(format) works for some sequences with instant pulses/gradients/readouts" begin
             for seq_orig in [
                 DWI(TE=80., bval=2., scanner=Siemens_Connectom),
+                DWI(TE=80., bval=2., scanner=Siemens_Connectom, fov=(10, 10, 10), voxel_size=2.),
                 DWI(TE=80., bval=2., scanner=Siemens_Terra, gradient=(type=:instant, )),
                 SpinEcho(TE=30.),
             ]
@@ -118,14 +119,30 @@
                 write_sequence(io, seq_orig; format=format)
                 seek(io, 0)
                 seq_json = read_sequence(io; format=format)
-                @test isapprox(variables.duration(seq_orig), variables.duration(seq_json), atol=1e-5)
-                @test length(seq_orig) == length(seq_json)
-                @test all(isapprox.(variables.duration.(seq_orig), variables.duration.(seq_json), atol=1e-5))
+                @test isapprox(variables.duration(seq_orig), variables.duration(seq_json), atol=1e-3)
+                @test length(iter_blocks(seq_orig)) == length(iter_blocks(seq_json))
+                @test all(isapprox.(variables.duration.(last.(iter_blocks(seq_orig))), variables.duration.(last.(iter_blocks(seq_json))), atol=1e-3))
                 @test length(iter_instant_gradients(seq_json)) == length(iter_instant_gradients(seq_json))
                 @test length(iter_instant_pulses(seq_json)) == length(iter_instant_pulses(seq_json))
                 @test length(variables.readout_times(seq_orig)) == length(variables.readout_times(seq_json))
-                @test all(isapprox.(variables.readout_times(seq_orig), variables.readout_times(seq_json), atol=1e-5))
+                @test all(isapprox.(variables.readout_times(seq_orig), variables.readout_times(seq_json), atol=1e-3))
             end
         end
     end
+    @testset "Test Pulseq shape compression" begin
+        import MRIBuilder.SequenceIO.PulseqIO.Parsers: compress, uncompress
+        import MRIBuilder.SequenceIO.PulseqIO.Types: PulseqShape
+
+        for (amplitudes, result) in [
+            ([0., 0.1, 0.25, 0.5, ones(7)..., 0.75, 0.5, 0.25, 0.], [0., 0.1, 0.15, 0.25, 0.5, 0., 0., 4, -0.25, -0.25, 2]),
+            (zeros(100), [0., 0., 98]),
+            (ones(100), [1., 0., 0., 97]),
+        ]
+            shape = PulseqShape(amplitudes)
+            res = compress(shape)
+            @test res.num == length(amplitudes)
+            @test all(res.samples == result)
+            @test all(uncompress(res).samples == amplitudes)
+        end
+    end
 end
\ No newline at end of file