Add Pulseq ROTATIONS support via extension hooks

KomaMRIBase/src/KomaMRIBase.jl

@@ -42,6 +42,7 @@ export Scanner, Sequence, Phantom
 export Grad, RF, ADC, Delay
 export dur, get_block_start_times, get_samples
 export RFuse, Excitation, Refocusing, Inversion, Saturation, Preparation, Other, Undefined
+export apply_rotations!
 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

KomaMRIBase/src/datatypes/sequence/EXT.jl

@@ -2,11 +2,13 @@ include("extensions/AdcLabels.jl")
 export AdcLabels
 
 abstract type Extension end
-# Supported extensions. To add a new extension, create a new file in the extensions folder and add it to the list below.
+
+# Supported extensions. To add a new extension, create a new file in the extensions folder and add it below.
 include("extensions/LabelInc.jl")
 include("extensions/LabelSet.jl")
 include("extensions/Trigger.jl")
+include("extensions/QuaternionRot.jl")
 
 get_EXT_type_from_symbol(::Val) = nothing
 
-export Extension, LabelInc, LabelSet, Trigger, SoftDelay
+export Extension, LabelInc, LabelSet, Trigger, SoftDelay, QuaternionRot, apply_rotations!

KomaMRIBase/src/datatypes/sequence/extensions/QuaternionRot.jl

