Pulse Designer

KomaMRIBase/Project.toml

@@ -11,11 +11,18 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
+[weakdeps]
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+
+[extensions]
+PulseDesignerUnitfulExt = "Unitful"
+
 [compat]
 Interpolations = "0.13, 0.14, 0.15"
 MAT = "0.10"
 MRIBase = "0.4"
 Parameters = "0.12"
 Pkg = "1.4"
 Reexport = "1"
+Unitful = "1.19"
 julia = "1.9"

KomaMRIBase/ext/PulseDesignerUnitfulExt.jl

@@ -0,0 +1,21 @@
+module PulseDesignerUnitfulExt
+
+using KomaMRIBase, Unitful
+
+Angle{T} = Union{Quantity{T, NoDims, typeof(u"rad")}, Quantity{T, NoDims, typeof(u"°")}} where T
+
+function PulseDesigner.block_pulse(flip_angle::Angle, duration::Unitful.Time;
+    phase_offset::Angle=0u"°", freq_offset::Unitful.Frequency=0u"Hz", delay::Unitful.Time=0u"s", sys=Scanner())
+    flip_angle, duration, phase_offset, freq_offset, delay = upreferred.((flip_angle, duration, phase_offset, freq_offset, delay)) .|> ustrip
+    flip_angle, duration = FlipAngle(flip_angle), Duration(duration)
+    return PulseDesigner.block_pulse(flip_angle, duration; phase_offset, freq_offset, delay, sys)
+end
+
+function PulseDesigner.block_pulse(flip_angle::Angle, bandwidth::Unitful.Frequency;
+    time_bw_product::DimensionlessQuantity=0.25u"Hz*s", phase_offset::Angle=0u"°", freq_offset::Unitful.Frequency=0u"Hz", delay::Unitful.Time=0u"s", sys=Scanner())
+    flip_angle, bandwidth, time_bw_product, phase_offset, freq_offset, delay = upreferred.((flip_angle, bandwidth, time_bw_product, phase_offset, freq_offset, delay)) .|> ustrip
+    flip_angle, bandwidth = FlipAngle(flip_angle), Bandwidth(bandwidth)
+    return PulseDesigner.block_pulse(flip_angle, bandwidth; time_bw_product, phase_offset, freq_offset, delay, sys)
+end
+
+end

KomaMRIBase/src/KomaMRIBase.jl

@@ -52,6 +52,13 @@ export get_flip_angles, is_RF_on, is_GR_on, is_ADC_on
 export get_M0, get_M1, get_M2, get_kspace
 
 # PulseDesigner submodule
+struct FlipAngle val end
+struct Duration val end
+struct Bandwidth val end
+struct TimeBwProduct val end
+struct SliceThickness val end
+export FlipAngle, Duration, Bandwidth, TimeBwProduct, SliceThickness
+
 include("sequences/PulseDesigner.jl")
 export PulseDesigner
 

KomaMRIBase/src/sequences/ADC/design.jl

@@ -0,0 +1,26 @@
+"""
+"""
+function make_adc(num_samles; Δtadc=0, duration=0, delay=0,
+    freq_offset=0, phase_offset=0, dead_time=0, sys=nothing)
+
+    if !isnothing(sys)
+        dead_time = sys.ADC_dead_time_T
+        Δtadc = sys.ADC_Δt
+    end
+
+    if (Δtadc == 0 && duration == 0) || (Δtadc > 0 && duration > 0)
+        @error "Either dwell or duration must be defined"
+    end
+    if duration > 0
+        Δtadc = duration / num_samles
+    elseif Δtadc > 0
+        duration = Δtadc * num_samles;
+    end
+    if dead_time > delay
+        delay = dead_time; # adcDeadTime is added before the actual sampling (and also second time after the sampling period)
+    end
+
+    adc = ADC(num_samles, duration, delay, freq_offset, phase_offset)
+
+    return adc
+end

KomaMRIBase/src/sequences/Grad/design.jl

