Add gradient echo sequence

src/MRIBuilder.jl

@@ -10,6 +10,7 @@ include("components/components.jl")
 include("containers/containers.jl")
 include("pathways.jl")
 include("parts/parts.jl")
+include("sequences/sequences.jl")
 include("printing.jl")
 
 import .BuildSequences: build_sequence, global_model, global_scanner, fixed
@@ -18,8 +19,8 @@ export build_sequence, global_model, global_scanner, fixed
 import .Scanners: Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 
-import .Variables: variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inverse_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, TR, Δ, get_gradient, get_pulse, get_readout
-export variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inversne_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, TR, Δ, get_gradient, get_pulse, get_readout
+import .Variables: variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inverse_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time
+export variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inversne_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time
 
 import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
 export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
@@ -33,6 +34,9 @@ export Pathway, duration_transverse, duration_dephase, bval, bmat, get_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 .Sequences: GradientEcho
+export GradientEcho
+
 import JuMP: @constraint, @objective, objective_function, value, Model
 export @constraint, @objective, objective_function, value, Model
 

src/containers/sequences.jl → src/containers/base_sequences.jl

 """
 Defines [`BaseSequence`](@ref) and [`Sequence`](@ref)
 """
-module Sequences
+module BaseSequences
 import StaticArrays: SVector
 import JuMP: @constraint
-import ...Variables: get_free_variable, TR, VariableType, duration, variables, VariableNotAvailable, Variables, set_simple_constraints!
+import ...Variables: get_free_variable, repetition_time, VariableType, duration, variables, VariableNotAvailable, Variables, set_simple_constraints!, TR
 import ...BuildSequences: global_model
 import ...Components: EventComponent
 import ..Abstract: ContainerBlock, start_time
@@ -22,8 +22,6 @@ Main interface:
 
 Sub-types need to implement:
 - `get_index_single_TR`: return the index assuming it is between 1 and N
-- [`nrepeat`](@ref): how often the sequence should repeat (if not implemented, this will be 1).
-- [`TR`](@ref): time scale on which to repeat (if not implemented, this will be the sum of the individual block durations).
 """
 abstract type BaseSequence{N} <: ContainerBlock end
 
@@ -49,7 +47,7 @@ function start_time(bs::BaseSequence{N}, index::Integer) where {N}
     if iszero(nTR)
         return base_time
     else
-        return nTR * TR(bs) + base_time
+        return nTR * repetition_time(bs) + base_time
     end
 end
 
@@ -79,28 +77,28 @@ How often sequence should be repeated.
 """
 nrepeat(bs::BaseSequence) = 1
 
-"""
-    TR(sequence)
-
-On what timescale the sequence should repeat itself in ms (will have no effect if [`nrepeat`](@ref) is one).
-
-By default this is set to the total duration of the sequence.
-"""
-TR(bs::BaseSequence) = duration(bs)
+repetition_time(bs::BaseSequence) = duration(bs)
 
 
 duration(bs::BaseSequence{0}) = 0.
 duration(bs::BaseSequence) = sum(duration.(bs); init=0.)
 
 """
-    Sequence(blocks; TR=:min, nrepeat=0)
-    Sequence(blocks...; TR=:min, nrepeat=0)
+    Sequence(blocks; name=:Sequence, variables...)
+    Sequence(blocks...; name=:Sequence, variables...)
 
 Defines an MRI sequence from a vector of building blocks.
 
 ## Arguments
-- [`nrepeat`](@ref): how often the sequence will repeat itself (keep at default of 0 to repeat indefinetely).
-- [`TR`](@ref): how long between repeats in ms (defaults to the duration of the sequence). Can be set to `nothing` to be a free variable.
+- `blocks`: The actual building blocks that will be played in sequence. All the building blocks must be of type [`ContainerBlock`](@ref), which means that they cannot only contain actual [`BaseBuildingBlock`](@ref) objects, but also other [`BaseSequence`](@ref) objects.
+    Objects of a different type are converted into a [`ContainerBlock`](@ref) internally:
+    - numbers/`nothing`/`:min`/`:max` : replaced with a [`Wait`](@ref) block with the appropriate constraint/objective added to its [`duration`](@ref).
+    - RF pulse or readout: will be embedded within a [`BuildingBlock`](@ref) of the appropriate length
+
+## Variables
+- [`repetition_time`](@ref)/[`TR`](@ref): how long between repeats in ms. This is always set to the total length of the sequence. If you want to add some down-time between repeats, you can simply add a [`Wait`](@ref) block of the appropriate length at the end of the sequence.
+
+Specific named sequences might define additional variables.
 """
 struct Sequence{S, N} <: BaseSequence{N}
     blocks :: SVector{N, Pair{<:Union{Symbol, Nothing}, <:ContainerBlock}}

src/containers/containers.jl

 module Containers
 include("abstract.jl")
 include("building_blocks.jl")
-include("sequences.jl")
+include("base_sequences.jl")
 include("alternatives.jl")
 
 import .Abstract: ContainerBlock, start_time, end_time, effective_time
 import .BuildingBlocks: BaseBuildingBlock, BuildingBlock, Wait, waveform, waveform_sequence, events
-import .Sequences: BaseSequence, Sequence, nrepeat, get_index_single_TR
+import .BaseSequences: BaseSequence, Sequence, nrepeat, get_index_single_TR
 import .Alternatives: AlternativeBlocks, match_blocks!
 
 end
\ No newline at end of file

src/scanners.jl

