Pulse Designer by beorostica · Pull Request #305 · JuliaHealth/KomaMRI.jl

beorostica · 2024-02-02T21:34:31Z

No description provided.

cncastillo · 2024-02-02T21:36:16Z

KomaMRIBase/src/sequences/ADC/design.jl

@@ -0,0 +1,26 @@
+"""
+"""
+function make_adc(num_samles; Δtadc=0, duration=0, delay=0,


typo num_samles to num_samples

This function will be completelly replaced

cncastillo · 2024-02-02T21:37:58Z

KomaMRIBase/src/sequences/ADC/design.jl

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


it should output a Sequence

This function will be completelly replaced, all the outpus will be of type Sequence

cncastillo · 2024-02-02T21:38:05Z

KomaMRIBase/src/sequences/Grad/design.jl

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


it should output a Sequence

This function will be completelly replaced, all the outpus will be of type Sequence

cncastillo · 2024-02-02T21:38:12Z

KomaMRIBase/src/sequences/RF/design.jl

+        dly = Delay(delay + duration + ring_down_time)     # I NEED a review
+    end
+
+    return RF(amplitude, duration, freq_offset, delay), dly


It should output a Sequence. The DUR variable is not set using Delay objects, we may need to change that.

The idea now is to add a duration to the ouput Sequence

seq.DUR needs to be a vector of Delay's. To be able to combine them as we talked #4

cncastillo · 2024-02-02T21:38:20Z

KomaMRIBase/src/sequences/RF/design.jl

+        dly = Delay(rf.delay + rf.T + ring_down_time)     # I NEED a review
+    end
+
+    return rf, gz, gzr, dly


It should output a Sequence. The Delay should be the DUR.

cncastillo · 2024-02-02T21:38:25Z

KomaMRIBase/src/sequences/RF/design.jl

+        dly = Delay(rf.delay + rf.T + ring_down_time)     # I NEED a review
+    end
+
+    return rf, gz, gzr, delay


it should output a Sequence

cncastillo · 2024-02-02T21:39:52Z

KomaMRIBase/test/runtests.jl

The tests should verify that the gradient/rf/adc satisfies the design inputs (area of the gradiet, duration, etc)

cncastillo · 2024-02-02T21:42:07Z

KomaMRIBase/src/sequences/ADC/design.jl

+    if !isnothing(sys)
+        dead_time = sys.ADC_dead_time_T
+        Δtadc = sys.ADC_Δt
+    end


This should not be required, set the default for optional parameter sys = Scanner()

cncastillo · 2024-02-02T21:42:48Z

KomaMRIBase/src/sequences/ADC/design.jl

+    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


use multiple dispatch

cncastillo · 2024-02-02T21:43:00Z

KomaMRIBase/src/sequences/Grad/design.jl

+    if !isnothing(sys)
+        max_grad = sys.Gmax
+        max_slew = sys.Smax
+        Δtgr = sys.GR_Δt
+    end


This should not be required, set the default for optional parameter sys = Scanner()

cncastillo · 2024-02-02T21:43:29Z

KomaMRIBase/src/sequences/Grad/design.jl

+    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


use multiple dispatch to simplify this, the code is too convoluted

cncastillo · 2024-02-02T21:43:53Z

KomaMRIBase/src/sequences/Grad/design.jl

+        end
+    end
+
+    @assert abs(amplitude) <= max_grad "trapezoid: invalid amplitude. Amplitude violation ($(round(abs(amplitude) / max_grad * 100))%)"


Use functions to check hardware limits for the input sys.

cncastillo · 2024-02-02T21:44:31Z

KomaMRIBase/src/sequences/Grad/design.jl

+    if !isnothing(sys)
+        max_grad = sys.Gmax
+        max_slew = sys.Smax
+        Δtgr = sys.GR_Δt
+    end


Not required

cncastillo · 2024-02-02T21:44:48Z

KomaMRIBase/src/sequences/Grad/design.jl

+    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)%)"


Use functions to check this.

cncastillo · 2024-02-02T21:45:58Z

KomaMRIBase/src/sequences/RF/design.jl

+    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


mulitiple dispatch

cncastillo · 2024-02-02T21:46:05Z

KomaMRIBase/src/sequences/RF/design.jl

+    if !isnothing(sys)
+        dead_time = sys.RF_dead_time_T
+        ring_down_time = sys.RF_ring_down_T
+    end


Not required

codecov · 2024-02-02T21:47:22Z

Codecov Report

Attention: 57 lines in your changes are missing coverage. Please review.

