Add magnetisation transfer sequence

src/MRIBuilder.jl

@@ -25,8 +25,8 @@ export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra, Default_Sc
 import .Variables: variables, make_generic, @defvar, get_pulse, get_readout, get_pathway, get_gradient, add_cost_function!, apply_simple_constraint!, set_simple_constraints!
 export variables, make_generic, @defvar, get_pulse, get_readout, get_pathway, get_gradient, add_cost_function!, apply_simple_constraint!, set_simple_constraints!
 
-import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
-export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
+import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times, BinomialPulse
+export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times, BinomialPulse
 
 import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AbstractAlternativeBlocks, get_alternatives_name, get_alternatives_options, AlternativeBlocks, match_blocks!, get_index_single_TR, iter_blocks, iter_instant_gradients, iter_instant_pulses
 export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AbstractAlternativeBlocks, get_alternatives_name, get_alternatives_options, AlternativeBlocks, match_blocks!, get_index_single_TR, iter_blocks, iter_instant_gradients, iter_instant_pulses
@@ -34,14 +34,14 @@ export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events
 import .Pathways: Pathway
 export Pathway
 
-import .Parts: dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
-export dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
+import .Parts: dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size, saturation_pulse_binomial
+export dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size, saturation_pulse_binomial
 
 import .PostHoc: adjust, merge_sequences
 export adjust, merge_sequences
 
-import .Sequences: GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI
-export GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI
+import .Sequences: GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI, MagnetisationTransfer
+export GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI, MagnetisationTransfer
 
 import .SequenceIO: read_sequence, write_sequence
 export read_sequence, write_sequence

src/components/components.jl

@@ -8,7 +8,7 @@ include("readouts/readouts.jl")
 import .AbstractTypes: BaseComponent, GradientWaveform, EventComponent, RFPulseComponent, ReadoutComponent, split_timestep, edge_times
 import .GradientWaveforms: ConstantGradient, ChangingGradient, NoGradient, split_gradient
 import .InstantGradients: InstantGradient, InstantGradient3D, InstantGradient1D
-import .Pulses: GenericPulse, InstantPulse, SincPulse, ConstantPulse, CompositePulse
+import .Pulses: GenericPulse, InstantPulse, SincPulse, ConstantPulse, CompositePulse, BinomialPulse
 import .Readouts: ADC, SingleReadout
 
 end
\ No newline at end of file

src/components/pulses/composite_pulses.jl

@@ -5,7 +5,7 @@ import ....Variables: VariableType, set_simple_constraints!, make_generic, get_f
 import ..GenericPulses: GenericPulse
 
 """
-    CompositePulse(; variables...)
+    CompositePulse(; base_pulse, nweights, variables...)
 
 A composite RF pulse formed by repeating a base RF pulse.
 
@@ -18,8 +18,8 @@ A composite RF pulse formed by repeating a base RF pulse.
 - `interpulse_delay`: Time between the center of the RF pulses. If not otherwise constrained, it will be minimised.
 - `scale_amplitude`: How strongly one should scale the amplitude versus the duration to achieve the desired weights. If set to 1 only the RF pulse amplitude will be scaled. If set to 0 only the RF pulse duration will be scaled.
 """
-struct CompositePulse <: RFPulseComponent
-    pulses :: Vector{RFPulseComponent}
+struct CompositePulse{T<:RFPulseComponent} <: RFPulseComponent
+    pulses :: Vector{T}
     interpulse_delay :: VariableType
     scale_amplitude :: VariableType
 end
@@ -42,7 +42,7 @@ function CompositePulse(; base_pulse::RFPulseComponent, nweights=nothing, weight
         )
         for w in weights
     ]
