Write Pulseq

.gitignore

@@ -34,4 +34,5 @@ examples/5.koma_paper/comparison_accuracy/**/*.mrd
 examples/5.koma_paper/mrf/**/*.mrd
 examples/5.koma_paper/comparison_accuracy/**/*.mat
 examples/5.koma_paper/comparison_accuracy/**/*.h5
-!examples/5.koma_paper/Manifest.toml
+!examples/5.koma_paper/Manifest.toml
+*_w.seq

KomaMRIBase/src/KomaMRIBase.jl

@@ -38,7 +38,7 @@ export Grad, RF, ADC, Delay
 export dur, get_block_start_times, get_samples
 export DiscreteSequence
 export discretize, get_adc_phase_compensation, get_adc_sampling_times
-export is_Gx_on, is_Gy_on, is_Gz_on, is_RF_on, is_ADC_on
+export is_Gx_on, is_Gy_on, is_Gz_on, is_RF_on, is_ADC_on, is_on
 export times, ampls, freqs
 # This are also used for simulation
 export kfoldperm, trapz, cumtrapz
@@ -55,8 +55,10 @@ export get_Mk, get_kspace, get_M0, get_M1, get_M2
 include("sequences/PulseDesigner.jl")
 export PulseDesigner
 
-#Package version, KomaMRIBase.__VERSION__
+# Package version, KomaMRIBase.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end # module KomaMRIBase

KomaMRIBase/src/datatypes/Sequence.jl

@@ -121,22 +121,35 @@
 -(x::Sequence) = Sequence(-x.GR, x.RF, x.ADC, x.DUR, x.DEF)
 *(x::Sequence, α::Real) = Sequence(α .* x.GR, x.RF, x.ADC, x.DUR, x.DEF)
 *(α::Real, x::Sequence) = Sequence(α .* x.GR, x.RF, x.ADC, x.DUR, x.DEF)
-*(x::Sequence, α::ComplexF64) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
-*(α::ComplexF64, x::Sequence) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
-*(x::Sequence, A::Matrix{Float64}) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
-*(A::Matrix{Float64}, x::Sequence) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
-/(x::Sequence, α::Real) = Sequence(x.GR/α, x.RF, x.ADC, x.DUR, x.DEF)
+*(x::Sequence, α::ComplexF64) = Sequence(x.GR, α .* x.RF, x.ADC, x.DUR, x.DEF)
+*(α::ComplexF64, x::Sequence) = Sequence(x.GR, α .* x.RF, x.ADC, x.DUR, x.DEF)
+*(x::Sequence, A::Matrix{Float64}) = Sequence(A * x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
+*(A::Matrix{Float64}, x::Sequence) = Sequence(A * x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
+/(x::Sequence, α::Real) = Sequence(x.GR / α, x.RF, x.ADC, x.DUR, x.DEF)
 #Grad operations
-+(s::Sequence, g::Grad) = s + Sequence(reshape([g],1,1)) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(g::Grad, s::Sequence) = Sequence(reshape([g],1,1)) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, g::Grad) = s + Sequence([g;;])
++(g::Grad, s::Sequence) = Sequence([g;;]) + s
 #RF operations
-+(s::Sequence, r::RF) = s + Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([r],1,1)) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(r::RF, s::Sequence) = Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([r],1,1)) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, r::RF) = s + Sequence([Grad(0.0, 0.0);;], [r;;])
++(r::RF, s::Sequence) = Sequence([Grad(0.0, 0.0);;], [r;;]) + s
 #ADC operations
-+(s::Sequence, a::ADC) = s + Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([RF(0.0,0.0)],1,1),[a]) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(a::ADC, s::Sequence) = Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([RF(0.0,0.0)],1,1),[a]) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, a::ADC) = s + Sequence([Grad(0.0, 0.0);;], [RF(0.0, 0.0);;], [a])
++(a::ADC, s::Sequence) = Sequence([Grad(0.0, 0.0);;], [RF(0.0, 0.0);;], [a]) + s
 #Sequence object functions
-size(x::Sequence) = size(x.GR[1,:])
+size(x::Sequence) = size(x.GR[1, :])
+
+"""
+For comparing two `Sequence`s custom types
+"""
+Base.isapprox(s1::Sequence, s2::Sequence) = begin
+    (length(s1) != length(s2)) && return false
+    return all([
+        s1.ADC[i] ≈ s2.ADC[i] &&
+        s1.RF[i] ≈ s2.RF[i] &&
+        s1.GR[i] ≈ s2.GR[i] &&
+        s1.DUR[i] ≈ s2.DUR[i] for i in 1:length(s1)
+    ])
+end
 
 """
     y = is_ADC_on(x::Sequence)

KomaMRIBase/src/datatypes/sequence/ADC.jl

@@ -39,10 +39,16 @@
     end
 end
 
-# ADC comparison
-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))
+"""
+For comparing two `ADC`s custom types
+"""
+function Base.isapprox(adc1::ADC, adc2::ADC)
+    for k in fieldnames(ADC)
+        if length(getfield(adc1, k)) != length(getfield(adc2, k))
+            return false
+        end
+    end
+    return all(getfield(adc1, k) ≈ getfield(adc2, k) for k in fieldnames(ADC))
 end
 
 """

KomaMRIBase/src/datatypes/sequence/Grad.jl