Comparison is base (6469a10) 92.60% compared to head (133a7b1) 90.77%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #305      +/-   ##
==========================================
- Coverage   92.60%   90.77%   -1.84%     
==========================================
  Files          33       36       +3     
  Lines        2245     2417     +172     
==========================================
+ Hits         2079     2194     +115     
- Misses        166      223      +57

Flag	Coverage Δ
base	`88.19% <66.66%> (-4.49%)`	⬇️
core	`90.47% <ø> (ø)`
files	`92.34% <ø> (ø)`
komamri	`93.89% <ø> (ø)`
plots	`93.28% <ø> (ø)`

Flags with carried forward coverage won't be shown. Click here to find out more.

Files	Coverage Δ
KomaMRIBase/src/sequences/PulseDesigner.jl	`94.11% <ø> (+0.05%)`	⬆️
KomaMRIBase/src/sequences/ADC/design.jl	`57.14% <57.14%> (ø)`
KomaMRIBase/src/sequences/Grad/design.jl	`63.76% <63.76%> (ø)`
KomaMRIBase/src/sequences/RF/design.jl	`70.45% <70.45%> (ø)`

cncastillo · 2024-02-02T21:48:03Z

KomaMRIBase/src/sequences/RF/design.jl

+    if dead_time > delay
+        delay = dead_time
+    end


delay = max(dead_time, delay)

cncastillo · 2024-02-02T21:49:38Z

KomaMRIBase/src/sequences/RF/design.jl

+    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


Not required.

cncastillo · 2024-02-02T21:51:42Z

KomaMRIBase/src/sequences/ADC/design.jl

Separate functions in different files. All of them, not only the ones in ADC.

cncastillo · 2024-02-02T21:52:26Z

KomaMRIBase/src/sequences/RF/design.jl

+    if dead_time > delay
+        delay = dead_time
+    end
+	rf = RF(signal, duration, freq_offset, delay)


Indentation problem.

cncastillo · 2024-02-02T21:52:45Z

KomaMRIBase/src/sequences/RF/design.jl

+    if dead_time > delay
+        delay = dead_time
+    end


delay = max(delay, dead_time)

cncastillo · 2024-02-02T21:53:35Z

KomaMRIBase/src/sequences/RF/design.jl

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


we need a function to calculate the area of a Grad (get_kspace almost does this)

cncastillo · 2024-02-02T21:57:13Z

KomaMRIBase/src/sequences/RF/design.jl

+
+    N =  length(signal)
+    duration = (N-1) * Δtrf
+    t = range(0, duration; length=N)


Is this even being used? ( t = range(0, duration; length=N))

cncastillo · 2024-02-02T21:59:13Z

KomaMRIBase/src/sequences/RF/design.jl

+    if rf.delay < gz.rise + gz.delay
+        rf.delay = gz.rise + gz.delay   # these are on the grad raster already which is coarser
+    end


delay = max(...)

cncastillo · 2024-02-02T21:59:54Z

KomaMRIBase/src/sequences/RF/design.jl

+    dly = Delay(0)
+    if ring_down_time > 0
+        dly = Delay(rf.delay + rf.T + ring_down_time)     # I NEED a review
+    end


if ring_down_time is not > 0, it must be zero, therefore dly = Delay(rf.delay + rf.T + ring_down_time) is always correct

cncastillo · 2024-02-02T22:00:47Z

KomaMRIBase/src/sequences/RF/design.jl

+    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π


Use get_flip_angle function, I think this assumes that signal is in a particular unit. Is missing the multiplication by gamma.

cncastillo · 2024-03-19T13:01:41Z

KomaMRIBase/ext/KomaMRIBaseUnitfulExt.jl