@@ -0,0 +1,116 @@
+"""
+"""
+function trapezoid(; duration=0, amplitude=0, area=0, flat_area=0,
+    flat_time=0, rise_time=0, fall_time=0, delay=0,
+    max_grad=Inf, max_slew=Inf, Δtgr=1e-3, sys=nothing)
+
+    if !isnothing(sys)
+        max_grad = sys.Gmax
+        max_slew = sys.Smax
+        Δtgr = sys.GR_Δt
+    end
+
+    if area == 0 && flat_area == 0 && amplitude == 0
+        @error "trapezoid: invalid keywords. Must supply either 'area', 'flat_area' or 'amplitude'"
+    end
+    if fall_time > 0 && rise_time == 0
+        @error "trapezoid: invalid keywords. Must always supply 'rise_time' if 'fall_time' is specified explicitly."
+    end
+
+    if flat_time > 0
+        if amplitude == 0
+            if flat_area == 0
+                @error "trapezoid: invalid keyworks. When 'flat_time' is provided either 'flat_area' or 'amplitude' must be provided as well; you may consider providing 'duration', 'area' and optionally ramp times instead."
+            end
+            amplitude = flat_area / flat_time
+        end
+        if rise_time == 0
+            rise_time = abs(amplitude) / max_slew;
+            rise_time = ceil(rise_time / Δtgr) * Δtgr;
+            if rise_time == 0
+                rise_time = Δtgr
+            end
+        end
+        if fall_time == 0
+            fall_time = rise_time
+        end
+    elseif duration > 0
+        if amplitude == 0
+            if rise_time == 0
+                dC = 1 / abs(2 * max_slew) + 1 / abs(2 * max_slew)
+                possible = duration^2 > 4 * abs(area) * dC;
+                @assert possible "Requested area is too large for this gradient. Minimum required duration (assuming triangle gradient can be realized) is $(round(sqrt(4 * abs(area) * dC) * 1e6)) us"
+                amplitude = (duration - sqrt(duration^2 - 4 * abs(area) * dC)) / (2 * dC)
+            else
+                if fall_time == 0
+                    fall_time = rise_time
+                end
+                amplitude = area / (duration - 0.5 * rise_time - 0.5 * fall_time)
+                possible = duration > (rise_time + fall_time) && abs(amplitude) < max_grad
+                @assert possible "Requested area is too large for this gradient. Probably amplitude is violated ($(round(abs(amplitude) / max_grad * 100))%)"
+            end
+        end
+        if rise_time == 0
+            rise_time = ceil(abs(amplitude) / max_slew / Δtgr) * Δtgr
+            if rise_time == 0
+                rise_time = Δtgr
+            end
+        end
+        if fall_time == 0
+            fall_time = rise_time
+        end
+        flat_time = duration - rise_time - fall_time
+        if amplitude == 0
+            # Adjust amplitude (after rounding) to achieve given area
+            amplitude = area / (rise_time / 2 + fall_time / 2 + flat_time)
+        end
+    else
+        if area == 0
+            @error "trapezoid: invalid keywords. Must supply area or duration"
+        else
+            # find the shortest possible duration
+            # first check if the area can be realized as a triangle
+            # if not we calculate a trapezoid
+            rise_time = ceil(sqrt(abs(area) / max_slew) / Δtgr) * Δtgr
+            if rise_time < Δtgr  # the "area" was probably 0 or almost 0 ...
+                rise_time = Δtgr;
+            end
+            amplitude = area / rise_time
+            t_eff = rise_time
+            if abs(amplitude) > max_grad
+                t_eff = ceil(abs(area) / max_grad / Δtgr) * Δtgr
+                amplitude = area / t_eff
+                if rise_time == 0
+                    rise_time = Δtgr
+                end
+            end
+            flat_time = t_eff - rise_time
+            fall_time = rise_time
+        end
+    end
+
+    @assert abs(amplitude) <= max_grad "trapezoid: invalid amplitude. Amplitude violation ($(round(abs(amplitude) / max_grad * 100))%)"
+
+    return Grad(amplitude, flat_time, rise_time, fall_time, delay)
+end
+
+"""
+"""
+function arbitrary_grad(waveform; delay=0,
+    max_grad=Inf, max_slew=Inf, Δtgr=1e-3, sys=nothing)
+
+    if !isnothing(sys)
+        max_grad = sys.Gmax
+        max_slew = sys.Smax
+        Δtgr = sys.GR_Δt
+    end
+
+    slew = (waveform[2:end] - waveform[1:end-1]) / Δtgr
+    if !isempty(slew)
+        @assert maximum(abs.(slew)) <= max_slew "Slew rate violation ($(maximum(abs.(slew)) / max_slew * 100)%)"
+    end
+    @assert maximum(abs.(waveform)) <= max_grad "Gradient amplitude violation ($(maximum(abs.(waveform)) / max_grad * 100)%)"
+
+    duration = (length(waveform)-1) * Δtgr
+    return Grad(waveform, duration, 0, 0, delay)
+end

KomaMRIBase/src/sequences/PulseDesigner.jl

@@ -6,6 +6,11 @@ A module to define different pulse sequences.
 module PulseDesigner
 using ..KomaMRIBase
 
+include("RF/block_pulse.jl")
+include("RF/sinc_pulse.jl")
+include("Grad/design.jl")
+include("ADC/design.jl")
+
 """
     seq = RF_hard(B1, T, sys; G=[0, 0, 0], Δf=0)
 

KomaMRIBase/src/sequences/RF/block_pulse.jl

@@ -0,0 +1,19 @@
+
+function block_pulse(flip_angle::FlipAngle, duration::Duration;
+    phase_offset=0, freq_offset=0, delay=0, sys=Scanner())
+
+    flip_angle, duration = flip_angle.val, duration.val
+    amplitude = flip_angle / (2π * γ * duration) * exp(im * phase_offset)
+    delay = max(delay, sys.RF_dead_time_T)
+    block_duration = delay + duration + sys.RF_ring_down_T
+
+    rf = RF(amplitude, duration, freq_offset, delay)
+    return Sequence([Grad(0, 0);;], [rf;;], [ADC(0, 0)], block_duration)
+end
+
+function block_pulse(flip_angle::FlipAngle, bandwidth::Bandwidth; time_bw_product=0.25, kwargs...)
+    duration = Duration(time_bw_product / bandwidth.val)
+    return block_pulse(flip_angle, duration; kwargs...)
+end
+
+export block_pulse

KomaMRIBase/src/sequences/RF/design.jl

@@ -0,0 +1,118 @@
+"""
+"""
+function sinc_pulse(flip_angle; duration=0, freq_offset=0, phase_offset=0,
+    time_bw_product=0, apodization=0.5, centerpos=0.5, delay=0, slice_thickness=0,
+    dead_time=0, ring_down_time=0, Δtrf=1e-5, Δtgr=1e-3, sys=nothing)
+
+    if !isnothing(sys)
+        dead_time = sys.RF_dead_time_T
+        ring_down_time = sys.RF_ring_down_T
+        Δtrf = sys.RF_Δt
+        Δtgr = sys.GR_Δt
+    end
+
+    if duration <= 0
+        @error "rf pulse duration must be positive"
+    end
+    if Δtrf <= 0
+        @error "the Δtrf gradient raster time must be positive"
+    end
+
+    BW = time_bw_product / duration
+    N = Integer(ceil(duration / Δtrf))
+    t = range(0, duration; length=N)
+    window = (1 - apodization) .+ apodization * cos.(2π * ((t  .- (centerpos * duration)) / duration))
+    signal = window .* sinc.(BW * (t  .- (centerpos * duration)))
+    flip = 0.5 * sum(signal[2:end] + signal[1:end-1]) * Δtrf * 2π
+    signal = signal * flip_angle / flip * cis(phase_offset)
+    if dead_time > delay
+        delay = dead_time
+    end
+	rf = RF(signal, duration, freq_offset, delay)
+
+    @assert slice_thickness > 0 "slice_thickness must be provided"
+
+    amplitude = BW / slice_thickness
+    area = amplitude * duration
+    gz = trapezoid(; flat_time=duration, flat_area=area, sys=sys);
+    gz_area = gz.A * (gz.T + gz.rise / 2 + gz.fall / 2)
+    gzr_area = -area*(1 - centerpos) - 0.5*(gz_area - area)
+    gzr = trapezoid(; sys=sys, area=gzr_area)
+    if rf.delay > gz.rise
+        gz.delay = ceil((rf.delay - gz.rise) / Δtgr) * Δtgr     # round-up to gradient raster
+    end
+    if rf.delay < gz.rise + gz.delay
+        rf.delay = gz.rise + gz.delay   # these are on the grad raster already which is coarser
+    end
+
+    dly = Delay(0)
+    if ring_down_time > 0
+        dly = Delay(rf.delay + rf.T + ring_down_time)     # I NEED a review
+    end
+
+    return rf, gz, gzr, dly
+end
+
+"""
+"""
+function arbitrary_rf(signal, flip; freq_offset=0, phase_offset=0,
+    time_bw_product=0, bandwidth=0, delay=0, slice_thickness=0,
+    dead_time=0, ring_down_time=0, Δtrf=1e-5, Δtgr=1e-3, sys=nothing)
+
+    if !isnothing(sys)
+        dead_time = sys.RF_dead_time_T
+        ring_down_time = sys.RF_ring_down_T
+        Δtrf = sys.RF_Δt
+        Δtgr = sys.GR_Δt
+    end
+
+    signal = signal / abs(sum(signal * Δtrf)) * flip / 2π * cis(phase_offset)
+
+    N =  length(signal)
+    duration = (N-1) * Δtrf
+    t = range(0, duration; length=N)
+
+    if dead_time > delay
+        delay = dead_time;
+    end
+
+    rf = RF(signal, duration, freq_offset, delay)
+
+    if time_bw_product > 0
+        if bandwidth > 0
+            @error "Both 'bandwidth' and 'time_bw_product' cannot be specified at the same time"
+        else
+            bandwidth = time_bw_product / duration
+        end
+    end
+
+    @assert slice_thickness > 0 "SliceThickness must be provided"
+    @assert bandwidth > 0 "Bandwidth of pulse must be provided"
+
+    BW = bandwidth
+    if time_bw_product > 0
+        BW = time_bw_product / duration
+    end
+
+    amplitude = BW / slice_thickness
+    area = amplitude * duration
+    gz = trapezoid(; flat_time=duration, flat_area=area, sys=sys)
+    gz_area = gz.A * (gz.T + gz.rise / 2 + gz.fall / 2)
+    gzr_area = -area*(1 - KomaMRIBase.get_RF_center(rf) / rf.T) - 0.5*(gz_area - area)
+    gzr = trapezoid(; sys=sys, area=gzr_area)
+
+
+    if rf.delay > gz.rise
+        gz.delay = ceil((rf.delay - gz.rise) / Δtgr) * Δtgr     # round-up to gradient raster
+    end
+    if rf.delay < gz.rise + gz.delay
+        rf.delay = gz.rise + gz.delay   # these are on the grad raster already which is coarser
+    end
+
+    dly = Delay(0)
+    if ring_down_time > 0
+        dly = Delay(rf.delay + rf.T + ring_down_time)     # I NEED a review
+    end
+
+    return rf, gz, gzr, delay
+end

KomaMRIBase/src/sequences/RF/sinc_pulse.jl

@@ -0,0 +1,24 @@
+function sinc_pulse(flip_angle::FlipAngle, duration::Duration, time_bw_product::TimeBwProduct;
+    phase_offset=0, freq_offset=0, delay=0, sys=Scanner(), apodization=0.5, centerpos=0.5)
+
+    flip_angle, duration, time_bw_product = flip_angle.val, duration.val, time_bw_product.val
+
+    dead_time = sys.RF_dead_time_T
+    ring_down_time = sys.RF_ring_down_T
+    Δtrf = sys.RF_Δt
+
+    BW = time_bw_product / duration
+    N = Integer(ceil(duration / Δtrf))
+    t = range(0, duration; length=N)
+    window = (1 - apodization) .+ apodization * cos.(2π * ((t  .- (centerpos * duration)) / duration))
+    signal = window .* sinc.(BW * (t  .- (centerpos * duration)))
+    flip = 0.5 * sum(signal[2:end] + signal[1:end-1]) * Δtrf * 2π
+    signal = signal * flip_angle / flip * exp(im * phase_offset)
+    delay = max(delay, dead_time)
+    block_duration = delay + duration + ring_down_time
+
+    rf = RF(signal, duration, freq_offset, delay)
+    return Sequence([Grad(0, 0);;], [rf;;], [ADC(0, 0)], block_duration)
+end
+
+export sinc_pulse