@@ -167,25 +167,31 @@
 - `y`: (`::Vector{Any}` or `::Matrix{Any}`) vector or matrix with the property defined
     by the symbol `f` for all elements of the Grad vector or matrix `x`
 """
-getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x,f)
+getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x, f)
 getproperty(x::Matrix{Grad}, f::Symbol) = begin
-	if f == :x
-		x[1,:]
-	elseif f == :y && size(x,1) >= 2
-		x[2,:]
-	elseif f == :z && size(x,1) >= 3
-		x[3,:]
-	elseif f == :dur
-		maximum(dur.(x), dims=1)[:]
-	else
-		getproperty.(x,f)
-	end
+    if f == :x
+        x[1, :]
+    elseif f == :y && size(x, 1) >= 2
+        x[2, :]
+    elseif f == :z && size(x, 1) >= 3
+        x[3, :]
+    elseif f == :dur
+        maximum(dur.(x); dims=1)[:]
+    else
+        getproperty.(x, f)
+    end
 end
 
-# Gradient comparison
-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))
+"""
+For comparing two `Grad`s custom types
+"""
+function Base.isapprox(gr1::Grad, gr2::Grad)
+    for k in fieldnames(Grad)
+        if length(getfield(gr1, k)) != length(getfield(gr2, k))
+            return false
+        end
+    end
+    return all(getfield(gr1, k) ≈ getfield(gr2, k) for k in fieldnames(Grad))
 end
 
 # Gradient operations

KomaMRIBase/src/datatypes/sequence/RF.jl

@@ -78,29 +78,35 @@
 - `y`: (`::Vector{Any}` or `::Matrix{Any}`) vector with the property defined by the
     symbol `f` for all elements of the RF vector or matrix `x`
 """
-getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x,f)
+getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x, f)
 getproperty(x::Matrix{RF}, f::Symbol) = begin
-	if     f == :Bx
-		real.(getproperty.(x,:A))
-	elseif f == :By
-		imag.(getproperty.(x,:A))
-	elseif f == :Δf
-		getproperty.(x,:Δf)
-	elseif f == :T || f == :delay
-		getproperty.(x[1,:],f)
-	elseif f == :dur
-		T, delay = x.T, x.delay
-		ΔT = [sum(t) for t=T] .+ delay
-		ΔT
-	else
-		getproperty.(x,f)
-	end
+    if f == :Bx
+        real.(getproperty.(x, :A))
+    elseif f == :By
+        imag.(getproperty.(x, :A))
+    elseif f == :Δf
+        getproperty.(x, :Δf)
+    elseif f == :T || f == :delay
+        getproperty.(x[1, :], f)
+    elseif f == :dur
+        T, delay = x.T, x.delay
+        ΔT = [sum(t) for t in T] .+ delay
+        ΔT
+    else
+        getproperty.(x, f)
+    end
 end
 
-# RF comparison
-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))
+"""
+For comparing two `RF`s custom types
+"""
+function Base.isapprox(rf1::RF, rf2::RF)
+    for k in fieldnames(RF)
+        if length(getfield(rf1, k)) != length(getfield(rf2, k))
+            return false
+        end
+    end
+    return all(≈(getfield(rf1, k), getfield(rf2, k), atol=1e-9) for k in fieldnames(RF))
 end
 
 # Properties

KomaMRICore/src/KomaMRICore.jl

@@ -29,8 +29,10 @@ export simulate, simulate_slice_profile
 # Spinors
 export Spinor, Rx, Ry, Rz, Q, Un
 
-#Package version, KomaMRICore.__VERSION__
+# Package version, KomaMRICore.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end

KomaMRIFiles/Project.toml

@@ -8,8 +8,10 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+MD5 = "6ac74813-4b46-53a4-afec-0b5dc9d7885c"
 MRIFiles = "5a6f062f-bf45-497d-b654-ad17aae2a530"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Scanf = "6ef1bc8b-493b-44e1-8d40-549aa65c4b41"
 
@@ -18,6 +20,7 @@ FileIO = "1"
 HDF5 = "0.16, 0.17"
 KomaMRIBase = "0.8"
 MAT = "0.10"
+MD5 = "0.2"
 MRIFiles = "0.1, 0.2, 0.3"
 Pkg = "1.4"
 Reexport = "1"

KomaMRIFiles/src/KomaMRIFiles.jl

@@ -2,22 +2,26 @@ module KomaMRIFiles
 
 using KomaMRIBase
 using Scanf, FileIO, HDF5, MAT    # IO related
+using Printf, MD5                 # For writing .seq files from Sequence struct
 
 using Reexport
 using MRIFiles
 import MRIFiles: insertNode
 @reexport using MRIFiles: ISMRMRDFile
 @reexport using FileIO: save
 
-include("Sequence/Pulseq.jl")
+include("Sequence/ReadPulseq.jl")
+include("Sequence/WritePulseq.jl")
 include("Phantom/JEMRIS.jl")
 include("Phantom/MRiLab.jl")
 
-export read_seq                                 # Pulseq
-export read_phantom_jemris, read_phantom_MRiLab # Phantom
+export read_seq, write_seq                       # Pulseq
+export read_phantom_jemris, read_phantom_MRiLab  # Phantom
 
 # Package version: KomaMRIFiles.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end # module KomaMRIFiles

KomaMRIFiles/src/Sequence/Pulseq.jl → KomaMRIFiles/src/Sequence/ReadPulseq.jl

@@ -520,7 +520,8 @@ function read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
             G = Grad(gA,gT,0.0,0.0,delay)
         end
     end
-    G
+    # If the amplitude of the gradient is zero, then return a zero gradient
+    return sum(abs.(G.A)) == 0 ? Grad(0.0, 0.0) : G
 end
 
 """