+#Angle{T} = Union{Quantity{T, NoDims, typeof(u"rad")}, Quantity{T, NoDims, typeof(u"°")}} where T
+
+function KomaMRIBase.block_pulse(x::Real)
+    return "Hola"
+end
+
+#function KomaMRIBase.block_pulse(duration::Unitful.Time;
+#    phase_offset::Angle=0u"°", freq_offset::Frequency=0u"Hz", delay=0u"s", sys=Scanner())
+#    duration, phase_offset, freq_offset, delay = upreferred.((duration, phase_offset, freq_offset, delay))
+#    return PulseDesigner.block_pulse(FlipAngle(π/2), Duration(duration); phase_offset, freq_offset, delay, sys)
+#end
+
+
+#function PulseDesigner.block_pulse(flip_angle::Angle, duration::Unitful.Time;
+#    phase_offset::Angle=0u"°", freq_offset::Frequency=0u"Hz", delay=0u"s", sys=Scanner())
+#    flip_angle, duration, phase_offset, freq_offset, delay = upreferred.((flip_angle, duration, phase_offset, freq_offset, delay))
+#    return PulseDesigner.block_pulse(FlipAngle(flip_angle), Duration(duration); phase_offset, freq_offset, delay, sys)
+#end
+#
+#function PulseDesigner.block_pulse(flip_angle::Angle, bandwidth::Unitful.Frequency;
+#    phase_offset::Angle=0u"°", freq_offset::Frequency=0u"Hz", delay=0u"s", sys=Scanner())
+#    flip_angle, bandwidth, phase_offset, freq_offset, delay = upreferred.((flip_angle, bandwidth, phase_offset, freq_offset, delay))
+#    return PulseDesigner.block_pulse(FlipAngle(flip_angle), Bandwidth(bandwidth); phase_offset, freq_offset, delay, sys)
+#end
+#
+#function PulseDesigner.block_pulse(flip_angle::Angle, bandwidth::Unitful.Frequency, time_bw_product::DimensionlessQuantity;
+#    phase_offset::Angle=0u"°", freq_offset::Frequency=0u"Hz", delay=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))
+#    return PulseDesigner.block_pulse(FlipAngle(flip_angle), Bandwidth(bandwidth), TimeBwProduct(time_bw_product); phase_offset, freq_offset, delay, sys)
+#end


cncastillo · 2024-03-19T13:02:17Z

KomaMRIBase/src/KomaMRIBase.jl

+struct Duration val end
+struct Bandwidth val end
+struct TimeBwProduct val end
+struct Angle end


Suggested change

struct Angle end

struct Angle val end

cncastillo · 2024-03-19T13:04:11Z

KomaMRIBase/src/sequences/RF/block_pulse.jl

+    dead_time = sys.RF_dead_time_T
+    ring_down_time = sys.RF_ring_down_T


why define these?

cncastillo · 2024-03-19T13:06:52Z

KomaMRIBase/src/sequences/RF/block_pulse.jl

+    block_duration = delay + duration + ring_down_time
+
+    rf = RF(amplitude, duration, freq_offset, delay)
+    return Sequence([Grad(0, 0);;], [rf;;], [ADC(0, 0)], block_duration)


we need a way to combine events so this is easier, seq += [rf, block_duration] where rf::RF and block_duration::Delay, issue #4.

Add simple high-level functions for designing sequences

133a7b1

cncastillo reviewed Feb 2, 2024

View reviewed changes

beorostica added 2 commits March 7, 2024 15:08

Merge branch 'master' into pulsedesigner

882eab9

Temporal changes for PulseDesigner ...

63fefd7

Create block_pulse.jl file

5898fdd

cncastillo mentioned this pull request Mar 13, 2024

Motion Phantom #184

Merged

14 tasks

beorostica added 2 commits March 18, 2024 10:28

Merge branch 'master' into pulsedesigner

3e4a339

Attempt for multidispatch and extensions for block_pulse function

844c3a5

cncastillo reviewed Mar 19, 2024

View reviewed changes

beorostica added 4 commits March 19, 2024 12:36

Remove Unitful extension for PulseDesigner

d1fcb58

Add sinc_pulse function without slice gradient

2d0749e

Add Unitful extensions for PulseDesigne using the block_pulse function

33f32fd

Polish block_pulse PulseDesigner function

94220ec

cncastillo mentioned this pull request Apr 1, 2024

Incorrect gradient interpretation when waveforms do not end in zero #321

Merged

beorostica mentioned this pull request Apr 3, 2024

Create types for MRISequenceEvents to enable multidispatch #362

Open

cncastillo mentioned this pull request Apr 5, 2024

Change all struct parameter types from <:Real to <:Number #363

Open

		dead_time = sys.RF_dead_time_T
		ring_down_time = sys.RF_ring_down_T

Conversation

beorostica commented Feb 2, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Mar 19, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cncastillo Feb 2, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

cncastillo Feb 2, 2024 •

edited

Loading

cncastillo Mar 19, 2024 •

edited

Loading

cncastillo Feb 2, 2024 •

edited

Loading

cncastillo Feb 2, 2024 •

edited

Loading

codecov bot commented Feb 2, 2024 •

edited

Loading

cncastillo Feb 2, 2024 •

edited

Loading

cncastillo Feb 2, 2024 •

edited

Loading

cncastillo Feb 2, 2024 •

edited

Loading