Load pulseq sequences faster and allow comparison

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

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}})

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)

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))

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)

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

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