diff --git a/src/containers/linearise.jl b/src/containers/linearise.jl deleted file mode 100644 index 3b06f6184980f487660689fab717099223141471..0000000000000000000000000000000000000000 --- a/src/containers/linearise.jl +++ /dev/null @@ -1,196 +0,0 @@ -module Linearise -import StaticArrays: SVector -import ...Components: GradientWaveform, split_timestep, InstantPulse, InstantGradient3D, InstantGradient1D, edge_times -import ...Variables: amplitude, phase, gradient_strength3, duration, qval -import ..Abstract: start_time, end_time, ContainerBlock, iter_instant_gradients, iter_instant_pulses -import ..BaseSequences: BaseSequence, Sequence -import ..BuildingBlocks: BaseBuildingBlock - - -struct LinearPart{T} - start_value :: T - end_value :: T -end - -Base.iszero(lp::LinearPart) = iszero(lp.start_value) && iszero(lp.end_value) -Base.iszero(lp::LinearPart{<:AbstractVector}) = all(iszero.(lp.start_value)) && all(iszero.(lp.end_value)) -(lp::LinearPart)(time::Number) = @. (1 - time) * lp.start_value + time * lp.end_value - -""" - SequencePart(sequence, time1, time2) - -Represents the time between `time1` and `time2` of a larger [`LinearSequence`](@ref) - -The gradient, RF amplitude, and RF phase are all be modeled as changing linearly during this time. -""" -struct SequencePart - gradient :: LinearPart{SVector{3, Float64}} - gradient_origin :: SVector{3, Float64} - rf_amplitude :: LinearPart{Float64} - rf_phase :: LinearPart{Float64} - duration :: Float64 -end - - -function SequencePart(sequence::BaseSequence{N}, time1::Number, time2::Number) where {N} - tmean = (time1 + time2) / 2 - nTR = div(tmean, duration(sequence), RoundDown) - time1 -= nTR * duration(sequence) - time2 -= nTR * duration(sequence) - tmean -= nTR * duration(sequence) - if -1e-9 < time1 < 0. - time1 = 0. - end - if time2 > duration(sequence) && time2 ≈ duration(sequence) - time2 = duration(sequence) - end - if !(0 <= time1 <= time2 <= duration(sequence)) - error("Initial time ($time1) and final time ($time2) are not part of the same TR.") - end - for key in 1:N - if (end_time(sequence, key) > tmean) - 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 -1e-9 < time1 < 0. - time1 = 0. - end - if time2 > duration(bb) && time2 ≈ duration(bb) - time2 = duration(bb) - end - if !(0 <= time1 <= time2 <= duration(bb)) - error("Sequence timings are out of bound") - end - - tmean = (time1 + time2) / 2 - - 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 SequencePart( - fit_linear_part(gradient_strength3), - zero(SVector{3, Float64}), - fit_linear_part(amplitude), - fit_linear_part(phase), - time2 - time1 - ) -end - -duration(sp::SequencePart) = sp.duration - -""" - split_times(sequence(s); precision=0.01, max_timestep=Inf) - -Suggests at what times to split one or more sequence into linear parts (see [`linearise`](@ref)). - -The split times will include any time when (for any of the provided sequences): -- 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)). -""" -split_times(sequence::BaseSequence, args...; kwargs...) = split_times([sequence], args...; kwargs...) -split_times(sequences::AbstractVector{<:BaseSequence}; kwargs...) = split_times(sequences, 0., maximum(duration.(sequences)); kwargs...) - -function split_times(sequences::AbstractVector{<:BaseSequence}, tstart::Number, tfinal::Number; precision=0.01, max_timestep=Inf) - edges = Float64.([tstart, tfinal]) - for sequence in sequences - raw_edges = edge_times(sequence) - nTR_start = Int(div(tstart, duration(sequence), RoundDown)) - nTR_final = Int(div(tfinal, duration(sequence), RoundUp)) - for nTR in nTR_start:nTR_final - for time in raw_edges .+ (nTR * duration(sequence)) - if tstart < time < tfinal - push!(edges, time) - end - end - end - end - edges = sort(unique(edges)) - splits = Float64[] - - for (t1, t2) in zip(edges[1:end-1], edges[2:end]) - tmean = (t1 + t2) / 2 - - timestep = Float64(max_timestep) - for sequence in sequences - (block_time, block) = sequence(mod(tmean, duration(sequence))) - for key in keys(block) - if !(start_time(block, key) < block_time < end_time(block, key)) - continue - end - new_timestep = split_timestep(block[key], precision) - if new_timestep < timestep - timestep = new_timestep - end - end - end - nsteps = isinf(timestep) ? 1 : Int(div(t2 - t1, timestep, RoundUp)) - append!(splits, range(t1, t2, length=nsteps+1)) - end - return unique(splits) -end - - -""" - LinearSequence(sequence, times) - LinearSequence(sequence; max_timestep=Inf, precision=0.01) - LinearSequence(sequence, time1, time2; max_timestep=Inf, precision=0.01) - -A piece-wise linear approximation of a sequence. - -By default it represents the sequence between `time1=0` and `time2=TR` - -The gradient, RF amplitude, and RF phase are all be modeled as changing linearly during this time. - -If multiple sequences are provided a vector of `LinearSequence` objects are returned, all of which are split at the same time. -""" -struct LinearSequence - finite_parts :: Vector{SequencePart} - instant_pulses :: Vector{Union{Nothing, InstantPulse}} - instant_gradient :: Vector{Union{Nothing, InstantGradient3D}} -end - -LinearSequence(sequence::BaseSequence; kwargs...) = LinearSequence(sequence, 0., duration(sequence); kwargs...) -LinearSequence(container::Union{BaseSequence, AbstractVector{<:BaseSequence}}, tstart::Number, tfinal::Number; kwargs...) = LinearSequence(container, split_times(container, tstart, tfinal; kwargs...)) -LinearSequence(containers::AbstractVector{<:BaseSequence}, times::AbstractVector{<:Number}) = map(c -> LinearSequence(c, times), containers) -function LinearSequence(container::BaseSequence, times::AbstractVector{<:Number}) - parts = [SequencePart(container, t1, t2) for (t1, t2) in zip(times[1:end-1], times[2:end])] - - tstart = times[1] - tfinal = times[end] - - pulses = Union{Nothing, InstantPulse}[nothing for _ in 1:length(times)] - gradients = Union{Nothing, InstantGradient3D}[nothing for _ in 1:length(times)] - for nTR in div(tstart, duration(container), RoundDown):div(tfinal, duration(container), RoundUp) - for (to_store, func) in [ - (pulses, iter_instant_pulses), - (gradients, iter_instant_gradients), - ] - for (time, pulse) in func(container) - real_time = time + nTR * duration(container) - if !(tstart <= real_time < tfinal) - continue - end - if pulse isa InstantGradient1D - pulse = InstantGradient3D(qvec(pulse) .* pulse.orientation, pulse.group) - end - index = findmin(t -> abs(t - real_time), times)[2] - to_store[index] = pulse - end - end - end - return LinearSequence(parts, pulses, gradients) -end - - -end \ No newline at end of file