From e7c6b33b85db8aea12e31b0fafac8a61ae28f8b3 Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
Date: Wed, 6 Mar 2024 17:08:37 +0000
Subject: [PATCH] Split `split_times` from `linearise`
---
src/MRIBuilder.jl | 4 +-
src/containers/containers.jl | 2 +-
src/containers/linearise.jl | 115 ++++++++++++++++++++++++-----------
src/pathways.jl | 4 --
src/variables.jl | 3 +-
5 files changed, 84 insertions(+), 44 deletions(-)
diff --git a/src/MRIBuilder.jl b/src/MRIBuilder.jl
index 32d4305..0621e5c 100644
--- a/src/MRIBuilder.jl
+++ b/src/MRIBuilder.jl
@@ -26,8 +26,8 @@ export variables, duration, effective_time, flip_angle, amplitude, phase, freque
import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
-import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, linearise
-export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, linearise
+import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, SequencePart, split_times, split_timestep, linearise
+export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, SequencePart, split_times, split_timestep, linearise
import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway
export Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway
diff --git a/src/containers/containers.jl b/src/containers/containers.jl
index 5164e74..b31fdb8 100644
--- a/src/containers/containers.jl
+++ b/src/containers/containers.jl
@@ -9,6 +9,6 @@ import .Abstract: ContainerBlock, start_time, end_time, effective_time, readout_
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 .Linearise: linearise
+import .Linearise: SequencePart, split_times, split_timestep, linearise
end
\ No newline at end of file
diff --git a/src/containers/linearise.jl b/src/containers/linearise.jl
index 86ac78c..6a65cdf 100644
--- a/src/containers/linearise.jl
+++ b/src/containers/linearise.jl
@@ -2,7 +2,7 @@ module Linearise
import StaticArrays: SVector
import ...Components: GradientWaveform, split_timestep
import ...Variables: amplitude, phase, gradient_strength3, duration
-import ..Abstract: edge_times, start_time, end_time
+import ..Abstract: edge_times, start_time, end_time, ContainerBlock
import ..BaseSequences: BaseSequence, Sequence
import ..BuildingBlocks: BaseBuildingBlock
@@ -14,7 +14,13 @@ end
"""
-Represents a small part of a larger [`Sequence`](@ref), where the gradient, RF amplitude, and RF phase can all be modeled as changing linearly.
+ SequencePart(sequence, time1, time2)
+
+Represents the time between `time1` and `time2` of a larger [`Sequence`](@ref)
+
+The gradient, RF amplitude, and RF phase are all be modeled as changing linearly during this time.
+
+See [`linearise`](@ref) to split a sequence into such linear parts.
"""
struct SequencePart
gradient :: LinearPart{SVector{3, Float64}}
@@ -24,22 +30,67 @@ struct SequencePart
duration :: Float64
end
-duration(sp::SequencePart) = sp.duration
+function SequencePart(sequence::BaseSequence{N}, time1::Number, time2::Number) where {N}
+ if !(0 <= time1 <= time2 <= duration(sequence))
+ error("Sequence timings are out of bound")
+ end
+ for key in 1:N
+ if (end_time(sequence, key) > time1)
+ return SequencePart(sequence[key], time1 - start_time(sequence, key), time2 - start_time(sequence, key))
+ end
+ end
+end
+
+function SequencePart(bb::BaseBuildingBlock, time1::Number, time2::Number)
+ if !(0 <= time1 <= time2 <= duration(bb))
+ error("Sequence timings are out of bound")
+ end
+
+ tmean = (time1 + time2) / 2
-function linearise(seq::BaseSequence; kwargs...)
- parts = SequencePart[]
- for block in seq
- append!(parts, linearise(block; kwargs...))
+ function fit_linear_part(fn)
+ start = fn(bb, time1)
+ final = fn(bb, time2)
+ middle = fn(bb, tmean)
+ mid_deviation = @. middle - (start + final) / 2
+ return LinearPart(start .+ mid_deviation ./ 3, final .+ mid_deviation ./ 3)
end
- return parts
+
+ return SequencePart(
+ fit_linear_part(gradient_strength3),
+ zero(SVector{3, Float64}),
+ fit_linear_part(amplitude),
+ fit_linear_part(phase),
+ time2 - time1
+ )
end
-function linearise(bb::BaseBuildingBlock; precision=0.01)
- parts = SequencePart[]
- if duration(bb) <= 0.
- return parts
+duration(sp::SequencePart) = sp.duration
+
+"""
+ split_times(sequence; precision=0.01)
+
+Suggests at what times to split a sequence into linear parts (see [`linearise`](@ref)).
+
+The split times will include any time when:
+- the starting/end points of building blocks, gradients, RF pulses, or ADC readouts.
+- a gradient or RF pulse is discontinuous in the first derivative
+- the time of any instantaneous gradients, RF pulses, or readouts.
+
+Continuous gradient waveforms or RF pulses might be split up further to ensure the linear approximations meet the required `precision` (see [`split_timestep`](@ref)).
+"""
+function split_times(sequence::BaseSequence; precision=0.01)
+ splits = Float64[]
+ for block in sequence
+ append!(splits, split_times(block; precision=precision))
end
+ return sort(unique(splits))
+end
+
+
+function split_times(bb::BaseBuildingBlock; precision=0.01)
+ splits = Float64[]
edges = edge_times(bb)
for (t1, t2) in zip(edges[1:end-1], edges[2:end])
tmean = (t1 + t2) / 2
@@ -50,36 +101,28 @@ function linearise(bb::BaseBuildingBlock; precision=0.01)
if !(start_time(bb, key) < tmean < end_time(bb, key))
continue
end
- new_timestep = split_timestep(block, precision)
+ new_timestep = split_timestep(bb[key], precision)
if new_timestep < timestep
timestep = new_timestep
end
end
- nsteps = div(t2 - t1, timestep, RoundUp)
- tsplits = range(t1, t2, length=nsteps)
-
- for (t1b, t2b) in zip(tsplits[1:end-1], tsplits[2:end])
- tmeanb = (t1b + t2b) / 2
-
- function fit_linear_part(fn)
- start = fn(bb, t1b)
- final = fn(bb, t2b)
- middle = fn(bb, tmeanb)
- mid_deviation = @. middle - (start + final) / 2
- return LinearPart(start .+ mid_deviation ./ 3, final .+ mid_deviation ./ 3)
- end
-
- push!(parts, SequencePart(
- fit_linear_part(gradient_strength3),
- zero(SVector{3, Float64}),
- fit_linear_part(amplitude),
- fit_linear_part(phase),
- t2b - t1b
- ))
- end
+ nsteps = (isinf(timestep) ? 1 : Int(div(t2 - t1, timestep, RoundUp))) + 1
+ append!(splits, range(t1, t2, length=nsteps))
end
- return parts
+ return sort(unique(splits))
end
+
+"""
+ linearise(sequence[, times])
+
+Splits any [`BaseSequence`](@ref) or [`BaseBuildingBlock`](@ref) into a series of [`SequencePart`](@ref) objects where the gradients/pulses are approximated to be linear.
+
+If the `times` are not explicitly set they will be obtained from [`split_times`](@ref).
+"""
+linearise(container::ContainerBlock) = linearise(container, split_times(container))
+linearise(container::ContainerBlock, times::AbstractVector{<:Number}) = [SequencePart(container, t1, t2) for (t1, t2) in zip(times[1:end-1], times[2:end])]
+
+
end
\ No newline at end of file
diff --git a/src/pathways.jl b/src/pathways.jl
index 7eb0aa9..f243b12 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -311,10 +311,6 @@ function walk_pathway!(block::BaseBuildingBlock, walker::PathwayWalker, pulse_ef
current_time = block_start_time
for (index_inter, interruption) in events(block)
- if interruption isa Tuple
- interruption = interruption[2]
- end
-
# determine if action should be taken
if interruption isa RFPulseComponent
if iszero(length(pulse_effects))
diff --git a/src/variables.jl b/src/variables.jl
index 9103625..0200fc1 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -13,6 +13,7 @@ In addition this defines:
"""
module Variables
import JuMP: @constraint, @variable, Model, @objective, objective_function, AbstractJuMPScalar
+import StaticArrays: SVector
import ..Scanners: gradient_strength, slew_rate, Scanner
import ..BuildSequences: global_model, global_scanner, fixed
@@ -290,7 +291,7 @@ for base_fn in [:qval, :gradient_strength, :slew_rate]
@eval function $fn3(bb::AbstractBlock, args...; kwargs...)
value = $base_fn(bb, args...; kwargs...)
if value isa Number && iszero(value)
- return value
+ return zero(SVector{3, Float64})
elseif value isa AbstractVector
return value
else
--
GitLab