@@ -0,0 +1,41 @@
+mutable struct QuaternionRot <: Extension
+    q0::Float64
+    qx::Float64
+    qy::Float64
+    qz::Float64
+end
+
+get_scale(::Type{QuaternionRot}) = [1.0 1.0 1.0 1.0]
+get_scanf_format(::Type{QuaternionRot}) = "%f %f %f %f"
+get_EXT_type_from_symbol(::Val{:ROTATIONS}) = QuaternionRot
+Base.transpose(q::QuaternionRot) = QuaternionRot(q.q0, -q.qx, -q.qy, -q.qz)
+Base.adjoint(q::QuaternionRot) = transpose(q)
+
+function Base.:*(q::QuaternionRot, grads)
+    @assert length(grads) == 3 "A sequence block must contain x, y and z gradient channels."
+
+    qnorm = sqrt(abs2(q.q0) + abs2(q.qx) + abs2(q.qy) + abs2(q.qz))
+    qnorm > 0 || throw(ArgumentError("Rotation quaternion must have non-zero norm."))
+
+    q0, qx, qy, qz = q.q0 / qnorm, q.qx / qnorm, q.qy / qnorm, q.qz / qnorm
+
+    rot = [
+        (1 - 2 * (qy^2 + qz^2))   2 * (qx * qy - q0 * qz)   2 * (qx * qz + q0 * qy)
+        2 * (qx * qy + q0 * qz)   (1 - 2 * (qx^2 + qz^2))   2 * (qy * qz - q0 * qx)
+        2 * (qx * qz - q0 * qy)   2 * (qy * qz + q0 * qx)   (1 - 2 * (qx^2 + qy^2))
+    ]
+    return rot * grads
+end
+
+function apply_rotations!(seq; reverse=false)
+    for b in eachindex(seq.EXT)
+        block_ext = seq.EXT[b]
+        isempty(block_ext) && continue
+
+        for ext in block_ext
+            ext isa QuaternionRot || continue
+            seq.GR[:, b] = (reverse ? ext' : ext) * seq.GR[:, b]
+        end
+    end
+    return seq
+end

KomaMRIBase/test/runtests.jl

@@ -282,6 +282,35 @@ using TestItems, TestItemRunner
         PHS_vec = [l[i].PHS for i in eachindex(l)] 
         @test PHS_vec == vec([0 0 0 2 2 2])
 
+        # Pulseq Rotation extension support
+        rot = QuaternionRot(cos(π / 8), 0.0, 0.0, sin(π / 8)) # 45 deg around z
+        grads_rot = rot * [Grad(1.0, 1e-3), Grad(0.0, 1e-3), Grad(0.0, 1e-3)]
+        @test grads_rot[1].A ≈ cos(π / 4)
+        @test grads_rot[2].A ≈ sin(π / 4)
+        @test grads_rot[3].A ≈ 0.0
+        grads_unrot = transpose(rot) * grads_rot
+        @test grads_unrot[1].A ≈ 1.0
+        @test grads_unrot[2].A ≈ 0.0
+        @test grads_unrot[3].A ≈ 0.0
+
+        ext90 = QuaternionRot(cos(π / 4), 0.0, 0.0, sin(π / 4)) # 90 deg around z
+
+        seq_rot = Sequence([Grad(1.0, 1e-3); Grad(0.0, 1e-3); Grad(0.0, 1e-3);;])
+        seq_rot.EXT = [[ext90]]
+        apply_rotations!(seq_rot)
+        @test seq_rot.GR[1, 1].A ≈ 0.0 atol = 1e-12
+        @test seq_rot.GR[2, 1].A ≈ 1.0 atol = 1e-12
+        @test seq_rot.GR[3, 1].A ≈ 0.0 atol = 1e-12
+        @test length(seq_rot.EXT[1]) == 1 && seq_rot.EXT[1][1] isa QuaternionRot
+
+        seq_rot_rev = Sequence([Grad(1.0, 1e-3); Grad(0.0, 1e-3); Grad(0.0, 1e-3);;])
+        seq_rot_rev.EXT = [[ext90]]
+        apply_rotations!(seq_rot_rev; reverse=true)
+        @test seq_rot_rev.GR[1, 1].A ≈ 0.0 atol = 1e-12
+        @test seq_rot_rev.GR[2, 1].A ≈ -1.0 atol = 1e-12
+        @test seq_rot_rev.GR[3, 1].A ≈ 0.0 atol = 1e-12
+        @test length(seq_rot_rev.EXT[1]) == 1 && seq_rot_rev.EXT[1][1] isa QuaternionRot
+
     end
 
     @testset "DiscreteSequence" begin

KomaMRIFiles/src/Sequence/Pulseq.jl

@@ -356,12 +356,13 @@ function simplify_waveforms(amp_shape, time_shape)
 end
 
 """
-    seq = read_seq(filename)
+    seq = read_seq(filename; apply_rotations=true)
 
 Returns the Sequence struct from a Pulseq file with `.seq` extension.
 
 # Arguments
 - `filename`: (`::String`) absolute or relative path of the sequence file `.seq`
+- `apply_rotations`: (`::Bool`) apply `ROTATIONS` extensions to gradients after loading
 
 # Returns
 - `seq`: (`::Sequence`) Sequence struct
@@ -375,7 +376,7 @@ julia> seq = read_seq(seq_file)
 julia> plot_seq(seq)
 ```
 """
-function read_seq(filename)
+function read_seq(filename; apply_rotations=true)
     @info "Loading sequence $(basename(filename)) ..."
     pulseq_version = v"0.0.0"
     gradLibrary = Dict()
@@ -509,6 +510,9 @@ function read_seq(filename)
     seq.DEF["Nx"] = trunc(Int64, get(seq.DEF, "Nx", maximum(adc.N for adc = seq.ADC)))
     seq.DEF["Nz"] = trunc(Int64, get(seq.DEF, "Nz", length(unique(seq.RF.Δf))))
     seq.DEF["Ny"] = trunc(Int64, get(seq.DEF, "Ny", sum(map(is_ADC_on, seq)) ÷ seq.DEF["Nz"]))
+    if apply_rotations
+        KomaMRIBase.apply_rotations!(seq)
+    end
     #Koma sequence
     return seq
 end

KomaMRIFiles/test/runtests.jl

@@ -88,6 +88,24 @@ end
         @test seq.DEF["PulseqVersion"] == v"1.2.1"
         @test seq.DEF["signature"][:hash] == "e827cfff4436b65a6341a4fa0f6deb07"
     end
+    @testset "Rotations" begin
+        pth = (@__DIR__ ) * "/test_files/pulseq/basic_tests/v1.5/lps.seq"
+        seq_no_rot = @suppress read_seq(pth; apply_rotations=false)
+        seq_rot    = @suppress read_seq(pth) # default apply_rotations=true
+
+        nrot_no_rot = count(ext -> ext isa QuaternionRot, Iterators.flatten(seq_no_rot.EXT))
+        nrot_rot    = count(ext -> ext isa QuaternionRot, Iterators.flatten(seq_rot.EXT))
+        @test nrot_no_rot > 0
+        @test nrot_no_rot == nrot_rot
+
+        seq_no_rot_applied = copy(seq_no_rot)
+        apply_rotations!(seq_no_rot_applied)
+        @test seq_no_rot_applied.GR.A ≈ seq_rot.GR.A
+
+        seq_rot_reversed = copy(seq_rot)
+        apply_rotations!(seq_rot_reversed; reverse=true)
+        @test seq_rot_reversed.GR.A ≈ seq_no_rot.GR.A
+    end
     @testset "Compression-Decompression" begin
         shape = ones(100)
         num_samples, compressed_data = KomaMRIFiles.compress_shape(shape)
@@ -134,4 +152,4 @@ end
             end
         end
     end
-end
+end