-    res = CompositePulse(
+    res = CompositePulse{typeof(base_pulse)}(
         pulses,
         get_free_variable(interpulse_delay),
         scale_amplitude
@@ -60,7 +60,7 @@ end
 @defvar duration(pulse::CompositePulse) = (
         0.5 * variables.duration(pulse.pulses[1]) + 
         0.5 * variables.duration(pulse.pulses[end]) + 
-        interpulse_delay(pulse) * (length(pulse) - 1)
+        variables.interpulse_delay(pulse) * (length(pulse) - 1)
     )
 
 @defvar begin
@@ -142,4 +142,18 @@ function adjust_internal(comp::CompositePulse; stretch=1., kwargs...)
     )
 end
 
+
+"""
+    BinomialPulse(; base_pulse, npulses, flip=true, variables...)
+
+Creates a [`CompositePulse`](@ref) of `npulses` repeats of `base_pulse`, where the `weights` are set by the biomial distribution.
+
+If `flip` is true (default) every other pulse will be flipped (so that the total excitation is canceled).
+The `variables` are defined in [`CompositePulse`](@ref).
+"""
+function BinomialPulse(; npulses, flip=true, kwargs...)
+    weights = [binomial(npulses-1, index) * (flip && isone(index % 2) ? -1 : 1) for index in 0:npulses-1]
+    return CompositePulse(; kwargs..., weights=weights)
+end
+
 end
\ No newline at end of file

src/components/pulses/constant_pulses.jl

@@ -35,6 +35,7 @@ function ConstantPulse(; amplitude=nothing, duration=nothing, phase=nothing, fre
     apply_simple_constraint!(res.amplitude, :>=, 0)
     set_simple_constraints!(res, kwargs)
     add_cost_function!(res.frequency^2 + res.phase^2)
+    add_cost_function!((res.flip_angle - 90)^2)
     return res
 end
 

src/components/pulses/instant_pulses.jl

@@ -29,6 +29,7 @@ function InstantPulse(; flip_angle=nothing, phase=nothing, group=nothing)
     )
     apply_simple_constraint!(res.flip_angle, :>=, 0)
     add_cost_function!(res.phase^2)
+    add_cost_function!((res.flip_angle - 90)^2)
     return res
 end
 

src/components/pulses/pulses.jl

@@ -10,6 +10,6 @@ import .GenericPulses: GenericPulse
 import .InstantPulses: InstantPulse
 import .ConstantPulses: ConstantPulse
 import .SincPulses: SincPulse
-import .CompositePulses: CompositePulse
+import .CompositePulses: CompositePulse, BinomialPulse
 
 end
\ No newline at end of file

src/components/pulses/sinc_pulses.jl

@@ -56,6 +56,7 @@ function SincPulse(;
     end
     set_simple_constraints!(res, kwargs)
     add_cost_function!(res.frequency^2 + res.phase^2)
+    add_cost_function!((res.flip_angle - 90)^2)
     return res
 end
 

src/parts/helper_functions.jl

@@ -5,10 +5,11 @@ import ..Trapezoids: Trapezoid, opposite_kspace_lines, SliceSelect
 import ..SpoiltSliceSelects: SpoiltSliceSelect
 import ..SliceSelectRephases: SliceSelectRephase
 import ..EPIReadouts: EPIReadout
+import ..Repeats: Repeat
 import ...BuildSequences: build_sequence, global_scanner
 import ...Containers: Sequence
-import ...Components: SincPulse, ConstantPulse, InstantPulse, SingleReadout, InstantGradient
-import ...Variables: variables, apply_simple_constraint!
+import ...Components: SincPulse, ConstantPulse, InstantPulse, SingleReadout, InstantGradient, BinomialPulse
+import ...Variables: variables, apply_simple_constraint!, set_simple_constraints!
 
 
 function _get_pulse(shape, flip_angle, phase, frequency, Nzeros, group, bandwidth, duration)
@@ -306,4 +307,40 @@ function gradient_spoiler(; optimise=false, scanner=nothing, orientation=[0, 0,
     end
 end
 
+
+
+
+"""
+    saturation_pulse(; pulse=(), binomial_order=1, variables...)
+
+Creates a saturation pulse consisting of repeating instances of `pulse_type`.
+
+## Parameters
+- `pulse`: Parameters/variables used to define the basic RF pulse used in the saturation pulse. These can be:
+    - `type`: [`ConstantPulse`](@ref) (default), [`SincPulse`](@ref), [`InstantPulse`](@ref), etc.
+    - Variables describing that type (e.g., `flip_angle`, `duration`, `frequency`)
+- `binomial_order`: How many repeats of `pulse_type` there should be in a single binomial pulse.
+
+## Variables
+- `nrepeat`: number of repeats of the binomial pulse. This should typically be set by the user.
+- `interblock_delay`: time between the start of each pulse. Should be short compared with the T1, but long enough to limit SAR.
+- `duration`: total duration of the saturation pulse.
+"""
+function saturation_pulse_binomial(; pulse=(), binomial_order=1, nrepeat, variables...)
+
+    function base_pulse(; type=ConstantPulse, variables...)
+        return type(; variables...)
+    end
+    base = base_pulse(; pulse...)
+
+    if binomial_order > 1
+        composite = BinomialPulse(; base_pulse=base, npulses=binomial_order)
+    else
+        composite = base
+    end
+    res = Repeat(composite, nrepeat=nrepeat)
+    set_simple_constraints!(res, variables)
+    return res
+end
+
 end
\ No newline at end of file

src/parts/parts.jl

@@ -3,13 +3,15 @@ include("trapezoids.jl")
 include("slice_select_rephases.jl")
 include("spoilt_slice_selects.jl")
 include("epi_readouts.jl")
+include("repeats.jl")
 include("helper_functions.jl")
 
 import .Trapezoids: Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines
 import .SpoiltSliceSelects: SpoiltSliceSelect
 import .SliceSelectRephases: SliceSelectRephase
 import .EPIReadouts: EPIReadout
-import .HelperFunctions: excitation_pulse, refocus_pulse, readout_event, dwi_gradients, interpret_image_size, gradient_spoiler
+import .Repeats: Repeat
+import .HelperFunctions: excitation_pulse, refocus_pulse, readout_event, dwi_gradients, interpret_image_size, gradient_spoiler, saturation_pulse_binomial
 
 
 end
\ No newline at end of file

src/parts/repeats.jl

+module Repeats
+import ...Containers: BaseSequence, get_index_single_TR, Wait
+import ...Variables: variables, VariableType, AbstractBlock, get_free_variable, set_simple_constraints!, apply_simple_constraint!, @defvar, get_pulse, get_gradient, get_readout
+
+"""
+    Repeat(block, variables...)
+
+# Parameters
+- `block`: RF pulse/gradient/ADC/sequence part that needs to be repeated.
+
+# Variables
+- `nrepeat`: number of times the block is repeated.
+- `wait_time`: time between end of one block and start of the next repeat in ms.
+- `interblock_delay`: time between start of one block and start of the next repeat in ms.
+- `frequency`: temporal frequency of the block repeats (i.e., `1/interblock_delay`) in kHz.
+- `duration`: total duration of this sequence part all repeats in ms.
+"""
+struct Repeat{T<:AbstractBlock} <: BaseSequence{1}
+    block :: T
+    wait_time :: VariableType
+    nrepeat :: VariableType
+end
+
+function Repeat(block::AbstractBlock; wait_time=nothing, nrepeat=nothing, variables...)
+    res = Repeat{typeof(block)}(block, get_free_variable(wait_time), get_free_variable(nrepeat, integer=true))
+    set_simple_constraints!(res, variables)
+    apply_simple_constraint!(res.wait_time, :>=, 0.)
+    apply_simple_constraint!(res.nrepeat, :>=, 1)
+    return res
+end
+
+function get_index_single_TR(r::Repeat, i::Integer)
+    if iszero(i % 2)
+        return r.block
+    else
+        return Wait(r.wait_time)
+    end
+end
+
+Base.length(r::Repeat) = r.nrepeat * 2 - 1
+
+@defvar begin
+    nrepeat(r::Repeat) = r.nrepeat
+    wait_time(r::Repeat) = r.wait_time
+    duration(r::Repeat) = r.interblock_delay * r.nrepeat - r.wait_time
+    interblock_delay(r::Repeat) = r.block.duration + r.wait_time
+    frequency(r::Repeat) = 1/(r.interblock_delay)
+end
+
+get_pulse(r::Repeat) = get_pulse(r.block)
+get_gradient(r::Repeat) = get_gradient(r.block)
+get_readout(r::Repeat) = get_readout(r.block)
+
+end
\ No newline at end of file

src/sequences/magnetisation_transfers.jl

+module MagnetisationTransfers
+import ...Containers: Sequence
+import ...Parts: excitation_pulse, readout_event, interpret_image_size, saturation_pulse_binomial
+import ...Variables: get_pulse, get_readout
+import ...BuildSequences: build_sequence
+import ...Scanners: Default_Scanner
+
+const MagnetisationTransfer = Sequence{:MagnetisationTransfer}
+
+
+"""
+    MagnetisationTransfer(; saturation=(), excitation=(), readout=(), optim=())
+"""
+function MagnetisationTransfer(; saturation=(), excitation=(), readout=(), optim=(), resolution=nothing, fov=nothing, voxel_size=nothing, slice_thickness=nothing, scanner=Default_Scanner, vars...)
+    build_sequence(scanner; optim...) do
+        (slice_thickness, _, extra_readout_params) = interpret_image_size(fov, resolution, voxel_size, slice_thickness)
+        parts = Any[
+            :saturation => saturation_pulse_binomial(; saturation...),
+            nothing,
+            :excitation => excitation_pulse(; slice_thickness=slice_thickness, excitation...),
+            nothing,
+            :readout => readout_event(; extra_readout_params..., readout...),
+            nothing,
+        ]
+        seq = Sequence(parts; name=:MagnetistionTransfer, vars...)
+        return seq
+    end
+end
+
+
+get_pulse(mt::MagnetisationTransfer) = (saturation=mt[:excitation].base_pulse, refocus=ge[:refocus])
+get_readout(mt::MagnetisationTransfer) = ge[:readout]
+
+
+end
\ No newline at end of file

src/sequences/sequences.jl

@@ -2,9 +2,11 @@ module Sequences
 include("gradient_echoes.jl")
 include("spin_echoes.jl")
 include("diffusion_spin_echoes.jl")
+include("magnetisation_transfers.jl")
 
 import .GradientEchoes: GradientEcho
 import .SpinEchoes: SpinEcho
 import .DiffusionSpinEchoes: DiffusionSpinEcho, DW_SE, DWI
+import .MagnetisationTransfers: MagnetisationTransfer
 
 end
\ No newline at end of file