From 0b8844f340886868d1df89cfa7ae0f0519817047 Mon Sep 17 00:00:00 2001 From: gabuzi <15203081+gabuzi@users.noreply.github.com> Date: Sun, 3 Dec 2023 20:28:03 +0000 Subject: [PATCH] Load pulseq sequences faster and allow comparison --- KomaMRICore/src/KomaMRICore.jl | 2 +- KomaMRICore/src/datatypes/Sequence.jl | 4 + KomaMRICore/src/datatypes/sequence/ADC.jl | 1 + KomaMRICore/src/datatypes/sequence/Delay.jl | 2 + KomaMRICore/src/datatypes/sequence/Grad.jl | 1 + KomaMRICore/src/datatypes/sequence/RF.jl | 1 + KomaMRICore/src/io/Pulseq.jl | 311 +++++++++++++++++++- 7 files changed, 320 insertions(+), 2 deletions(-) diff --git a/KomaMRICore/src/KomaMRICore.jl b/KomaMRICore/src/KomaMRICore.jl index efac073da..be7ffff10 100644 --- a/KomaMRICore/src/KomaMRICore.jl +++ b/KomaMRICore/src/KomaMRICore.jl @@ -51,7 +51,7 @@ export Scanner, Sequence, Phantom export Grad, RF, ADC, Delay export Mag, dur #Pulseq -export read_seq +export read_seq, read_seq_via_blocks_as_int_array #ISMRMRD export signal_to_raw_data #Phantom diff --git a/KomaMRICore/src/datatypes/Sequence.jl b/KomaMRICore/src/datatypes/Sequence.jl index 916750277..5348c1f1d 100644 --- a/KomaMRICore/src/datatypes/Sequence.jl +++ b/KomaMRICore/src/datatypes/Sequence.jl @@ -138,6 +138,10 @@ recursive_merge(x...) = x[end] #Sequence object functions size(x::Sequence) = size(x.GR[1,:]) +# concatenation (+ for collections at once) +Base.hcat(v::Vector{Sequence}) = Sequence(reduce(hcat, s.GR for s in v), reduce(hcat, s.RF for s in v), reduce(vcat, s.ADC for s in v), reduce(vcat, s.DUR for s in v), reduce(recursive_merge, s.DEF for s in v)) +Base.:(==)(x::Sequence, y::Sequence) = x.GR == y.GR && x.RF == y.RF && x.ADC == y.ADC && x.DUR == y.DUR && x.DEF == y.DEF + """ y = is_ADC_on(x::Sequence) y = is_ADC_on(x::Sequence, t::Union{Array{Float64,1}, Array{Float64,2}}) diff --git a/KomaMRICore/src/datatypes/sequence/ADC.jl b/KomaMRICore/src/datatypes/sequence/ADC.jl index 6f533443d..b1a70326d 100644 --- a/KomaMRICore/src/datatypes/sequence/ADC.jl +++ b/KomaMRICore/src/datatypes/sequence/ADC.jl @@ -44,6 +44,7 @@ Base.isapprox(adc1::ADC, adc2::ADC) = begin return all(length(getfield(adc1, k)) ≈ length(getfield(adc2, k)) for k ∈ fieldnames(ADC)) all(getfield(adc1, k) ≈ getfield(adc2, k) for k ∈ fieldnames(ADC)) end +Base.:(==)(x::ADC, y::ADC) = all(getfield(x, k) == getfield(y, k) for k ∈ fieldnames(ADC)) """ y = getproperty(x::Vector{ADC}, f::Symbol) diff --git a/KomaMRICore/src/datatypes/sequence/Delay.jl b/KomaMRICore/src/datatypes/sequence/Delay.jl index 03bcd778a..6ce9cbb46 100644 --- a/KomaMRICore/src/datatypes/sequence/Delay.jl +++ b/KomaMRICore/src/datatypes/sequence/Delay.jl @@ -56,3 +56,5 @@ sequence. """ +(s::Sequence, d::Delay) = s + Sequence([Grad(0.,d.T)]) +(d::Delay, s::Sequence) = Sequence([Grad(0.,d.T)]) + s + +Base.:(==)(x::Delay, y::Delay) = all(getfield(x, k) == getfield(y, k) for k ∈ fieldnames(Delay)) \ No newline at end of file diff --git a/KomaMRICore/src/datatypes/sequence/Grad.jl b/KomaMRICore/src/datatypes/sequence/Grad.jl index df9d73eba..d4a0ee3b2 100644 --- a/KomaMRICore/src/datatypes/sequence/Grad.jl +++ b/KomaMRICore/src/datatypes/sequence/Grad.jl @@ -187,6 +187,7 @@ Base.isapprox(gr1::Grad, gr2::Grad) = begin return all(length(getfield(gr1, k)) ≈ length(getfield(gr2, k)) for k ∈ fieldnames(Grad)) && all(getfield(gr1, k) ≈ getfield(gr2, k) for k ∈ fieldnames(Grad)) end +Base.:(==)(x::Grad, y::Grad) = all(getfield(x, k) == getfield(y, k) for k ∈ fieldnames(Grad)) # Gradient operations *(x::Grad,α::Real) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay) diff --git a/KomaMRICore/src/datatypes/sequence/RF.jl b/KomaMRICore/src/datatypes/sequence/RF.jl index 185e7decc..f68e743a2 100644 --- a/KomaMRICore/src/datatypes/sequence/RF.jl +++ b/KomaMRICore/src/datatypes/sequence/RF.jl @@ -102,6 +102,7 @@ Base.isapprox(rf1::RF, rf2::RF) = begin return all(length(getfield(rf1, k)) == length(getfield(rf2, k)) for k ∈ fieldnames(RF)) all(≈(getfield(rf1, k), getfield(rf2, k), atol=1e-9) for k ∈ fieldnames(RF)) end +Base.:(==)(x::RF, y::RF) = all(getfield(x, k) == getfield(y, k) for k ∈ fieldnames(RF)) # Properties size(r::RF, i::Int64) = 1 #To fix [r;r;;] concatenation of Julia 1.7.3 diff --git a/KomaMRICore/src/io/Pulseq.jl b/KomaMRICore/src/io/Pulseq.jl index a18a2bcb1..cf9cb7aba 100644 --- a/KomaMRICore/src/io/Pulseq.jl +++ b/KomaMRICore/src/io/Pulseq.jl @@ -82,11 +82,104 @@ function read_blocks(io, blockDurationRaster, version_combined) end end - r == NumberBlockEvents || break #Break on white space + r == NumberBlockEvents || break #Break on white space # fixme: pulseq doesn't technically forbid empty lines, does it? or comments... this would trip! although p8 or 9 of spec says "each subsequent line..." + end + sort(eventTable), sort(blockDurations), sort(delayIDs_tmp) # this is not correct! see docstring comment! +end +""" +read_blocks Read the [BLOCKS] section of a sequence file. + library=read_blocks(fid) Read blocks from file identifier of an + open MR sequence file and return the event table. + produces same output as read_blocks, but does not use scanf, + rather reads everything as int +""" +function read_blocks_split_instead_of_scanf(io, blockDurationRaster, version_combined) + eventTable = Dict() + blockDurations = Dict() + delayIDs_tmp = Dict() + if version_combined <= 1002001 + NumberBlockEvents = 7 + else + NumberBlockEvents = 8 + end + while true + blockEvents = [parse(Int, d) for d in eachsplit(readline(io))] + if length(blockEvents) == 0 + # empty line + break + end + if length(blockEvents) != NumberBlockEvents + error("Expected $NumberBlockEvents events but got $(length(blockEvents)).") + end + if blockEvents[1] != 0 + if version_combined <= 1002001 + eventTable[blockEvents[1]] = [0 blockEvents[3:end]... 0] + else + eventTable[blockEvents[1]] = [0 blockEvents[3:end]...] + end + + if version_combined >= 1004000 + blockDurations[blockEvents[1]] = blockEvents[2]*blockDurationRaster + else + delayIDs_tmp[blockEvents[1]] = blockEvents[2] + end + end end sort(eventTable), sort(blockDurations), sort(delayIDs_tmp) end +""" +read_blocks Read the [BLOCKS] section of a sequence file. + library=read_blocks(fid) Read blocks from file identifier of an + open MR sequence file and returns blockIDs, durations and delayIDs_tmp + as arrays or vectors, rather than dicts + +""" +function read_blocks_and_durs_as_arrays(io, blockDurationRaster, version_combined) + if version_combined <= 1002001 + NumberBlockEvents = 7 + else + NumberBlockEvents = 8 + end + # we'll collect everything into a vector and then reshape later + # we assume that we have at least 1000 blocks and allocate for that + blocks = empty!(Vector{Int}(undef, NumberBlockEvents * 1000)) # go through empty! to preallocate + blockDurations = empty!(Vector{Float64}(undef, 1000)) # go through empty! to preallocate + delayIDs_tmp = empty!(Vector{Float64}(undef, 1000)) # go through empty! to preallocate + num_lines = 0 + while true + blockEvents = [parse(Int, d) for d in eachsplit(readline(io))] + # todo: maybe better read chunks!? but then when do we stop? Maybe we should read entire file to ram before anyway... + if length(blockEvents) == 0 + # empty line. break (fixme: pulseq doesn't forbid that, I believe, but other funcs here rely on that too for detecting end of blocks block...) + # or maybe it forbids. p9 says 'subsequent lines...' + break + end + if length(blockEvents) != NumberBlockEvents + error("Expected $NumberBlockEvents events but got $(length(blockEvents)).") + end + if blockEvents[1] != 0 # id is not 0. Only then we care! + if version_combined <= 1002001 + # add an extra 0 at the end to make equal to other version + append!(blocks, blockEvents[3:end], 0) # only from 3 to end (discard id and duration, duration we treat below) + else + append!(blocks, blockEvents[3:end]) # only from 3 to end (discard id and duration, duration we treat below) + end + + if version_combined >= 1004000 + push!(blockDurations, blockEvents[2] * blockDurationRaster) + else + push!(delayIDs_tmp, blockEvents[2]) + end + end + num_lines += 1 + end + reshaped = reshape(blocks, :, num_lines) + # we need 6 vals because we drop IDs and durations. earlier versions we added a 0 -> for all it's 8 - 2 = 6 + @assert size(reshaped, 1) == 6 "unexpected number of fields per block" # for all versions! + reshaped, blockDurations, delayIDs_tmp +end + """ read_events Read an event section of a sequence file. library=read_events(fid) Read event data from file identifier of @@ -665,3 +758,219 @@ function get_block(obj, i) s = Sequence(G,R,A,D,E) s end + +""" +does what get_blocks does, but for many blocks at once +and works on array that contains block ids directly +it will use all blocks (2nd dim in blockEvents) to create the resulting sequence +""" +function get_seq_from_blocks(libraries, blockEvents::Array{Int, 2}, blockDurations::Vector{Float64}) + # get some views for blocks things + ids_rf = blockEvents[1, :] + ids_gradx = blockEvents[2, :] + ids_grady = blockEvents[3, :] + ids_gradz = blockEvents[4, :] + ids_adc = blockEvents[5, :] + ids_ext = blockEvents[6, :] + + # grads first + Δt_gr = libraries["definitions"]["GradientRasterTime"] + # allocate grads array space + num_blocks = size(blockEvents, 2) + GR = Array{Grad}(undef, 3, num_blocks) + GR[1, :] = collect(read_Grad(libraries["gradLibrary"], libraries["shapeLibrary"], Δt_gr, id) for id in ids_gradx) + GR[2, :] = collect(read_Grad(libraries["gradLibrary"], libraries["shapeLibrary"], Δt_gr, id) for id in ids_grady) + GR[3, :] = collect(read_Grad(libraries["gradLibrary"], libraries["shapeLibrary"], Δt_gr, id) for id in ids_gradz) + + #RFs + Δt_rf = libraries["definitions"]["RadiofrequencyRasterTime"] + # result of read_RF is a 1-element matrix... + RFs = empty!(Vector{RF}(undef, num_blocks)) + append!(RFs, read_RF(libraries["rfLibrary"], libraries["shapeLibrary"], Δt_rf, id)[1, 1] for id in ids_rf) + RFs = reshape(RFs, 1, num_blocks) + + # ADC + ADCs = empty!(Vector{ADC}(undef, num_blocks)) + # result of adc is a 1-element vector... + append!(ADCs, (read_ADC(libraries["adcLibrary"], id)[1] for id in ids_adc)) + + #DEFs which here are only Extensions. we can collect them in a vector + # and then set the dict key + DEFs = Dict("extension"=>ids_ext) + Sequence(GR, RFs, ADCs, blockDurations, DEFs) +end + +""" + + seq = read_seq_via_blocks_as_int_array(filename) + +Returns the Sequence struct from a Pulseq file with `.seq` extension. + +ALTERNATIVE IMPLEMENTATION OF read_seq. +This reads blocks as a big int array. Does not sort blocks by ID! +(pulseq specification says not to order IDs, although it's done in matlab!) +Calls costly Sequence constructor only once. + +# Arguments +- `filename`: (`::String`) absolute or relative path of the sequence file `.seq` + +# Returns +- `seq`: (`::Sequence`) Sequence struct + +# Examples +```julia-repl +julia> seq_file = joinpath(dirname(pathof(KomaMRI)), "../examples/1.sequences/spiral.seq") + +julia> seq = read_seq_via_blocks_as_int_array(seq_file) + +julia> plot_seq(seq) +``` +""" +function read_seq_via_blocks_as_int_array(filename) + @info "Loading sequence $(basename(filename)) ..." + version_combined = 0 + version_major = 0 + version_minor = 0 + gradLibrary = Dict() + def = Dict() + rfLibrary = Dict() + adcLibrary = Dict() + tmp_delayLibrary = Dict() + shapeLibrary = Dict() + extensionLibrary = Dict() + triggerLibrary = Dict() + blockEvents = Array{Int, 2}(undef, 0, 0) + blockDurations = Vector{Float64}(undef, 0) + delayInd_tmpVec = Vector{Int}(undef, 0) + #Reading file and storing data + open(filename) do io + while !eof(io) + section = readline(io) + if section == "[DEFINITIONS]" + def = read_definitions(io) + elseif section == "[VERSION]" + version_major, version_minor, _, version_combined = read_version(io) + elseif section == "[BLOCKS]" + if version_combined == 0 + @error "Pulseq file MUST include [VERSION] section prior to [BLOCKS] section" + end + blockEvents, blockDurations, delayInd_tmpVec = read_blocks_and_durs_as_arrays(io, def["BlockDurationRaster"], version_combined) + elseif section == "[RF]" + if version_combined >= 1004000 + rfLibrary = read_events(io, [1/γ 1 1 1 1e-6 1 1]) # this is 1.4.x format + else + rfLibrary = read_events(io, [1/γ 1 1 1e-6 1 1]) # this is 1.3.x and below + # we will have to scan through the library later after all the shapes have been loaded + end + elseif section == "[GRADIENTS]" + if version_combined >= 1004000 + gradLibrary = read_events(io, [1/γ 1 1 1e-6]; type='g', eventLibrary=gradLibrary) # this is 1.4.x format + else + gradLibrary = read_events(io, [1/γ 1 1e-6]; type='g', eventLibrary=gradLibrary) # this is 1.3.x and below + end + elseif section == "[TRAP]" + gradLibrary = read_events(io, [1/γ 1e-6 1e-6 1e-6 1e-6]; type='t', eventLibrary=gradLibrary); + elseif section == "[ADC]" + adcLibrary = read_events(io, [1 1e-9 1e-6 1 1]) + elseif section == "[DELAYS]" + if version_combined >= 1004000 + @error "Pulseq file revision 1.4.0 and above MUST NOT contain [DELAYS] section" + end + tmp_delayLibrary = read_events(io, 1e-6); + elseif section == "[SHAPES]" + shapeLibrary = read_shapes(io, (version_major==1 && version_minor<4)) + elseif section == "[EXTENSIONS]" + extensionLibrary = read_events(io,[1 1 1]) #For now, it reads the extensions but it does not take it them into account + elseif section == "extension TRIGGERS 1" + triggerLibrary = read_events(io,[1 1 1e-6 1e-6]) + elseif section == "[SIGNATURE]" + #Not implemented yet + end + + end + end + # fix blocks, gradients and RF objects imported from older versions + if version_combined < 1004000 + # scan through the RF objects + for i = 0:length(rfLibrary)-1 + rfLibrary[i]["data"] = [rfLibrary[i]["data"][1:3]' 0 rfLibrary[i]["data"][4:end]'] + end + # scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients) + for i = 0:length(gradLibrary)-1 + if gradLibrary[i]["type"] == 't' + #(1)amplitude (2)rise (2)flat (3)fall (4)delay + if gradLibrary[i]["data"][2] == 0 #rise + if abs(gradLibrary[i]["data"][1]) == 0 && gradLibrary[i]["data"][3] > 0 + gradLibrary[i]["data"][3] -= def["gradRasterTime"] + gradLibrary[i]["data"][2] = def["gradRasterTime"] + end + end + if gradLibrary[i]["data"][4] == 0 #delay + if abs(gradLibrary[i]["data"][1]) == 0 && gradLibrary[i]["data"][3] > 0 + gradLibrary[i]["data"][3] -= def["gradRasterTime"] + gradLibrary[i]["data"][4] = def["gradRasterTime"] + end + end + end + if gradLibrary[i]["type"] == 'g' + #(1)amplitude (2)amp_shape_id (3)time_shape_id (4)delay + gradLibrary[i]["data"] = [gradLibrary[i]["data"][1:2]; 0; gradLibrary[i]["data"][3:end]] + end + end + # for versions prior to 1.4.0 blockDurations have not been initialized + resize!(blockDurations, size(blockEvents, 2)) + for i = 1:size(blockEvents, 1) + idelay = delayInd_tmpVec[i] + if idelay > 0 + delay = tmp_delayLibrary[idelay]["data"][1] + blockDurations[i] = delay + end + end + end + #Sequence + libraries = Dict( + "gradLibrary"=>gradLibrary, + "rfLibrary"=>rfLibrary, + "adcLibrary"=>adcLibrary, + "tmp_delayLibrary"=>tmp_delayLibrary, + "shapeLibrary"=>shapeLibrary, + "extensionLibrary"=>extensionLibrary, + "triggerLibrary"=>triggerLibrary, + "definitions"=>def) + #Transforming Dictionary to Sequence object + #This should only work for Pulseq files >=1.4.0 + # in this most optimized implementation, the following lines still take + # 45% of all the time. These have to be optimized, as we still need + # 2s for a 2d 18point mrf seq with 64 PE lines. + # a 3d seq with 300 points and say 30x30 = 900 lines + # would take 234x longer -> 8min! + + seq = get_seq_from_blocks(libraries, blockEvents, blockDurations) + # old code + # seq = Sequence() + # for i = 1:length(blockEvents) + # seq += get_block(obj,i) + # end + # Final details + # Remove dummy seq block at the start, Issue #203 + # not needed anymore + # seq = seq[2:end] + # Hack for including extension and triggers + seq.DEF["additional_text"] = read_Extension(extensionLibrary, triggerLibrary) #Temporary hack + seq.DEF = recursive_merge(libraries["definitions"], seq.DEF) + # Koma specific details for reconstrucion + seq.DEF["FileName"] = basename(filename) + seq.DEF["PulseqVersion"] = version_combined + if !haskey(seq.DEF,"Nx") + Nx = maximum(seq.ADC.N) + RF_ex = (get_flip_angles(seq) .<= 90.01) .* is_RF_on.(seq) + Nz = max(length(unique(seq.RF[RF_ex].Δf)), 1) + Ny = sum(is_ADC_on.(seq)) / Nz |> x->floor(Int,x) + + seq.DEF["Nx"] = Nx #Number of samples per ADC + seq.DEF["Ny"] = Ny #Number of ADC events + seq.DEF["Nz"] = Nz #Number of unique RF frequencies, in a 3D acquisition this should not work + end + #Koma sequence + return seq +end \ No newline at end of file