diff --git a/docs/src/adjust_sequences.md b/docs/src/adjust_sequences.md
index 82ca76a6c85cd5295ba95294a9e0cbbe203de66d..2464ab35d0b92fc9dea8797d38f08ad031aeaa6e 100644
--- a/docs/src/adjust_sequences.md
+++ b/docs/src/adjust_sequences.md
@@ -11,6 +11,7 @@ Post-hoc alterations can be applied to gradients or RF pulses with a specific la
Some example usages are:
- Reduce the RF pulse amplitude by 20% (e.g., to model the effect of transmit bias field): `adjust(sequence, pulse=(scale=0.8, ))`
- Repeat sequence 2 times with different diffusion-weighted gradient orientations (x- and y-direction) and gradient strength reduced by 30%: `adjust(sequence, diffusion=(orientation=[[1., 0., 0], [0., 1., 0.]], scale=0.7))`
+- Repeat the sequence by shifting the excited slice by the given number of millimetres in the slice-select direction: `adjust(sequence, FOV=(shift=[-7.5, -2.5, 2.5, 7.5, -5., 0., 5., 10.]))`. These shifts represent an interleaved acquisition scheme, where the acquired slices/bands are 2.5 mm apart.
- Rotations defined using the [`Rotations.jl`](https://github.com/JuliaGeometry/Rotations.jl) package can be applied to gradient orientations or the field of view. For example, to rotate the field of view by 45 degrees around the y-axis:
```julia
using Rotations
diff --git a/docs/src/sequence_optimisation.md b/docs/src/sequence_optimisation.md
index 0ab0233bdcc45e99a55330c3622fd5980e74b5f1..11bd984e7cbcb1429b27abe21735dc04e705e807 100644
--- a/docs/src/sequence_optimisation.md
+++ b/docs/src/sequence_optimisation.md
@@ -64,7 +64,7 @@ variables.rise_time
variables.slew_rate
variables.slew_rate_norm
variables.slice_thickness
-variables.spoiler_scale
+variables.spoiler
variables.time_to_center
variables.voxel_size
variables.Δ
diff --git a/src/MRIBuilder.jl b/src/MRIBuilder.jl
index 6d55ef6c5abc749a18f8deea5f1d7b3917c316f0..f49f5f46786e2dc0a2233875b89fd9002bd2d08d 100644
--- a/src/MRIBuilder.jl
+++ b/src/MRIBuilder.jl
@@ -28,8 +28,8 @@ export variables, make_generic, @defvar, get_pulse, get_readout, get_pathway, ge
import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
-import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, 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, AlternativeBlocks, match_blocks!, get_index_single_TR, iter_blocks, iter_instant_gradients, iter_instant_pulses
+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
import .Pathways: Pathway
export Pathway
diff --git a/src/build_sequences.jl b/src/build_sequences.jl
index 64eba7a869562a17d95586202f95302ace2d11f8..c9a715a7a2daff57d504b2df2b7882e89582c54b 100644
--- a/src/build_sequences.jl
+++ b/src/build_sequences.jl
@@ -1,5 +1,5 @@
module BuildSequences
-import JuMP: Model, optimizer_with_attributes, optimize!, AbstractJuMPScalar, value, solution_summary, termination_status, LOCALLY_SOLVED, OPTIMAL, num_variables, all_variables, set_start_value, ALMOST_LOCALLY_SOLVED, objective_value, INVALID_MODEL, @variable, @objective, @constraint
+using JuMP
import Ipopt
import Juniper
import ..Scanners: Scanner, gradient_strength, Default_Scanner
@@ -95,6 +95,10 @@ function build_sequence(f::Function, scanner::Union{Nothing, Scanner}, model::Tu
end
end
+function number_equality_constraints(model::Model)
+ sum([num_constraints(model, expr, comp) for (expr, comp) in JuMP.list_of_constraint_types(model) if comp <: MOI.EqualTo])
+end
+
function optimise_with_cost_func(jump_model::Model, cost_func, n_attempts)
@objective jump_model Min cost_func
min_objective = Inf
@@ -107,11 +111,10 @@ function optimise_with_cost_func(jump_model::Model, cost_func, n_attempts)
set_start_value(var, v)
end
end
- optimize!(jump_model)
- while termination_status(jump_model) == INVALID_MODEL
+ for _ in num_variables(jump_model):number_equality_constraints(jump_model)
@variable(jump_model)
- optimize!(jump_model)
end
+ optimize!(jump_model)
if termination_status(jump_model) in (LOCALLY_SOLVED, OPTIMAL)
if objective_value(jump_model) < min_objective
min_objective = objective_value(jump_model)
@@ -126,7 +129,7 @@ function optimise_with_cost_func(jump_model::Model, cost_func, n_attempts)
end
end
-function build_sequence(f::Function, scanner::Union{Nothing, Scanner}, optimiser_constructor; optimise=true, n_attempts=100, kwargs...)
+function build_sequence(f::Function, scanner::Union{Nothing, Scanner}, optimiser_constructor; optimise=true, n_attempts=30, kwargs...)
if optimise || GLOBAL_MODEL[] == IGNORE_MODEL
model = (
Model(optimizer_with_attributes(optimiser_constructor, [string(k) => v for (k, v) in kwargs]...)),
diff --git a/src/components/abstract_types.jl b/src/components/abstract_types.jl
index bcbfe7dd50e461108f94682bd1715d81eb0d4d3c..88862b66cab4e5b35f6a60ebe2daccd8feac39ed 100644
--- a/src/components/abstract_types.jl
+++ b/src/components/abstract_types.jl
@@ -1,5 +1,5 @@
module AbstractTypes
-import ...Variables: AbstractBlock, variables, adjustable, gradient_orientation, @defvar
+import ...Variables: AbstractBlock, variables, adjust_groups, gradient_orientation, @defvar
"""
Super-type for all individual components that form an MRI sequence (i.e., RF pulse, gradient waveform, or readout event).
@@ -110,8 +110,8 @@ It should be infinite if the component is linear.
split_timestep(comp_tuple::Tuple{<:Number, <:EventComponent}, precision::Number) = split_timestep(comp_tuple[2], precision)
-adjustable(::RFPulseComponent) = :pulse
-adjustable(::GradientWaveform) = :gradient
+adjust_groups(p::RFPulseComponent) = [p.group, :pulse]
+adjust_groups(g::GradientWaveform) = [g.group, :gradient]
gradient_orientation(gw::GradientWaveform{1}) = gw.orientation
diff --git a/src/components/gradient_waveforms/no_gradient_blocks.jl b/src/components/gradient_waveforms/no_gradient_blocks.jl
index 31a534ed79c512bd1cd008efe71148df261cbb7e..3b0eec892f90603cbdea2ac0dbc459c374f2a4ed 100644
--- a/src/components/gradient_waveforms/no_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/no_gradient_blocks.jl
@@ -1,6 +1,6 @@
module NoGradientBlocks
import StaticArrays: SVector, SMatrix
-import ....Variables: VariableType, get_free_variable, adjustable, variables, @defvar
+import ....Variables: VariableType, get_free_variable, adjust_groups, variables, @defvar
import ...AbstractTypes: GradientWaveform
import ..ChangingGradientBlocks: split_gradient
@@ -35,6 +35,6 @@ function split_gradient(ngb::NoGradient, times::VariableType...)
return [NoGradient(d) for d in durations]
end
-adjustable(::NoGradient) = :false
+adjust_groups(::NoGradient) = Symbol[]
end
\ No newline at end of file
diff --git a/src/components/instant_gradients.jl b/src/components/instant_gradients.jl
index da550bfabdb2f9598a9d5d1e36a38ac5beda8816..ddca7f0d0a77f16ca1067ac36f59d4ebff640b32 100644
--- a/src/components/instant_gradients.jl
+++ b/src/components/instant_gradients.jl
@@ -1,7 +1,7 @@
module InstantGradients
import StaticArrays: SVector, SMatrix
import JuMP: @constraint
-import ...Variables: @defvar, VariableType, variables, get_free_variable, set_simple_constraints!, make_generic, adjust_internal, adjustable, gradient_orientation, apply_simple_constraint!
+import ...Variables: @defvar, VariableType, variables, get_free_variable, set_simple_constraints!, make_generic, adjust_internal, adjust_groups, gradient_orientation, apply_simple_constraint!
import ..AbstractTypes: EventComponent, GradientWaveform
"""
@@ -15,7 +15,7 @@ If the `orientation` is set an [`InstantGradient1D`](@ref) is returned, otherwis
## Variables
- [`variables.qvec`](@ref): Spatial frequency on which spins will be dephased due to this pulsed gradient in rad/um.
-- [`variables.spoiler_scale`](@ref): Length-scale on which spins will be dephased by exactly 2Ï€ in mm.
+- [`variables.spoiler`](@ref): Length-scale on which spins will be dephased by exactly 2Ï€ in mm.
"""
abstract type InstantGradient{N} <: EventComponent end
@@ -62,7 +62,7 @@ end
make_generic(ig::InstantGradient) = ig
-adjustable(::InstantGradient) = :gradient
+adjust_groups(g::InstantGradient) = [g.group, :gradient]
gradient_orientation(ig::InstantGradient{1}) = ig.orientation
diff --git a/src/containers/alternatives.jl b/src/containers/alternatives.jl
index 8847d5f3e8f911943a812fbca9d093c97a2bae63..791e3fb1d811200cf495bdd26614be9c99f21375 100644
--- a/src/containers/alternatives.jl
+++ b/src/containers/alternatives.jl
@@ -4,6 +4,17 @@ import ..Abstract: ContainerBlock
import ...BuildSequences: fixed
import ...Variables: @defvar, make_generic, apply_simple_constraint!
+"""
+Parent type for all blocks that can take different MR sequence components between multiple repetitions of the sequence.
+
+They can be extended into their individual components using `adjust(<name>=:all)`.
+
+Each subtype of [`AbstractAlternativeBlock`](@ref) needs to implement two methods:
+- [`get_alternatives_name`](@ref): returns the `name` used to identify this block in `adjust`
+- [`get_alternatives_options`](@ref): returns a dictionary mapping the name of the different options to a [`ContainerBlock`](@ref) with the actual sequence building block.
+"""
+abstract type AbstractAlternativeBlocks <: ContainerBlock end
+
"""
AlternativeBlocks(name, blocks)
@@ -11,19 +22,35 @@ 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.
+The `name` is a symbol that is used to identify this `AlternativeBlocks` in the broader sequence (as in `adjust`).
"""
-struct AlternativeBlocks <: ContainerBlock
+struct AlternativeBlocks <: AbstractAlternativeBlocks
name :: Symbol
options :: Dict{Any, <:ContainerBlock}
end
AlternativeBlocks(name::Symbol, options_vector::AbstractVector) = AlternativeBlocks(name, Dict(index => value for (index, value) in enumerate(options_vector)))
-Base.getindex(alt::AlternativeBlocks, index) = alt.options[index]
-Base.length(alt::AlternativeBlocks) = length(alt.options)
+"""
+ get_alternatives_name(alternative_block)
+
+Get the name with which any [`AbstractAlternativeBlocks`](@ref) will be identified in a call to `adjust`.
+"""
+get_alternatives_name(alt::AlternativeBlocks) = alt.name
-@defvar duration(alt::AlternativeBlocks) = maximum(variables.duration.(values(alt.options)))
+"""
+ get_alternatives_options(alternative_block)
+
+Get the options available for a [`AbstractAlternativeBlocks`](@ref).
+"""
+get_alternatives_options(alt::AlternativeBlocks) = alt.options
+
+Base.getindex(alt::AbstractAlternativeBlocks, index) = get_alternatives_options(alt)[index]
+Base.length(alt::AbstractAlternativeBlocks) = length(get_alternatives_options(alt))
+Base.keys(alt::AbstractAlternativeBlocks) = keys(get_alternatives_options(alt))
+Base.values(alt::AbstractAlternativeBlocks) = values(get_alternatives_options(alt))
+
+@defvar duration(alt::AbstractAlternativeBlocks) = maximum(variables.duration.(values(alt.options)))
"""
match_blocks!(alternatives, function)
@@ -32,8 +59,11 @@ Matches the outcome of given `function` on each of the building blocks in [`Alte
For example, `match_blocks!(alternatives, duration)` will ensure that all the alternative building blocks have the same duration.
"""
-function match_blocks!(alternatives::AlternativeBlocks, func)
- options = [values(alternatives.options)...]
+function match_blocks!(alternatives::AbstractAlternativeBlocks, func)
+ options = [values(get_alternatives_options(alternatives.options))...]
+ if length(options) <= 1
+ return
+ end
baseline = func(options[1])
for other_block in options[2:end]
apply_simple_constraint!(func(other_block), baseline)
@@ -41,7 +71,7 @@ function match_blocks!(alternatives::AlternativeBlocks, func)
end
fixed(alt::AlternativeBlocks) = AlternativeBlocks(alt.name, Dict(key=>fixed(value) for (key, value) in alt.options))
-make_generic(alt::AlternativeBlocks) = AlternativeBlocks(alt.name, Dict(key=>make_generic(value) for (key, value) in alt.options))
+make_generic(alt::AbstractAlternativeBlocks) = AlternativeBlocks(get_alternatives_name(alt), get_alternatives_options(alt))
end
\ No newline at end of file
diff --git a/src/containers/containers.jl b/src/containers/containers.jl
index 96198e5d98ad698ee789af1dacbcc23a85f1555e..55c92a1679bc47674301ea34abbb5227e153f97f 100644
--- a/src/containers/containers.jl
+++ b/src/containers/containers.jl
@@ -7,6 +7,6 @@ include("alternatives.jl")
import .Abstract: ContainerBlock, start_time, end_time, iter_blocks, iter_instant_gradients, iter_instant_pulses
import .BuildingBlocks: BaseBuildingBlock, BuildingBlock, Wait, waveform, waveform_sequence, events, ndim_grad
import .BaseSequences: BaseSequence, Sequence, nrepeat, get_index_single_TR
-import .Alternatives: AlternativeBlocks, match_blocks!
+import .Alternatives: AbstractAlternativeBlocks, get_alternatives_name, get_alternatives_options, AlternativeBlocks, match_blocks!
end
\ No newline at end of file
diff --git a/src/parts/epi_readouts.jl b/src/parts/epi_readouts.jl
index ed89c70cc9701dd972bfa8119bf9d728db6c54bd..0757d2e311cacc89f9835f0951f7df0d0a960754 100644
--- a/src/parts/epi_readouts.jl
+++ b/src/parts/epi_readouts.jl
@@ -2,7 +2,7 @@ module EPIReadouts
import ...Containers: BaseSequence, get_index_single_TR
import ..Trapezoids: Trapezoid, opposite_kspace_lines, LineReadout
import ...Components: ADC
-import ...Variables: get_free_variable, VariableType, set_simple_constraints!, get_readout, apply_simple_constraint!, variables, @defvar
+import ...Variables: get_free_variable, VariableType, set_simple_constraints!, get_readout, apply_simple_constraint!, variables, @defvar, add_cost_function!
import ...Pathways: PathwayWalker, update_walker_till_time!, walk_pathway!
"""
@@ -58,6 +58,12 @@ function EPIReadout(; resolution::AbstractVector{<:Integer}, recenter=false, gro
res.blips[shift] = Trapezoid(qvec=[0., shift * res.ky_step, 0.], group=group)
end
set_simple_constraints!(res, vars)
+ add_cost_function!(sum(variables.duration.([
+ res.start_gradient,
+ res.positive_line,
+ res.negative_line,
+ values(res.blips)...,
+ ])))
return res
end
diff --git a/src/parts/helper_functions.jl b/src/parts/helper_functions.jl
index 66a307aa81ae6aeb6595fb0e73b364a8f5706276..b89f7bfb4925de84dadf24634a2445056ccccb03 100644
--- a/src/parts/helper_functions.jl
+++ b/src/parts/helper_functions.jl
@@ -135,7 +135,7 @@ function refocus_pulse(; flip_angle=180, phase=0., frequency=0., shape=nothing,
return SliceSelect(pulse; duration=:min, slice_thickness=slice_thickness, orientation=orientation, group=:FOV)
end
else
- res = SpoiltSliceSelect(pulse; orientation=orientation, duration=:min, group=:FOV, slice_thickness=slice_thickness, spoiler_scale=spoiler)
+ res = SpoiltSliceSelect(pulse; orientation=orientation, duration=:min, group=:FOV, slice_thickness=slice_thickness, spoiler=spoiler)
return res
end
end
@@ -285,7 +285,7 @@ function interpret_image_size(fov, resolution, voxel_size, slice_thickness)
end
"""
- gradient_spoiler(; optimise=false, orientation=[0, 0, 1], rotate=:FOV, scale=:spoiler, spoiler_scale=1., duration=:min, variables...)
+ gradient_spoiler(; optimise=false, orientation=[0, 0, 1], rotate=:FOV, scale=:spoiler, spoiler=1., duration=:min, variables...)
Returns two DWI gradients that are guaranteed to cancel each other out.
@@ -297,12 +297,12 @@ Returns two DWI gradients that are guaranteed to cancel each other out.
- `scanner`: Used for testing. Do not set this parameter at this level (instead set it for the total sequence using [`build_sequence`](@ref)).
## Variables
-- [`variables.spoiler_scale`](@ref): maximum spoiler scale (before applying any reductions due to `scale`).
+- [`variables.spoiler`](@ref): maximum spoiler scale (before applying any reductions due to `scale`).
- Any other parameters expected by [`Trapezoid`](@ref).
"""
-function gradient_spoiler(; optimise=false, scanner=nothing, orientation=[0, 0, 1], group=:FOV, spoiler_scale=1., duration=:min, variables...)
+function gradient_spoiler(; optimise=false, scanner=nothing, orientation=[0, 0, 1], group=:FOV, spoiler=1., duration=:min, variables...)
build_sequence(scanner; optimise=optimise) do
- Trapezoid(; orientation=orientation, group=group, spoiler_scale=spoiler_scale, duration=duration, variables...)
+ Trapezoid(; orientation=orientation, group=group, spoiler=spoiler, duration=duration, variables...)
end
end
diff --git a/src/parts/spoilt_slice_selects.jl b/src/parts/spoilt_slice_selects.jl
index adef2102f4d17c0e28312183a7d11977cd810a5f..b529123c4581d99079e5666143e7b6a86c2ff99d 100644
--- a/src/parts/spoilt_slice_selects.jl
+++ b/src/parts/spoilt_slice_selects.jl
@@ -4,13 +4,13 @@ import LinearAlgebra: norm
import StaticArrays: SVector
import JuMP: @constraint, @objective, objective_function
import ...BuildSequences: global_scanner
-import ...Variables: VariableType, get_pulse, set_simple_constraints!, variables, @defvar, gradient_orientation
+import ...Variables: VariableType, get_pulse, set_simple_constraints!, variables, @defvar, gradient_orientation, adjust, adjust_groups, adjust_internal, get_free_variable, apply_simple_constraint!
import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent
import ...Containers: BaseBuildingBlock
"""
- SpoiltSliceSelect(pulse; parameters, variables...)
+ SpoiltSliceSelect(pulse; parameters..., variables...)
Adds slice selection to the `pulse` and surrounds it with spoiler gradients.
@@ -21,7 +21,7 @@ Adds slice selection to the `pulse` and surrounds it with spoiler gradients.
## 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. This sets the minimum spoiling. If this spoiling level is not achieved by the slice-select gradient alone, then there will be additional gradients added.
+- `spoiler`: length scale on which the spoilers achieve 2Ï€ dephasing in mm. This sets the minimum spoiling. If this spoiling level is not achieved by the slice-select gradient alone, then there will be additional gradients added.
"""
struct SpoiltSliceSelect <: BaseBuildingBlock
orientation :: SVector{3, Float64}
@@ -35,7 +35,7 @@ struct SpoiltSliceSelect <: BaseBuildingBlock
slew_rate :: Number
end
-function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group=nothing, slice_thickness=nothing, kwargs...)
+function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group=:FOV, slice_thickness=nothing, spoiler, kwargs...)
res = nothing
if slice_thickness isa Number && isinf(slice_thickness)
rise_time_var = get_free_variable(nothing)
@@ -49,7 +49,7 @@ function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group
flat_time_var,
rise_time_var,
group,
- slew_rate(global_scanner()),
+ global_scanner().slew_rate,
)
for time_var in (rise_time_var, flat_time_var)
apply_simple_constraint!(time_var, :>=, 0)
@@ -64,21 +64,23 @@ function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group
get_free_variable(nothing; start=0.1),
get_free_variable(nothing; start=0.1),
group,
- slew_rate(global_scanner()),
+ global_scanner().slew_rate,
)
for time_var in (res.rise_time1, res.flat_time1, res.diff_time, res.flat_time2, res.fall_time2)
apply_simple_constraint!(time_var, :>=, 0)
end
apply_simple_constraint!(res.diff_time, :<=, res.rise_time1)
- apply_simple_constraint!(res.diff_time, :<=, res.rise_time2)
- apply_simple_constraint!(qvec(res, nothing, :pulse), qvec(res, :pulse, nothing))
+ apply_simple_constraint!(res.diff_time, :<=, res.fall_time2)
+ apply_simple_constraint!(variables.qvec(res, nothing, :pulse), variables.qvec(res, :pulse, nothing))
+ apply_simple_constraint!(variables.inverse_slice_thickness(res), 1/slice_thickness)
end
set_simple_constraints!(res, kwargs)
+ apply_simple_constraint!(variables.spoiler(res), :<=, spoiler)
- max_time = gradient_strength(global_scanner()) / res.slew_rate
- apply_simple_constraint!(rise_time(res)[1], :<=, max_time)
- apply_simple_constraint!(fall_time(res)[2], :<=, max_time)
+ max_time = global_scanner().gradient / res.slew_rate
+ apply_simple_constraint!(variables.rise_time(res)[1], :<=, max_time)
+ apply_simple_constraint!(variables.fall_time(res)[2], :<=, max_time)
return res
end
@@ -87,36 +89,37 @@ gradient_orientation(spoilt::SpoiltSliceSelect) = spoilt.orientation
@defvar begin
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
+ effective_time(spoilt::SpoiltSliceSelect) = variables.effective_time(spoilt, Val(:pulse))
end
Base.keys(::SpoiltSliceSelect) = Val.((:rise1, :flat1, :fall1, :flat_pulse, :pulse, :rise2, :flat2, :fall2))
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise1}) = ChangingGradient(0., variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.rise_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat1}) = ConstantGradient(variables.slew_rate(spoilt) * variables.rise_time(spoilt)[1], gradient_orientation(spoilt), variables.flat_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall1}) = ChangingGradient(variables.slew_rate(spoilt) * variables.rise_time(spoilt)[1], -variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.fall_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat_pulse}) = ConstantGradient(variables.slew_rate(spoilt) * spoilt.diff_time, gradient_orientation(spoilt), variables.duration(spoilt.pulse), spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:pulse}) = spoilt.pulse
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise2}) = ChangingGradient(variables.slew_rate(spoilt) * spoilt.diff_time, variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.rise_time(spoilt)[2], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat2}) = ConstantGradient(variables.slew_rate(spoilt) * variables.fall_time(spoilt)[2], gradient_orientation(spoilt), variables.flat_time(spoilt)[2], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall2}) = ChangingGradient(variables.slew_rate(spoilt) * variables.fall_time(spoilt)[2], -variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.fall_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise1}) = ChangingGradient(zero(SVector{3, Float64}), variables.slew_rate(spoilt), variables.rise_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat1}) = ConstantGradient(variables.slew_rate(spoilt) .* variables.rise_time(spoilt)[1], variables.flat_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall1}) = ChangingGradient(variables.slew_rate(spoilt) .* variables.rise_time(spoilt)[1], -variables.slew_rate(spoilt), variables.fall_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat_pulse}) = ConstantGradient(variables.slew_rate(spoilt) .* spoilt.diff_time, variables.duration(spoilt.pulse), spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:pulse}) = (0., spoilt.pulse)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise2}) = ChangingGradient(variables.slew_rate(spoilt) .* spoilt.diff_time, variables.slew_rate(spoilt), variables.rise_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat2}) = ConstantGradient(variables.slew_rate(spoilt) .* variables.fall_time(spoilt)[2], variables.flat_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall2}) = ChangingGradient(variables.slew_rate(spoilt) .* variables.fall_time(spoilt)[2], -variables.slew_rate(spoilt), variables.fall_time(spoilt)[2], spoilt.group)
@defvar begin
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)
- slew_rate(spoilt::SpoiltSliceSelect) = spoilt.slew_rate
+ slew_rate(spoilt::SpoiltSliceSelect) = spoilt.slew_rate .* spoilt.orientation
end
@defvar begin
- duration(spoilt::SpoiltSliceSelect) = sum(variables.rise_time(spoilt)) + sum(variables.flat_time(spoilt)) + sum(variables.flat_time(spoilt)) + variables.duration(spoilt.pulse)
- inverse_slice_thickness(spoilt::SpoiltSliceSelect) = variables.slew_rate(spoilt) * spoilt.diff_time * variables.duration(spoilt.pulse) * 1e3
- gradient_strength(spoilt::SpoiltSliceSelect) = variables.slew_rate(spoilt) * max(spoilt.rise_time1, spoilt.fall_time2)
+ duration(spoilt::SpoiltSliceSelect) = sum(variables.rise_time(spoilt)) + sum(variables.flat_time(spoilt)) + sum(variables.fall_time(spoilt)) + variables.duration(spoilt.pulse)
+ inverse_slice_thickness(spoilt::SpoiltSliceSelect) = spoilt.slew_rate * spoilt.diff_time * 1e3 * variables.inverse_bandwidth(spoilt.pulse)
+ gradient_strength(spoilt::SpoiltSliceSelect) = variables.slew_rate(spoilt) .* max(spoilt.rise_time1, spoilt.fall_time2)
end
get_pulse(spoilt::SpoiltSliceSelect) = spoilt.pulse
@defvar function all_gradient_strengths(spoilt::SpoiltSliceSelect)
grad1 = spoilt.slew_rate * variables.rise_time(spoilt)[1]
- grad2 = grad1 - spoilt.slew_rate * variables.flat_time(spoilt)[1]
+ grad2 = grad1 - spoilt.slew_rate * variables.fall_time(spoilt)[1]
grad3 = spoilt.slew_rate * variables.fall_time(spoilt)[2]
return [grad1, grad2, grad3]
end
@@ -134,6 +137,25 @@ variables.all_gradient_strengths
Returns the time of the [`SpoiltSliceSelect`](@ref) to return to zero.
"""
-variables.all_gradient_strengths
+variables.fall_time
+
+adjust_groups(ss::SpoiltSliceSelect) = [ss.group, :slice]
+
+function adjust_internal(ss::SpoiltSliceSelect; shift=0.)
+ return SpoiltSliceSelect(
+ ss.orientation,
+ ss.rise_time1,
+ ss.flat_time1,
+ ss.diff_time,
+ adjust(
+ ss.pulse;
+ pulse=(frequency=variables.all_gradient_strengths(ss)[2] * shift * 1e3, )
+ ),
+ ss.flat_time2,
+ ss.fall_time2,
+ ss.group,
+ ss.slew_rate,
+ )
+end
end
\ No newline at end of file
diff --git a/src/parts/trapezoids.jl b/src/parts/trapezoids.jl
index 8a7eebb8b618e93a0ee1e1a5bfa75ccdd0bab5e3..084accce7285f0b85656a0ceb4bbbb9926262412 100644
--- a/src/parts/trapezoids.jl
+++ b/src/parts/trapezoids.jl
@@ -6,7 +6,7 @@ module Trapezoids
import JuMP: @constraint
import StaticArrays: SVector
import LinearAlgebra: norm
-import ...Variables: variables, @defvar, scanner_constraints!, get_free_variable, set_simple_constraints!, gradient_orientation, VariableType, get_gradient, get_pulse, get_readout, adjustable, adjust_internal, apply_simple_constraint!, add_cost_function!
+import ...Variables: variables, @defvar, scanner_constraints!, get_free_variable, set_simple_constraints!, gradient_orientation, VariableType, get_gradient, get_pulse, get_readout, adjust, adjust_groups, adjust_internal, apply_simple_constraint!, add_cost_function!
import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent, ADC
import ...Containers: BaseBuildingBlock
@@ -146,7 +146,7 @@ variables.δ
@defvar qvec(g::BaseTrapezoid, ::Nothing, ::Nothing) = variables.δ(g) .* variables.gradient_strength(g) .* 2π
-adjustable(::BaseTrapezoid) = :gradient
+adjust_groups(t::Trapezoid) = [t.group, :gradient]
function adjust_internal(trap::Trapezoid1D; orientation=nothing, scale=1., rotation=nothing)
if !isnothing(orientation) && !isnothing(rotation)
@@ -192,12 +192,14 @@ Parameters and variables are identical as for [`Trapezoid`](@ref) with the addit
struct SliceSelect{N} <: BaseTrapezoid{N}
trapezoid :: Trapezoid{N}
pulse :: RFPulseComponent
+ group :: Symbol
end
-function SliceSelect(pulse::RFPulseComponent; orientation=nothing, rise_time=nothing, group=nothing, slew_rate=nothing, vars...)
+function SliceSelect(pulse::RFPulseComponent; orientation=nothing, rise_time=nothing, group=:FOV, slew_rate=nothing, vars...)
res = SliceSelect(
- Trapezoid(; orientation=orientation, rise_time=rise_time, flat_time=variables.duration(pulse), group=group, slew_rate=slew_rate),
- pulse
+ Trapezoid(; orientation=orientation, rise_time=rise_time, flat_time=variables.duration(pulse), slew_rate=slew_rate, group=nothing),
+ pulse,
+ group
)
set_simple_constraints!(res, vars)
return res
@@ -221,6 +223,19 @@ get_pulse(ss::SliceSelect) = ss.pulse
get_gradient(ss::SliceSelect) = ss.trapezoid
@defvar effective_time(ss::SliceSelect) = variables.effective_time(ss, :pulse)
+adjust_groups(ss::SliceSelect) = [ss.group, :slice]
+function adjust_internal(ss::SliceSelect; shift)
+ return SliceSelect(
+ ss.trapezoid,
+ adjust(
+ ss.pulse;
+ pulse=(frequency=variables.gradient_strength_norm(ss) * shift * 1e3,)
+ ),
+ ss.group,
+ )
+end
+
+
"""
LineReadout(adc; ramp_overlap=1., orientation=nothing, group=nothing, variables...)
diff --git a/src/post_hoc.jl b/src/post_hoc.jl
index 9a8d7f14c7f44f20d4cc5c21a8989adc436bbde3..b071ad3250226b16d9ce530a0d0c15d24d23edbc 100644
--- a/src/post_hoc.jl
+++ b/src/post_hoc.jl
@@ -3,10 +3,12 @@ Define post-fitting adjustments of the sequences
"""
module PostHoc
-import ..Variables: AbstractBlock, adjust_internal, adjustable
+import ..Variables: AbstractBlock, adjust_internal, adjust_groups, adjust
import ..Components: GradientWaveform, RFPulseComponent, BaseComponent, NoGradient
import ..Containers: ContainerBlock, Sequence, Wait
+const UsedNamesType = Dict{Symbol, Set{Symbol}}
+
"""
adjust(block; kwargs...)
@@ -34,7 +36,11 @@ To affect all gradients or pulses, use `gradient=` or `pulse`, e.g.
will divide the amplitude of all RV pulses by two.
"""
function adjust(block::AbstractBlock; merge=true, kwargs...)
- used_names = Set{Symbol}()
+ invalid_type = Set(key for (key, value) in pairs(kwargs) if !(value isa NamedTuple))
+ if length(invalid_type) > 0
+ error("All `adjust` keywords except for merge should be a NamedTuple, like (scale=3, ). This is not the case for: $(invalid_type).")
+ end
+ used_names = UsedNamesType()
n_adjust, kwargs_list = adjust_kwargs_list(; kwargs...)
if isnothing(n_adjust)
res = adjust_helper(block, used_names; kwargs_list[1]...)
@@ -48,11 +54,22 @@ function adjust(block::AbstractBlock; merge=true, kwargs...)
end
end
unused_names = filter(keys(kwargs)) do key
- !(key in used_names)
+ !(key in keys(used_names))
end
if length(unused_names) > 0
@warn "Some group/type names were not used in call to `MRIBuilder.adjust`, namely: $(unused_names)."
end
+ for group_name in keys(kwargs)
+ if group_name in unused_names
+ continue
+ end
+ unused_keys = filter(keys(kwargs[group_name])) do key
+ !(key in used_names[group_name])
+ end
+ if length(unused_keys) > 0
+ @warn "Some keywords provided for group `$(group_name)` were not used, namely: $(unused_keys)."
+ end
+ end
res
end
@@ -90,34 +107,35 @@ function adjust_kwargs_list(; kwargs...)
return (n_adjust, kwargs_list)
end
-
-function adjust_helper(block::AbstractBlock, used_names::Set{Symbol}; gradient=(), pulse=(), kwargs...)
+function adjust_helper(block::AbstractBlock, used_names::UsedNamesType; kwargs...)
params = []
- adjust_type = adjustable(block)
- if adjust_type == :false
- for prop_name in propertynames(block)
- push!(params, adjust_helper(getproperty(block, prop_name), used_names; gradient=gradient, pulse=pulse, kwargs...))
- end
- return typeof(block)(params...)
- else
- if !isnothing(block.group) && (block.group in keys(kwargs))
- push!(used_names, block.group)
- return adjust_internal(block; kwargs[block.group]...)
- elseif adjust_type == :gradient
- push!(used_names, :gradient)
- return adjust_internal(block; gradient...)
- elseif adjust_type == :pulse
- push!(used_names, :pulse)
- return adjust_internal(block; pulse...)
+
+ for prop_name in propertynames(block)
+ push!(params, adjust_helper(getproperty(block, prop_name), used_names; kwargs...))
+ end
+ new_block = typeof(block)(params...)
+
+ for group in adjust_groups(new_block)
+ if group in keys(kwargs)
+ if !(group in keys(used_names))
+ used_names[group] = Set{Symbol}()
+ end
+ all_available_kwargs = kwargs[group]
+ use_kwargs = reduce(vcat, Base.kwarg_decl.(methods(adjust_internal, (typeof(new_block), ))))
+ @assert length(use_kwargs) > 0 "Invalid definition of `internal_kwargs` for $(typeof(new_block))"
+ internal_kwargs = Dict(key => value for (key, value) in pairs(all_available_kwargs) if key in use_kwargs)
+ union!(used_names[group], keys(internal_kwargs))
+ return adjust_internal(block; internal_kwargs...)
end
end
+ return new_block
end
-adjust_helper(some_value, used_names::Set{Symbol}; kwargs...) = some_value
-adjust_helper(array_variable::AbstractArray, used_names::Set{Symbol}; kwargs...) = map(array_variable) do v adjust_helper(v, used_names; kwargs...) end
-adjust_helper(dict_variable::AbstractDict, used_names::Set{Symbol}; kwargs...) = typeof(dict_variable)(k => adjust_helper(v, used_names; kwargs...) for (k, v) in pairs(dict_variable))
-adjust_helper(tuple_variable::Tuple, used_names::Set{Symbol}; kwargs...) = map(tuple_variable) do v adjust_helper(v, used_names; kwargs...) end
-adjust_helper(pair:: Pair, used_names::Set{Symbol}; kwargs...) = adjust_helper(pair[1], used_names; kwargs...) => adjust_helper(pair[2], used_names; kwargs...)
+adjust_helper(some_value, used_names::UsedNamesType; kwargs...) = some_value
+adjust_helper(array_variable::AbstractArray, used_names::UsedNamesType; kwargs...) = map(array_variable) do v adjust_helper(v, used_names; kwargs...) end
+adjust_helper(dict_variable::AbstractDict, used_names::UsedNamesType; kwargs...) = typeof(dict_variable)(k => adjust_helper(v, used_names; kwargs...) for (k, v) in pairs(dict_variable))
+adjust_helper(tuple_variable::Tuple, used_names::UsedNamesType; kwargs...) = map(tuple_variable) do v adjust_helper(v, used_names; kwargs...) end
+adjust_helper(pair:: Pair, used_names::UsedNamesType; kwargs...) = adjust_helper(pair[1], used_names; kwargs...) => adjust_helper(pair[2], used_names; kwargs...)
"""
diff --git a/src/variables.jl b/src/variables.jl
index 3e7cac0d8f3ffb38a60af723373265676e06549b..0551d2bf3f7c961cb6165a86262a635c453bf1fc 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -37,21 +37,26 @@ end
Returns the adjusted blocks and add any keywords used in the process to `names_used`.
-This is a helper function used by `adjust`.
+This is a helper function used by [`adjust`](@ref).
+It should be defined for any block that is adjustable (as defined by [`adjust_groups`](@ref)).
"""
function adjust_internal end
"""
- adjustable(block)
+ adjust_groups(block)
-Returns whether a sequence, building block, or component can be adjusted
+Returns an array of keywords in [`adjust`](@ref) that should affect a specfic block.
-Can return one of:
-- `:false`: not adjustable
+If any of these keywords are present in [`adjust`](@ref), then [`adjust_internal`](@ref) will be called.
+
+Some standard keywords are:
- `:gradient`: expects gradient adjustment parameters
- `:pulse`: expects RF pulse adjustment parameters
"""
-adjustable(::AbstractBlock) = :false
+adjust_groups(::AbstractBlock) = Symbol[]
+
+# further defined in post_hoc.jl
+function adjust end
abstract type AnyVariable end
@@ -245,7 +250,7 @@ function def_alternate_variable!(name::Symbol, other_var::Symbol, from_other::Fu
setproperty!(variables, name, AlternateVariable(name, other_var, from_other, to_other, inverse))
end
-def_alternate_variable!(:spoiler_scale, :qval, q->1e-3 * 2Ï€/q, l->1e-3 * 2Ï€/l, true)
+def_alternate_variable!(:spoiler, :qval, q->1e-3 * 2Ï€/q, l->1e-3 * 2Ï€/l, true)
def_alternate_variable!(:qval, :qval_square, sqrt, q -> q * q, false)
def_alternate_variable!(:qval_square, :qvec, qv -> sum(q -> q * q, qv), nothing, false)
@@ -257,13 +262,13 @@ The norm of the [`variables.qvec`](@ref).
variables.qval
"""
- spoiler_scale(gradient)
+ spoiler(gradient)
Spatial scale in mm over which the spoiler gradient will dephase by 2Ï€.
Automatically computed based on [`variables.qvec`](@ref).
"""
-variables.spoiler_scale
+variables.spoiler
for vec_variable in [:gradient_strength, :slew_rate]
diff --git a/test/test_post_hoc.jl b/test/test_post_hoc.jl
index 657d6b0e986ec25d87e7730ea5ba1a0a2136ee7b..072f3a27872f73df92b416d988ddeee9d5b328ed 100644
--- a/test/test_post_hoc.jl
+++ b/test/test_post_hoc.jl
@@ -53,5 +53,15 @@
@test all(variables.qvec(new_dwi[:gradient]) .≈ [qval_orig/√2, qval_orig/√2, 0.])
end
end
+ @testset "Slice selection" begin
+ dwi = DWI(bval=1., slice_thickness=2., refocus=(spoiler=0.5, ))
+ @test variables.frequency(dwi) == (excitation=0., refocus=0.)
+ new_dwi = adjust(dwi, slice=(shift=2., ))
+ @test variables.frequency(new_dwi).excitation ≈ variables.bandwidth(dwi).excitation
+ @test variables.frequency(new_dwi).refocus ≈ variables.bandwidth(dwi).refocus
+ another_dwi = adjust(new_dwi, slice=(shift=-3., ))
+ @test variables.frequency(another_dwi).excitation ≈ -variables.bandwidth(dwi).excitation * 0.5
+ @test variables.frequency(another_dwi).refocus ≈ -variables.bandwidth(dwi).refocus * 0.5
+ end
end
end