Add spoilt slice select gradient waveform

src/MRIBuilder.jl

@@ -40,8 +40,8 @@ export InstantGradientBlock
 import .Readouts: InstantReadout
 export InstantReadout
 
-import .Overlapping: TrapezoidGradient
-export TrapezoidGradient
+import .Overlapping: TrapezoidGradient, SpoiltSliceSelect, interruptions, waveform
+export TrapezoidGradient, SpoiltSliceSelect, interruptions, waveform
 
 import .Sequences: Sequence
 export Sequence

src/build_sequences.jl

@@ -55,7 +55,7 @@ end
 function build_sequence(f::Function, scanner::Scanner)
     ipopt_opt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3)
     juniper_opt = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt_opt)
-    model = Model(juniper_opt)
+    model = Model(ipopt_opt)
     build_sequence(f, scanner, model)
 end
 

src/building_blocks.jl

@@ -145,6 +145,9 @@ function _robust_value(possible_vector::AbstractVector)
     return result
 end
 
+_robust_value(possible_tuple::Tuple) = _robust_value([possible_tuple...])
+
+
 function Base.show(io::IO, printer::BuildingBlockPrinter)
     block = printer.bb
     print(io, string(typeof(block)), "(")

src/helper_functions.jl

@@ -3,7 +3,7 @@ import JuMP: @constraint
 import ..BuildSequences: get_global_model
 import ..Sequences: Sequence
 import ..Pulses: SincPulse, ConstantPulse, InstantRFPulseBlock
-import ..Overlapping: TrapezoidGradient
+import ..Overlapping: TrapezoidGradient, SpoiltSliceSelect
 import ..Variables: qvec
 
 
@@ -40,23 +40,23 @@ For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_an
 - `scale`: name of the parameter with which the RF pulse amplitude can be modulated after sequence optimisation (default: `:transmit_B1`).
 
 ### Slice selection
-- `min_slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
+- `slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
 - `rephase`: set to false to disable the spin rephasing after the RF pulse.
 - `rotate_grad`: name of the parameter with which the slice selection gradient will be rotated after sequence optimisation (default: `:FOV`).
 """
-function excitation_pulse(; flip_angle=90, phase=0., frequency=0., shape=:sinc, min_slice_thickness=Inf, rephase=true, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing)
+function excitation_pulse(; flip_angle=90, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, rephase=true, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing)
     pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
     if pulse isa InstantRFPulseBlock
-        if !isinf(min_slice_thickness)
-            error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `min_slice_thickness`.")
+        if !isinf(slice_thickness)
+            error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
         end
         return pulse
     end
-    if isinf(min_slice_thickness)
+    if isinf(slice_thickness)
         return pulse
     end
 
-    grad = TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, min_slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+    grad = TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
     if !rephase
         return grad
     end
@@ -77,8 +77,8 @@ end
 Create an excitation RF pulse.
 
 By default there is no slice-selective gradient. 
-To enable slice selection `min_slice_thickness` has to be set to a number or to :min.
-If `min_slice_thickness` is not set or is set to `:min`, then either `bandwidth` or `duration` should be set, otherwise the optimisation might be unconstrained (ignore this for `shape=:instant`).
+To enable slice selection `slice_thickness` has to be set to a number or to :min.
+If `slice_thickness` is not set or is set to `:min`, then either `bandwidth` or `duration` should be set, otherwise the optimisation might be unconstrained (ignore this for `shape=:instant`).
 
 ## Parameters
 For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_angle`, `phase`, and `scale` will be used. All other parameters are ignored.
@@ -92,31 +92,32 @@ For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_an
 - `Nzero`: sets the number of zero crossings for a [`SincPulse`](@ref) (default: 3). Can be set to a tuple of two numbers to set a different number of zero crossings before and after the pulse maximum.
 - `scale`: name of the parameter with which the RF pulse amplitude can be modulated after sequence optimisation (default: `:transmit_B1`).
 
-### Slice selection
-- `min_slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
-- `rotate_grad`: name of the parameter with which the slice selection gradient will be rotated after sequence optimisation (default: `:FOV`).
-
-### Spoilers
+### Slice selection and spoilers
+- `slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
 - `spoiler`: set to the spatial scale on which the spins should be dephased in mm. For rotating spoilers, this does include the contribution from the slice select gradient as well.
-- `rotate_spoiler`: set to true to rotate the spoilers with the slice select gradients. Otherwise, the spoilers will have a fixed orientation.
+- `rotate_grad`: name of the parameter with which the slice selection and spoiler gradient will be rotated after sequence optimisation (default: `:FOV`).
 """
-function refocus_pulse(; flip_angle=180, phase=0., frequency=0., shape=:sinc, min_slice_thickness=Inf, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing, spoiler=Inf, rotate_spoiler=false)
+function refocus_pulse(; flip_angle=180, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing, spoiler=Inf)
     pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
-    if pulse isa InstantRFPulseBlock && !isinf(min_slice_thickness)
-        error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `min_slice_thickness`.")
-    end
-    if !isinf(min_slice_thickness)
-        pulse = TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, min_slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+    if pulse isa InstantRFPulseBlock && !isinf(slice_thickness)
+        error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
     end
     if isinf(spoiler)
-        return pulse
-    end
-    if rotate_spoiler
-        error("spoiler rotation not implemented yet.")
+        if isinf(slice_thickness)
+            return pulse
+        else
+            return TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+        end
     else
-        grad = TrapezoidGradient(orientation=[1, 1, 1], duration=:min)
-        @constraint get_global_model() qvec(grad)[1] == 2π * 1e-3 / spoiler
-        return Sequence(grad, pulse, grad; TR=Inf)
+        orientation=isnothing(rotate_grad) ? [1, 1, 1] : [0, 0, 1]
+        if isinf(slice_thickness)
+            grad = TrapezoidGradient(orientation=orientation, duration=:min, rotate=rotate_grad)
+            @constraint get_global_model() qvec(grad)[3] == 2π * 1e-3 / spoiler
+            return Sequence(grad, pulse, grad; TR=Inf)
+        else
+            res = SpoiltSliceSelect(pulse; orientation=orientation, duration=:min, rotate=rotate_grad, slice_thickness=slice_thickness, spoiler_scale=spoiler)
+            return res
+        end
     end
 end
 

src/overlapping/overlapping.jl

@@ -2,9 +2,11 @@ module Overlapping
 include("abstract.jl")
 include("generic.jl")
 include("trapezoid_gradients.jl")
+include("spoilt_slice_selects.jl")
 
 import .Abstract: AbstractOverlapping, interruptions, waveform
 import .Generic: GenericOverlapping
 import .TrapezoidGradients: TrapezoidGradient
+import .SpoiltSliceSelects: SpoiltSliceSelect
 
 end
\ No newline at end of file

src/overlapping/spoilt_slice_selects.jl

+module SpoiltSliceSelects
+
+import LinearAlgebra: norm
+import StaticArrays: SVector
+import JuMP: @variable, @constraint, @objective, objective_function
+import ...BuildingBlocks: RFPulseBlock, set_simple_constraints!
+import ...BuildSequences: GLOBAL_SCANNER, get_global_model
+import ...Variables: VariableType, variables, duration, rise_time, flat_time, effective_time, qvec, gradient_strength, slew_rate
+import ...Gradients: ChangingGradientBlock, ConstantGradientBlock
+import ..Abstract: interruptions, waveform, AbstractOverlapping
+
+
+"""
+    SpoiltSliceSelect(pulse; parameters, variables...)
+
+Adds slice selection to the `pulse` and surrounds it with spoiler gradients.
+
+## Parameters
+- `orientation`: vector with orientation of the slice selection and the spoilers (default: [0, 0, 1])
+- `rotate`: with which user-set parameter will this gradient be rotated (e.g., :FOV). Default: no rotation.
+
+## Variables
+- `duration`: total duration of the block in ms.
+- `slice_thickness`: slice thickness in mm.
+- `spoiler_scale`: length scale on which the spoilers achieve 2π dephasing in mm.
+"""
+struct SpoiltSliceSelect <: AbstractOverlapping
+    orientation :: SVector{3, Float64}
+    rise_time1 :: VariableType
+    flat_time1 :: VariableType
+    diff_time :: VariableType
+    pulse :: RFPulseBlock
+    flat_time2 :: VariableType
+    fall_time2 :: VariableType
+    rotate :: Union{Nothing, Symbol}
+    slew_rate :: Number
+end
+
+function SpoiltSliceSelect(pulse; orientation=[0, 0, 1], rotate=nothing, spoiler_scale=nothing, kwargs...)
+    model = get_global_model()
+    res = SpoiltSliceSelect(
+        orientation ./ norm(orientation),
+        @variable(model, start=0.1),
+        @variable(model, start=0.1),
+        @variable(model, start=0.05),
+        pulse,
+        @variable(model, start=0.1),
+        @variable(model, start=0.1),
+        rotate,
+        slew_rate(GLOBAL_SCANNER[])
+    )
+    for time_var in (res.rise_time1, res.flat_time1, res.diff_time, res.flat_time2, res.fall_time2)
+        @constraint model time_var >= 0
+    end
+    @constraint model res.diff_time <= res.rise_time1
+    @constraint model res.diff_time <= res.fall_time2
+
+    dim = findfirst(v -> !iszero(v), res.orientation)
+    @constraint model qvec(res, nothing, 1)[dim] == qvec(res, 1, nothing)[dim]
+    if !isnothing(spoiler_scale)
+        if spoiler_scale == :min # note that spoiler_scale is inverse of qvec
+            @objective model Min objective_function(model) - qvec(res, nothing, 1)[dim]
+        elseif spoiler_scale == :max
+            @objective model Min objective_function(model) + qvec(res, nothing, 1)[dim]
+        else
+            @constraint model (qvec(res, nothing, 1)[dim] ./ res.orientation[dim]) >= 2π * 1e-3 / spoiler_scale
+        end
+    end
+    set_simple_constraints!(model, res, kwargs)
+
+    max_time = gradient_strength(GLOBAL_SCANNER[]) / slew_rate(GLOBAL_SCANNER[])
+    @constraint model rise_time(res)[1] <= max_time
+    @constraint model fall_time(res)[2] <= max_time
+    return res
+end
+
+duration_trap1(spoilt::SpoiltSliceSelect) = 2 * spoilt.rise_time1 + spoilt.flat_time1 - spoilt.diff_time
+duration_trap2(spoilt::SpoiltSliceSelect) = 2 * spoilt.fall_time2 + spoilt.flat_time2 - spoilt.diff_time
+
+function waveform(spoilt::SpoiltSliceSelect)
+    slew = slew_rate(spoilt)
+    (grad1, grad2, grad3) = [g .* slew_rate(spoilt) for g in all_gradient_strengths(spoilt)]
+    return [
+        ChangingGradientBlock(zeros(3), slew, rise_time(spoilt)[1], spoilt.rotate, nothing),
+        ConstantGradientBlock(grad1, flat_time(spoilt)[1], spoilt.rotate, nothing),
+        ChangingGradientBlock(grad1, -slew, fall_time(spoilt)[1], spoilt.rotate, nothing),
+        ConstantGradientBlock(grad2, duration(spoilt.pulse), spoilt.rotate, nothing),
+        ChangingGradientBlock(grad2, slew, rise_time(spoilt)[2], spoilt.rotate, nothing),
+        ConstantGradientBlock(grad3, flat_time(spoilt)[2], spoilt.rotate, nothing),
+        ChangingGradientBlock(grad3, -slew, fall_time(spoilt)[2], spoilt.rotate, nothing),
+    ]
+end
+interruptions(spoilt::SpoiltSliceSelect) = [(index=4, time=duration_trap1(spoilt) + effective_time(spoilt.pulse), object=spoilt.pulse)]
+
+rise_time(spoilt::SpoiltSliceSelect) = (spoilt.rise_time1, spoilt.fall_time2 - spoilt.diff_time)
+flat_time(spoilt::SpoiltSliceSelect) = (spoilt.flat_time1, spoilt.flat_time2)
+fall_time(spoilt::SpoiltSliceSelect) = (spoilt.rise_time1 - spoilt.diff_time, spoilt.fall_time2)
+duration(spoilt::SpoiltSliceSelect) = sum(rise_time(spoilt)) + sum(flat_time(spoilt)) + sum(flat_time(spoilt)) + duration(spoilt.pulse)
+slew_rate(spoilt::SpoiltSliceSelect) = spoilt.orientation .* spoilt.slew_rate
+inverse_slice_thickness(spoilt::SpoiltSliceSelect) = spoilt.slew_rate * spoilt.diff_time * duration(spoilt.pulse) * 1e3
+function all_gradient_strengths(spoilt::SpoiltSliceSelect)
+    grad1 = spoilt.slew_rate .* rise_time(spoilt)[1]
+    grad2 = grad1 .- spoilt.slew_rate .* flat_time(spoilt)[1]
+    grad3 = spoilt.slew_rate .* fall_time(spoilt)[2]
+    return [grad1, grad2, grad3]
+end
+
+variables(::Type{<:SpoiltSliceSelect}) = [duration, qvec, slew_rate, inverse_slice_thickness, all_gradient_strengths, rise_time, flat_time, fall_time]
+
+
+end
\ No newline at end of file