Skip to content
Snippets Groups Projects
Verified Commit 6ae11cdf authored by Michiel Cottaar's avatar Michiel Cottaar
Browse files

Refactor abstract containers away from concrete containers

parent f66b7fd1
No related branches found
No related tags found
No related merge requests found
Showing
with 220 additions and 233 deletions
......@@ -7,12 +7,9 @@ include("scanners.jl")
include("build_sequences.jl")
include("variables.jl")
include("components/components.jl")
include("container_blocks.jl")
include("all_building_blocks/all_building_blocks.jl")
include("all_sequences/all_sequences.jl")
include("alternatives.jl")
include("containers/containers.jl")
include("pathways.jl")
include("helper_functions.jl")
include("parts/parts.jl")
#include("printing.jl")
import .BuildSequences: build_sequence, global_model, global_scanner
......@@ -24,26 +21,17 @@ 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
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
import .ContainerBlocks: ContainerBlock, start_time, end_time
export start_time, end_time
import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
import .AllBuildingBlocks: waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Trapezoid, SliceSelect, LineReadout, SpoiltSliceSelect, Wait
export waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Trapezoid, SliceSelect, LineReadout, SpoiltSliceSelect, Wait
import .AllSequences: BaseSequence, nrepeat, Sequence, EPIReadout
export BaseSequence, nrepeat, Sequence, EPIReadout
import .Alternatives: AlternativeBlocks, match_blocks!
export AlternativeBlocks, match_blocks!
import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR
export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR
import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat
export Pathway, duration_transverse, duration_dephase, bval, bmat
import .HelperFunctions: excitation_pulse, refocus_pulse, epi_readout, single_line_readout
export excitation_pulse, refocus_pulse, epi_readout, single_line_readout
import .Parts: excitation_pulse, refocus_pulse, epi_readout, single_line_readout, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout
export excitation_pulse, refocus_pulse, epi_readout, single_line_readout, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout
import JuMP: @constraint, @objective, objective_function, value, Model
export @constraint, @objective, objective_function, value, Model
......
module BuildingBlocks
import LinearAlgebra: norm
import ..BaseBuildingBlocks: BaseBuildingBlock, events, waveform_sequence
import ...Variables: VariableType, duration, make_generic, get_pulse, get_readout, scanner_constraints!
import ...Components: BaseComponent, DelayedEvent, RFPulseComponent, ReadoutComponent
"""
BuildingBlock(waveform, events; duration=nothing, orientation=nothing, group)
Generic [`BaseBuildingBlock`](@ref) that can capture any overlapping gradients, RF pulses, and/or readouts.
The gradients cannot contain any free variables.
## Arguments
- `waveform`: Sequence of 2-element tuples with (time, (Gx, Gy, Gz)). If `orientation` is set then the tuple is expected to look like (time, G). This cannot contain any free variables.
- `events`: Sequence of 2-element tuples with (index, pulse/readout). The start time of the pulse/readout at the start of the gradient waveform element with index `index` (use [`DelayedEvent`](@ref) to make this earlier or later).
- `duration`: duration of this `BuildingBlock`. If not set then it will be assumed to be the time of the last element in `waveform`.
- `orientation`: orientation of the gradients in the waveform. If not set, then the full gradient vector should be given explicitly.
- `group`: group of the gradient waveform
"""
struct BuildingBlock <: BaseBuildingBlock
parts :: Vector{<:BaseComponent}
function BuildingBlock(parts::AbstractVector{<:BaseComponent})
res = new(duration, parts)
for (_, part) in waveform_sequence(parts)
scanner_constraints!(part)
end
return res
end
end
function BuildingBlock(waveform::AbstractVector, events::AbstractVector; duration=nothing, orientation=nothing, group=nothing)
events = Any[events...]
waveform = Any[waveform...]
ndim = isnothing(orientation) ? 1 : 3
zero_grad = isnothing(orientation) ? zeros(3) : 0.
if length(waveform) == 0 || waveform[1][1] > 0.
pushfirst!(waveform, (0., zero_grad))
events = [(i+1, e) for (i, e) in events]
end
if isnothing(duration)
duration = waveform[end][1]
end
if !(duration waveform[end][1])
@assert duration > waveform[end][1]
push!(waveform, (duration, zero_grad))
end
components = BaseComponent[]
for (index_grad, ((prev_time, prev_grad), (time, grad))) in enumerate(zip(waveform[1:end-1], waveform[2:end]))
duration = time - prev_time
if norm(prev_grad) <= 1e-12 && norm(grad) <= 1e-12
push!(components, NoGradient{ndim}(duration))
elseif norm(prev_grad) norm(grad)
push!(components, ConstantGradient(prev_grad, orientation, duration, group))
else
push!(components, ChangingGradient(prev_grad, (grad .- prev_grad) ./ duration, orientation, duration, group))
end
while length(events) > 0 && index_grad == events[1][1]
(_, event) = popfirst!(events)
push!(components, event)
end
end
return components
end
make_generic(other_block::BaseBuildingBlock) = BuildingBlock(duration(other_block), [other_block...])
Base.keys(bb::BuildingBlock) = 1:length(bb.parts)
Base.getindex(bb::BuildingBlock, i::Integer) = bb.parts[i]
duration(bb::BuildingBlock) = sum(duration, waveform_sequence(bb); init=0.)
function get_pulse(bb::BuildingBlock)
pulses = [p for p in events(bb) if p isa RFPulseComponent]
if length(pulses) == 0
error("BuildingBlock does not contain any pulses.")
end
if length(pulses) == 1
return pulses[1]
end
error("BuildingBlock contains more than one pulse. Not sure which one to return.")
end
function get_readout(bb::BuildingBlock)
readouts = [r for r in events(bb) if r isa ReadoutComponent]
if length(readouts) == 0
error("BuildingBlock does not contain any readouts.")
end
if length(readouts) == 1
return readouts[1]
end
error("BuildingBlock contains more than one readout. Not sure which one to return.")
end
end
\ No newline at end of file
module WaitBlocks
import JuMP: @constraint
import ...BuildSequences: global_model
import ...Variables: get_free_variable, VariableType, duration
import ...Components: NoGradient
import ..BaseBuildingBlocks: BaseBuildingBlock
struct Wait <: BaseBuildingBlock
duration :: VariableType
function Wait(var)
res = new(get_free_variable(var))
if !(res.duration isa Number)
@constraint global_model() res.duration >= 0
end
return res
end
end
duration(wb::Wait) = wb.duration
Base.keys(::Wait) = (Val(:empty),)
Base.getindex(wb::Wait, ::Val{:empty}) = NoGradient{1}(wb.duration)
end
\ No newline at end of file
module AllSequences
include("base_sequences.jl")
include("sequences.jl")
include("slice_select_rephases.jl")
include("epi_readouts.jl")
import .BaseSequences: BaseSequence, nrepeat
import .Sequences: Sequence
import .SliceSelectRephases: SliceSelectRephase
import .EPIReadouts: EPIReadout
end
\ No newline at end of file
module Sequences
import StaticArrays: SVector
import JuMP: @constraint
import ...Variables: get_free_variable, TR, VariableType, duration
import ...AllBuildingBlocks: Wait, BuildingBlock
import ...BuildSequences: global_model
import ...Components: EventComponent
import ..BaseSequences: ContainerBlock, BaseSequence, nrepeat, get_index_single_TR
"""
Sequence(blocks; TR=:min, nrepeat=0)
Sequence(blocks...; TR=:min, nrepeat=0)
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.
"""
struct Sequence{N} <: BaseSequence{N}
blocks :: SVector{N, <:ContainerBlock}
TR :: VariableType
end
function Sequence(blocks::AbstractVector; TR=:min)
blocks = to_block.(blocks)
actual_duration = sum(duration, blocks; init=0.)
if TR == :min
TR = actual_duration
end
res = Sequence{length(blocks)}(SVector{length(blocks)}(blocks), get_free_variable(TR))
if !(res.TR isa Number) || !(duration(res) isa Number)
@constraint global_model() res.TR >= actual_duration
end
return res
end
Sequence(blocks...; kwargs...) = Sequence([blocks...]; kwargs...)
get_index_single_TR(s::Sequence, i::Integer) = s.blocks[i]
nrepeat(::Sequence) = 0
"""
to_block(block_like)
Converst object into something that can be included in the sequence:
- :min/:max/number/variable/nothing => [`Wait`](@ref).
- `building_block` or `sequence` => no change.
- RF pulse/readout => will be embedded within a [`BuildingBlock`](@ref).
"""
to_block(cb::ContainerBlock) = cb
to_block(s::Symbol) = to_block(Val(s))
to_block(s::Union{VariableType, Nothing, Val{:min}, Val{:max}}) = Wait(s)
to_block(ec::EventComponent) = BuildingBlock([], [(0, ec)]; duration=duration(ec))
end
\ No newline at end of file
module ContainerBlocks
import ..Variables: AbstractBlock, duration, effective_time
module Abstract
import ...Variables: AbstractBlock, duration, effective_time
"""
Parent type for `BuildingBlock` or `BaseSequence`, i.e., any building block that contains other MRI components/blocks.
......
module Alternatives
import JuMP: @constraint
import ..ContainerBlocks: ContainerBlock
import ..BuildSequences: global_model, fixed
import ..Variables: duration, make_generic
import ..Abstract: ContainerBlock
import ...BuildSequences: global_model, fixed
import ...Variables: duration, make_generic
"""
AlternativeBlocks(name, blocks)
......
module BaseBuildingBlocks
import ...ContainerBlocks: ContainerBlock
import ...Components: BaseComponent, GradientWaveform, EventComponent, NoGradient, ChangingGradient, ConstantGradient, split_gradient
"""
Defines [`BaseBuildingBlock`](@ref), [`BuildingBlock`](@ref) and [`Wait`](@ref).
"""
module BuildingBlocks
import LinearAlgebra: norm
import JuMP: @constraint
import ..Abstract: ContainerBlock
import ...Components: BaseComponent, GradientWaveform, EventComponent, NoGradient, ChangingGradient, ConstantGradient, split_gradient, DelayedEvent, RFPulseComponent, ReadoutComponent
import ...Variables: qval, bmat_gradient, effective_time
import ...Variables: VariableType, duration, make_generic, get_pulse, get_readout, scanner_constraints!
"""
Basic BuildingBlock, which can consist of a gradient waveforms with any number of RF pulses/readouts overlaid
......@@ -175,4 +181,106 @@ function bmat_gradient(bb::BaseBuildingBlock, qstart, index1, index2)
end
bmat_gradient(bb::BaseBuildingBlock, qstart) = bmat_gradient(bb, qstart, nothing, nothing)
"""
BuildingBlock(waveform, events; duration=nothing, orientation=nothing, group)
Generic [`BaseBuildingBlock`](@ref) that can capture any overlapping gradients, RF pulses, and/or readouts.
The gradients cannot contain any free variables.
## Arguments
- `waveform`: Sequence of 2-element tuples with (time, (Gx, Gy, Gz)). If `orientation` is set then the tuple is expected to look like (time, G). This cannot contain any free variables.
- `events`: Sequence of 2-element tuples with (index, pulse/readout). The start time of the pulse/readout at the start of the gradient waveform element with index `index` (use [`DelayedEvent`](@ref) to make this earlier or later).
- `duration`: duration of this `BuildingBlock`. If not set then it will be assumed to be the time of the last element in `waveform`.
- `orientation`: orientation of the gradients in the waveform. If not set, then the full gradient vector should be given explicitly.
- `group`: group of the gradient waveform
"""
struct BuildingBlock <: BaseBuildingBlock
parts :: Vector{<:BaseComponent}
function BuildingBlock(parts::AbstractVector{<:BaseComponent})
res = new(duration, parts)
for (_, part) in waveform_sequence(parts)
scanner_constraints!(part)
end
return res
end
end
function BuildingBlock(waveform::AbstractVector, events::AbstractVector; duration=nothing, orientation=nothing, group=nothing)
events = Any[events...]
waveform = Any[waveform...]
ndim = isnothing(orientation) ? 1 : 3
zero_grad = isnothing(orientation) ? zeros(3) : 0.
if length(waveform) == 0 || waveform[1][1] > 0.
pushfirst!(waveform, (0., zero_grad))
events = [(i+1, e) for (i, e) in events]
end
if isnothing(duration)
duration = waveform[end][1]
end
if !(duration waveform[end][1])
@assert duration > waveform[end][1]
push!(waveform, (duration, zero_grad))
end
components = BaseComponent[]
for (index_grad, ((prev_time, prev_grad), (time, grad))) in enumerate(zip(waveform[1:end-1], waveform[2:end]))
duration = time - prev_time
if norm(prev_grad) <= 1e-12 && norm(grad) <= 1e-12
push!(components, NoGradient{ndim}(duration))
elseif norm(prev_grad) norm(grad)
push!(components, ConstantGradient(prev_grad, orientation, duration, group))
else
push!(components, ChangingGradient(prev_grad, (grad .- prev_grad) ./ duration, orientation, duration, group))
end
while length(events) > 0 && index_grad == events[1][1]
(_, event) = popfirst!(events)
push!(components, event)
end
end
return components
end
make_generic(other_block::BaseBuildingBlock) = BuildingBlock(duration(other_block), [other_block...])
Base.keys(bb::BuildingBlock) = 1:length(bb.parts)
Base.getindex(bb::BuildingBlock, i::Integer) = bb.parts[i]
duration(bb::BuildingBlock) = sum(duration, waveform_sequence(bb); init=0.)
function get_pulse(bb::BuildingBlock)
pulses = [p for p in events(bb) if p isa RFPulseComponent]
if length(pulses) == 0
error("BuildingBlock does not contain any pulses.")
end
if length(pulses) == 1
return pulses[1]
end
error("BuildingBlock contains more than one pulse. Not sure which one to return.")
end
function get_readout(bb::BuildingBlock)
readouts = [r for r in events(bb) if r isa ReadoutComponent]
if length(readouts) == 0
error("BuildingBlock does not contain any readouts.")
end
if length(readouts) == 1
return readouts[1]
end
error("BuildingBlock contains more than one readout. Not sure which one to return.")
end
struct Wait <: BaseBuildingBlock
duration :: VariableType
function Wait(var)
res = new(get_free_variable(var))
if !(res.duration isa Number)
@constraint global_model() res.duration >= 0
end
return res
end
end
duration(wb::Wait) = wb.duration
Base.keys(::Wait) = (Val(:empty),)
Base.getindex(wb::Wait, ::Val{:empty}) = NoGradient{1}(wb.duration)
end
\ No newline at end of file
module Containers
include("abstract.jl")
include("building_blocks.jl")
include("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 .Alternatives: AlternativeBlocks, match_blocks!
end
\ No newline at end of file
module BaseSequences
import ...ContainerBlocks: ContainerBlock, start_time
import ...AllBuildingBlocks: BaseBuildingBlock
import ...Variables: TR, duration
"""
Defines [`BaseSequence`](@ref) and [`Sequence`](@ref)
"""
module Sequences
import StaticArrays: SVector
import JuMP: @constraint
import ...Variables: get_free_variable, TR, VariableType, duration
import ...BuildSequences: global_model
import ...Components: EventComponent
import ..Abstract: ContainerBlock, start_time
import ..BuildingBlocks: Wait, BuildingBlock, BaseBuildingBlock
"""
Super-type of any sequence of non-overlapping building blocks that should be played after each other.
......@@ -77,5 +84,50 @@ TR(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)
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.
"""
struct Sequence{N} <: BaseSequence{N}
blocks :: SVector{N, <:ContainerBlock}
TR :: VariableType
end
function Sequence(blocks::AbstractVector; TR=:min)
blocks = to_block.(blocks)
actual_duration = sum(duration, blocks; init=0.)
if TR == :min
TR = actual_duration
end
res = Sequence{length(blocks)}(SVector{length(blocks)}(blocks), get_free_variable(TR))
if !(res.TR isa Number) || !(duration(res) isa Number)
@constraint global_model() res.TR >= actual_duration
end
return res
end
Sequence(blocks...; kwargs...) = Sequence([blocks...]; kwargs...)
get_index_single_TR(s::Sequence, i::Integer) = s.blocks[i]
nrepeat(::Sequence) = 0
"""
to_block(block_like)
Converst object into something that can be included in the sequence:
- :min/:max/number/variable/nothing => [`Wait`](@ref).
- `building_block` or `sequence` => no change.
- RF pulse/readout => will be embedded within a [`BuildingBlock`](@ref).
"""
to_block(cb::ContainerBlock) = cb
to_block(s::Symbol) = to_block(Val(s))
to_block(s::Union{VariableType, Nothing, Val{:min}, Val{:max}}) = Wait(s)
to_block(ec::EventComponent) = BuildingBlock([], [(0, ec)]; duration=duration(ec))
end
module EPIReadouts
import ...AllBuildingBlocks: LineReadout, Trapezoid, opposite_kspace_lines
import ...Containers: BaseSequence, get_index_single_TR
import ..Trapezoids: Trapezoid, opposite_kspace_lines, LineReadout
import ...Components: ADC
import ...Variables: get_free_variable, VariableType, qval, qvec, set_simple_constraints!, resolution, inverse_voxel_size, inverse_fov, resolution, get_readout, apply_simple_constraint!
import ..BaseSequences: BaseSequence, get_index_single_TR
"""
EPIReadout(resolution; ky_lines=-resolution[2]:resolution[2], recenter=false, group=:FOV, variables...)
......
module HelperFunctions
import JuMP: @constraint
import ..AllBuildingBlocks: BuildingBlock, Trapezoid, SpoiltSliceSelect, opposite_kspace_lines, SliceSelect
import ..BuildSequences: global_model, build_sequence
import ..AllSequences: Sequence, SliceSelectRephase
import ..Components: SincPulse, ConstantPulse, InstantPulse
import ..Variables: qvec, flat_time, rise_time
import ..Alternatives: AlternativeBlocks, match_blocks!
import ...Containers: AlternativeBlocks, match_blocks!, BuildingBlock
import ..Trapezoids: Trapezoid, opposite_kspace_lines, SliceSelect
import ..SpoiltSliceSelects: SpoiltSliceSelect
import ..SliceSelectRephases: SliceSelectRephase
import ...BuildSequences: global_model, build_sequence
import ...Containers: Sequence
import ...Components: SincPulse, ConstantPulse, InstantPulse
import ...Variables: qvec, flat_time, rise_time
function _get_pulse(shape, flip_angle, phase, frequency, Nzeros, group, bandwidth, duration)
......
module AllBuildingBlocks
include("base_building_blocks.jl")
include("building_blocks.jl")
module Parts
include("trapezoids.jl")
include("slice_select_rephases.jl")
include("spoilt_slice_selects.jl")
include("wait_blocks.jl")
include("epi_readouts.jl")
include("helper_functions.jl")
import .BaseBuildingBlocks: waveform, waveform_sequence, events, BaseBuildingBlock
import .BuildingBlocks: BuildingBlock
import .Trapezoids: Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines
import .SpoiltSliceSelects: SpoiltSliceSelect
import .WaitBlocks: Wait
import .SliceSelectRephases: SliceSelectRephase
import .EPIReadouts: EPIReadout
import .HelperFunctions: excitation_pulse, refocus_pulse, epi_readout, single_line_readout
end
\ No newline at end of file
module SliceSelectRephases
import ..BaseSequences: BaseSequence, get_index_single_TR
import ...AllBuildingBlocks: SliceSelect, Trapezoid
import ...Containers: BaseSequence, get_index_single_TR
import ..Trapezoids: SliceSelect, Trapezoid
import ...Variables: get_pulse, qval, apply_simple_constraint!
import ...Components: RFPulseComponent
......
......@@ -6,7 +6,7 @@ import JuMP: @constraint, @objective, objective_function
import ...BuildSequences: global_model, global_scanner
import ...Variables: VariableType, duration, rise_time, flat_time, effective_time, qval, gradient_strength, slew_rate, inverse_slice_thickness, get_free_variable, get_pulse, set_simple_constraints!, gradient_orientation
import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent
import ..BaseBuildingBlocks: BaseBuildingBlock
import ...Containers: BaseBuildingBlock
"""
......
......@@ -10,7 +10,7 @@ import ...Variables: qval, rise_time, flat_time, slew_rate, gradient_strength, v
import ...Variables: Variables, all_variables_symbols, dwell_time, inverse_fov, inverse_voxel_size, fov, voxel_size, get_gradient, get_pulse, get_readout, gradient_orientation, ramp_overlap
import ...BuildSequences: global_model
import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent, ADC
import ..BaseBuildingBlocks: BaseBuildingBlock
import ...Containers: BaseBuildingBlock
"""
......
......@@ -2,11 +2,8 @@ module Pathways
import LinearAlgebra: norm, tr
import StaticArrays: SVector, SMatrix
import ..Components: NoGradient, RFPulseComponent, ReadoutComponent, InstantGradient, GradientWaveform, DelayedEvent
import ..AllSequences: BaseSequence, Sequence
import ..AllBuildingBlocks: BaseBuildingBlock, waveform, events, waveform_sequence
import ..Containers: BaseSequence, Sequence, BaseBuildingBlock, waveform, events, waveform_sequence, start_time, AlternativeBlocks
import ..Variables: qvec, qval, bmat_gradient, VariableType, effective_time, duration, TR
import ..ContainerBlocks: start_time
import ..Alternatives: AlternativeBlocks
"""
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment