Pulse Designer

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,10 @@ A module to define different pulse sequences.
 module PulseDesigner
 using ..KomaMRIBase
 
+include("Grad/design.jl")
+include("RF/design.jl")
+include("ADC/design.jl")
+
 """
     seq = RF_hard(B1, T, sys; G=[0, 0, 0], Δf=0)
 

KomaMRIBase/src/sequences/RF/design.jl

@@ -0,0 +1,152 @@
+"""
+"""
+function block_pulse(flip_angle; duration=0, freq_offset=0, phase_offset=0,
+    time_bw_product=0, bandwidth=0, delay=0,
+    dead_time=0, ring_down_time=0, sys=nothing)
+
+    if !isnothing(sys)
+        dead_time = sys.RF_dead_time_T
+        ring_down_time = sys.RF_ring_down_T
+    end
+
+    if duration == 0
+        if time_bw_product > 0
+            duration = time_bw_product / bandwidth
+        elseif bandwidth > 0
+            duration = 1 / (4 * bandwidth)
+        else
+            @error "Either bandwidth or duration must be defined"
+        end
+    end
+    amplitude = flip_angle / (2π * γ * duration) * cis(phase_offset)
+
+    if dead_time > delay
+        delay = dead_time
+    end
+
+    dly = Delay(0)
+    if ring_down_time > 0
+        dly = Delay(delay + duration + ring_down_time)     # I NEED a review
+    end
+
+    return RF(amplitude, duration, freq_offset, delay), dly
+end
+
+"""
+"""
+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/test/runtests.jl

@@ -347,6 +347,60 @@ end
         spoke = PulseDesigner.radial_base(FOV, N, sys)
         @test spoke.DEF["Δθ"] ≈ π / Nspokes
     end
+    @testset "design" begin
+        # trapezoid
+        duration = 10
+        amplitude = 5
+        grad = PulseDesigner.trapezoid(; duration, amplitude)
+        seq = Sequence()
+        seq += grad
+        @test true
+
+        # arbitrary_grad
+        t = range(0, 1; length=20)
+        waveform = sin.(2π * t)
+        grad = PulseDesigner.arbitrary_grad(waveform)
+        seq = Sequence()
+        seq += grad
+        @test true
+
+        # block_pulse
+        flip_angle = π / 2
+        duration = 0.001
+        rf, delay = PulseDesigner.block_pulse(flip_angle; duration)
+        seq = Sequence()
+        seq += rf
+        @test true
+
+        # sinc_pulse
+        flip_angle = π / 2
+        duration = 0.001
+        slice_thickness = 0.001
+        time_bw_product = 4
+        rf, gz, gzr, dly = PulseDesigner.sinc_pulse(flip_angle; duration, time_bw_product, slice_thickness)
+        seq = Sequence([Grad(0,0); Grad(0,0); gz;;], [rf;;])
+        seq += Sequence([Grad(0,0); Grad(0,0); gzr;;])
+        @test true
+
+        # arbitrary_rf
+        flip = π / 2
+        slice_thickness = 0.001
+        time_bw_product = 4
+        t = range(0, 0.5; length=20)
+        signal = sin.(2π * t)
+        rf, gz, gzr, dly = PulseDesigner.arbitrary_rf(signal, flip; time_bw_product, slice_thickness)
+        seq = Sequence([Grad(0,0); Grad(0,0); gz;;], [rf;;])
+        seq += Sequence([Grad(0,0); Grad(0,0); gzr;;])
+        @test true
+
+        # make_adc
+        num_samles = 10
+        duration = 0.001
+        adc = PulseDesigner.make_adc(num_samles; duration)
+        seq = Sequence()
+        seq += adc
+        @test true
+    end
 end
 
 @testitem "Phantom" tags=[:base] begin