@@ -63,6 +63,13 @@ By setting `units` to :Tesla, the slew rate can be returned in T/m/s instead.
 slew_rate(scanner::Scanner, units=:kHz) = units == :kHz ? scanner.slew_rate : scanner.slew_rate / (gyromagnetic_ratio * 1e-9)
 
 
+"""
+A default 1.5T scanner.
+
+Matches the one used in `pulseq` (https://github.com/pulseq/pulseq/blob/master/matlab/%2Bmr/opts.m).
+"""
+Default_Scanner = Scanner(B0=1.5, gradient=40, slew_rate=170, units=:Tesla)
+
 """
 Siemens MAGNETOM 3T Prisma MRI scanner (https://www.siemens-healthineers.com/en-uk/magnetic-resonance-imaging/3t-mri-scanner/magnetom-prisma).
 """

src/sequences/gradient_echos.jl

+module GradientEchos
+import ...Containers: Sequence
+import ...Parts: excitation_pulse, readout_event, interpret_image_size, Trapezoid, gradient_spoiler
+import ...Variables: get_pulse, get_readout, echo_time, duration_transverse
+import ...Pathways: Pathway, get_pathway
+import ...BuildSequences: build_sequence
+import ...Scanners: Default_Scanner
+
+const GradientEcho = Sequence{:GradientEcho}
+
+"""
+    GradientEcho(; echo_time, excitation=(), readout=(), resolution/fov/voxel_size/slice_thickness)
+
+Defines a gradient echo sequence with a single readout event.
+
+By default, an instant excitation pulse and readout event are used.
+If image parameters are provided, this will switch to a sinc pulse and EPI readout.
+
+## Parameters
+- [`excitation`](@ref): properties of the excitation pulse as described in [`excitation_pulse`](@ref).
+- [`readout`](@ref): properties of the readout as described in [`readout_event`](@ref).
+- [`spoiler`](@ref): if set adds a spoiler [`gradient_spoiler`](@ref) gradient after the readout (e.g., `spoiler=(spoiler_scale=1, orientation=[0, 0, 1], group=:FOV)` to add a gradient in the z-direction of the `FOV` coordinate system that fully dephases spins over 1 mm).
+- Image parameters ([`resolution`](@ref)/[`fov`](@ref)/[`voxel_size`](@ref)/[`slice_thickness`](@ref)): describe the properties of the resulting image. See [`interpret_image_size`](@ref) for details.
+
+## Variables
+- [`TE`](@ref)/[`echo_time`](@ref): echo time between excitation pulse and readout in ms (required).
+- [`TR`](@ref)/[`repetition_time`](@ref)/[`duration`](@ref): total duration of the sequence from start of excitation pulse to end of readout or spoiler in ms.
+"""
+function GradientEcho(; excitation=(), readout=(), spoiler=nothing, resolution=nothing, fov=nothing, voxel_size=nothing, slice_thickness=nothing, scanner=Default_Scanner, variables...)
+    build_sequence(scanner) do
+        (slice_thickness, _, extra_readout_params) = interpret_image_size(fov, resolution, voxel_size, slice_thickness)
+        parts = Any[
+            :excitation => excitation_pulse(; slice_thickness=slice_thickness, excitation...),
+            nothing,
+            :readout => readout_event(; extra_readout_params..., readout...),
+        ]
+        if !isnothing(spoiler)
+            push!(parts, gradient_spoiler(; spoiler...))
+        end
+        return Sequence(parts; name=:GradientEcho, variables...)
+    end
+end
+
+
+get_pulse(ge::GradientEcho) = ge.excitation
+get_readout(ge::GradientEcho) = ge.readout
+get_pathway(ge::GradientEcho) = Pathway(ge, [90], 1)
+echo_time(ge::GradientEcho) = duration_transverse(ge)
+
+
+
+end
\ No newline at end of file

src/sequences/sequences.jl

+module Sequences
+include("gradient_echos.jl")
+
+import .GradientEchos: GradientEcho
+
+end
\ No newline at end of file

src/variables.jl

@@ -34,7 +34,10 @@ all_variables_symbols = [
         :duration => "duration of the building block in ms.",
     ],
     :sequence => [
-        :TR => "Time on which an MRI sequence repeats itself in ms.",
+        :TR => "Time on which an MRI sequence repeats itself in ms. Defaults to the result of [`repetition_time`](@ref).",
+        :repetition_time => "Time on which an MRI sequence repeats itself in ms.",
+        :TE => "Echo time of the sequence in ms. Defaults to the result of [`echo_time`](@ref).",
+        :echo_time => "Echo time of the sequence in ms.",
         :Δ => "Diffusion time in ms (i.e., time between start of the diffusion-weighted gradients).",
         :qvec => "Net dephasing due to gradients in rad/um.",
         :area_under_curve => "Net dephasing due to gradients in rad/um (same as [`qvec`](@ref)).",
@@ -98,6 +101,9 @@ for (block_symbol, all_functions) in all_variables_symbols
 end
 
 
+TE(ab::AbstractBlock) = echo_time(ab)
+TR(ab::AbstractBlock) = repetition_time(ab)
+
 """
 Dictionary with alternative versions of specific function.
     
@@ -228,7 +234,7 @@ end
 
 for (target_name, all_vars) in all_variables_symbols
     for (variable_func, _) in all_vars
-        if variable_func == :qval3
+        if variable_func in [:qval3, :TR, :TE]
             continue
         end
         get_func = Symbol("get_" * string(target_name))