diff --git a/src/sequence_io/pulseq.jl b/src/sequence_io/pulseq.jl index 5f84e23d561ebb6776ce5496f36c395a5cf5633c..2aeb46cf9c907b344e2f732ffdcb30f627b72302 100644 --- a/src/sequence_io/pulseq.jl +++ b/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., )) diff --git a/src/sequence_io/pulseq_io/parsers/shapes.jl b/src/sequence_io/pulseq_io/parsers/shapes.jl index ef0a30f79365855f2f143906a85909e1016a3fda..3a062f7f6251e7b979d3fc84c68ad35afcbcc731 100644 --- a/src/sequence_io/pulseq_io/parsers/shapes.jl +++ b/src/sequence_io/pulseq_io/parsers/shapes.jl @@ -1,6 +1,6 @@ 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)) diff --git a/test/test_IO.jl b/test/test_IO.jl index 12c10efb74ca70c5869fd50059097fc728ec3628..32bcfa8fc1c7b5eb55b6b2ce43a3e3adf0a84ab1 100644 --- a/test/test_IO.jl +++ b/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