Allow for multiple alternative blocks in part of sequence

src/MRIBuilder.jl

@@ -13,6 +13,7 @@ include("pulses/pulses.jl")
 include("readouts/readouts.jl")
 include("overlapping/overlapping.jl")
 include("sequences.jl")
+include("alternatives.jl")
 include("pathways.jl")
 include("helper_functions.jl")
 
@@ -46,6 +47,9 @@ export TrapezoidGradient, SpoiltSliceSelect, interruptions, waveform, SingleLine
 import .Sequences: Sequence
 export Sequence
 
+import .Alternatives: AlternativeBlocks
+export AlternativeBlocks
+
 import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat
 export Pathway, duration_transverse, duration_dephase, bval, bmat
 

src/alternatives.jl

+module Alternatives
+import JuMP: @constraint
+import ..BuildingBlocks: BuildingBlock, match_blocks!
+import ..BuildSequences: global_model
+import ..Variables: duration
+
+"""
+    AlternativeBlocks(name, blocks)
+
+Represents a part of the sequence where there are multiple possible alternatives.
+
+Variables can be matched across these alternatives using [`match_blocks!`](@ref).
+
+The `name` is a symbol that is used to identify this `AlternativeBlocks` in the broader sequence.
+"""
+struct AlternativeBlocks <: BuildingBlock
+    name :: Symbol
+    options :: Vector{<:BuildingBlock}
+end
+
+Base.getindex(alt::AlternativeBlocks, index::Int) = alt.options[index]
+
+duration(alt::AlternativeBlocks) = maximum(duration.(alt.options))
+
+"""
+    match_blocks!(alternatives, function)
+
+Matches the outcome of given `function` on each of the building blocks in [`AlternativeBlocks`](@ref).
+
+For example, `match_blocks!(alternatives, duration)` will ensure that all the alternative building blocks have the same duration.
+"""
+function match_blocks!(alternatives::AlternativeBlocks, func)
+    baseline = func(alternatives[1])
+    for other_block in alternatives.options[2:end]
+        if baseline isa AbstractVector
+            @constraint global_model() baseline == func(other_block)
+        else
+            @constraint global_model() baseline .== func(other_block)
+        end
+    end
+end
+
+
+
+end
\ No newline at end of file

src/helper_functions.jl

@@ -5,7 +5,8 @@ import ..BuildSequences: global_model, build_sequence
 import ..Sequences: Sequence
 import ..Pulses: SincPulse, ConstantPulse, InstantRFPulseBlock
 import ..Overlapping: TrapezoidGradient, SpoiltSliceSelect, opposite_kspace_lines
-import ..Variables: qvec
+import ..Variables: qvec, flat_time, rise_time
+import ..Alternatives: AlternativeBlocks, match_blocks!
 
 
 function _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
@@ -198,10 +199,18 @@ Helper function used to build the readout for any Cartesian readout, i.e.:
 function cartesian_readout(start_lines, readout_lines, fov, resolution_x; scanner=nothing, optimise=false, kwargs...)
     @assert iszero(readout_lines[1])
     build_sequence(scanner; optimise=optimise) do
-        max_ky = maximum(abs.(start_lines)) * 1e-3 / fov[2]
-        to_scale_ky = length(start_lines) == 1 ? nothing : (nothing, :phase_encode, nothing)
-        prepare_kspace = TrapezoidGradient(scale=to_scale_ky, rotate=:FOV, duration=:min)
         (pos_line, neg_line) = opposite_kspace_lines(; rotate=:FOV, fov=fov[1], resolution=resolution_x, kwargs...)
+        ky = @. start_lines * 1e-3 / fov[2]
+        if length(start_lines) == 1
+            prepare_kspace = TrapezoidGradient(rotate=:FOV, duration=:min, qvec=[-qvec(pos_line, nothing, 1)[1], ky[1], 0.])
+        else
+            prepare_kspace = AlternativeBlocks(
+                :readout_segment,
+                [TrapezoidGradient(rotate=:FOV, duration=:min, qvec=[-qvec(pos_line, nothing, 1)[1], ky_prep, 0.]) for ky_prep in ky]
+            )
+            match_blocks!(prepare_kspace, flat_time)
+            match_blocks!(prepare_kspace, rise_time)
+        end
 
         steps = (readout_lines[2:end] - readout_lines[1:end-1])
         blips = Dict(
@@ -218,9 +227,6 @@ function cartesian_readout(start_lines, readout_lines, fov, resolution_x; scanne
             append!(result, [blips[s], next_line])
             next_line = next_line == pos_line ? neg_line : pos_line
         end
-        @constraint global_model() qvec(prepare_kspace)[1] == -qvec(pos_line, nothing, 1)[1]
-        @constraint global_model() qvec(prepare_kspace)[2] == max_ky
-        @constraint global_model() qvec(prepare_kspace)[3] == 0.
         return Sequence(result...; TR=Inf)
     end
 end