Skip to content
Snippets Groups Projects

Add writing to Pulseq files

Merged Add writing to Pulseq files
pulseq-refactor into main
Merged Michiel Cottaar requested to merge pulseq-refactor into main
Viewing commit f095a94f
Prev Next
Show latest version
2 files
+ 10
− 5
Compare changes
  • Side-by-side
  • Inline
Files
2
src/sequence_io/pulseq.jl
+ 143
− 278
"""
Module converting MRIBuilder sequences to and from sequences recognised by [`PulseqIO`](@ref).
"""
module Pulseq
import ...Variables: duration
import Interpolations: linear_interpolation
import ..PulseqIO.Types: PulseqSequence, PulseqBlock, PulseqTrapezoid, PulseqGradient, PulseqRFPulse, PulseqShape, PulseqADC, PulseqExtension
import ...Containers: Sequence, BuildingBlock, BaseBuildingBlock, events, waveform, iter_blocks
import ...Components: GenericPulse, ADC, RFPulseComponent
import ...Scanners: Scanner
import ...Components: GenericPulse, GradientWaveform, ADC
import ...Containers: BuildingBlock, Sequence, Wait
import DataStructures: OrderedDict
import Interpolations: linear_interpolation, Flat
import StaticArrays: SVector
struct PulseqSection
title :: String
content :: Vector{String}
end
import ...Variables: duration, nsamples, dwell_time, make_generic
struct CompressedPulseqShape
id :: Int
num :: Int
samples :: Vector{Float64}
function Sequence(pulseq::PulseqSequence; scanner=nothing, B0=nothing)
if isnothing(scanner)
use_B0 = isnothing(B0) ? get(pulseq.definitions, :B0, 3.) : B0
scanner = Scanner(B0=use_B0)
end
blocks = BuildingBlock.(pulseq.blocks; pulseq.definitions..., version=pulseq.version)
return Sequence(blocks; name=Symbol(get(pulseq.definitions, :Name, "from_pulseq")), scanner=scanner)
end
"""
control_points(shape::CompressedPulseqShape)
function BuildingBlock(pulseq::PulseqBlock; version, BlockDurationRaster, RadiofrequencyRasterTime, GradientRasterTime, kwargs...)
stated_duration = pulseq.duration * BlockDurationRaster * 1e3
Returns a tuple with:
- vector of times
- vector of amplitudes
"""
function control_points(pulseq::CompressedPulseqShape)
compressed = length(pulseq.samples) != pulseq.num
if !compressed
times = range(0, 1, length=pulseq.num)
return (times, pulseq.samples)
end
times = [zero(Float64)]
amplitudes = [pulseq.samples[1]]
repeating = false
time_norm = 1 / (pulseq.num - 1)
prev_sample = pulseq.samples[1]
prev_applied = true
for sample in pulseq.samples[2:end]
if repeating
nrepeats = Int(sample) + 2
if prev_applied
nrepeats -= 1
end
push!(times, times[end] + nrepeats * time_norm)
push!(amplitudes, amplitudes[end] + nrepeats * prev_sample)
repeating = false
prev_applied = true
elseif sample == prev_sample
repeating = true
events = []
if !isnothing(pulseq.rf)
f(samples) = isnothing(pulseq.rf.time) ? [samples[1], samples..., samples[end]] : samples
if isnothing(pulseq.rf.time)
time = [0., ((1:length(pulseq.rf.magnitude.samples)) .- 0.5)..., length(pulseq.rf.magnitude.samples)] .* RadiofrequencyRasterTime .* 1e3
else
if !prev_applied
push!(times, times[end] + time_norm)
push!(amplitudes, amplitudes[end] + prev_sample)
end
prev_sample = sample
prev_applied = false
time = pulseq.rf.time.samples .* 1e3 * RadiofrequencyRasterTime
end
push!(events, (
pulseq.rf.delay * 1e-3,
GenericPulse(
time,
f(pulseq.rf.magnitude.samples) * pulseq.rf.amplitude * 1e-3,
rad2deg.(f(pulseq.rf.phase.samples) .+ pulseq.rf.phase_offset .+ pulseq.rf.frequency .* time .* 2π)
)
))
end
if !prev_applied
push!(times, times[end] + time_norm)
push!(amplitudes, amplitudes[end] + prev_sample)
if !isnothing(pulseq.adc)
dwell_time = pulseq.adc.dwell * 1e-6
push!(events, (
pulseq.adc.delay * 1e-3,
ADC(
pulseq.adc.num,
dwell_time,
dwell_time * pulseq.adc.num / 2,
1.
)
))
end
return (times, amplitudes)
end
"""
read_pulseq(IO; scanner=Scanner(B0=B0), B0=3., TR=<sequence duration>)
Reads a sequence from a pulseq file (http://pulseq.github.io/).
Pulseq files can be produced using matlab (http://pulseq.github.io/) or python (https://pypulseq.readthedocs.io/en/master/).
"""
function read_pulseq(io::IO; kwargs...)
keywords = read_pulseq_sections(io)
return build_sequence(; kwargs..., keywords...)
end
grads = [pulseq.gx, pulseq.gy, pulseq.gz]
min_duration = max(
maximum(e[1] + duration(e[2]) for e in events; init=0.),
maximum(vcat(_control_times.(grads, GradientRasterTime)...); init=0.)
)
function read_pulseq_sections(io::IO)
sections = PulseqSection[]
version = nothing
title = ""
for line in readlines(io)
line = strip(line)
if length(line) == 0 || line[1] == '#'
continue # ignore comments
end
if line[1] == '[' && line[end] == ']'
if title == "VERSION"
version = parse_pulseq_section(sections[end]).second
end
# new section starts
title = line[2:end-1]
push!(sections, PulseqSection(title, String[]))
elseif length(sections) > 0
push!(sections[end].content, line)
if min_duration > stated_duration
if version == v"1.3.1"
stated_duration = min_duration
else
error("Content found in pulseq file before first section")
error("Minimum duration to play all RF/gradient/ADC events exceeds stated duration.")
end
end
Dict(filter(pair -> !isnothing(pair), parse_pulseq_section.(sections, version))...)
end
function parse_pulseq_ordered_dict(strings::Vector{<:AbstractString}, names, dtypes)
as_dict = OrderedDict{Int, NamedTuple{Tuple(names), Tuple{dtypes...}}}()
for line in strings
parts = split(line)
@assert length(parts) == length(names)
values = parse.(dtypes, split(line))
@assert names[1] == :id
as_dict[values[1]] = (; zip(names, values)...)
end
return as_dict
times = sort(unique(vcat([0., stated_duration], _control_times.(grads, GradientRasterTime)...)))
waveform = [(t, _get_amplitude.(grads, t, GradientRasterTime)) for t in times]
return BuildingBlock(waveform, events)
end
function parse_pulseq_properties(strings::Vector{<:AbstractString})
result = Dict{String, Any}()
for s in strings
(name, value) = split(s, limit=2)
result[name] = value
_control_times(::Nothing, ::Number) = Float64[]
_control_times(trap::PulseqTrapezoid, ::Number) = cumsum([trap.delay, trap.rise, trap.flat, trap.fall]) * 1e-3
function _control_times(grad::PulseqGradient, raster::Number)
if isnothing(grad.time)
return ((1:length(grad.shape.samples)) .- 0.5) .* (1e3 * raster)
else
return grad.time.samples .* (1e3 * raster)
end
return result
end
section_headers = Dict(
"BLOCKS" => ([:id, :duration, :rf, :gx, :gy, :gz, :adc, :ext], fill(Int, 8)),
("BLOCKS", v"1.3.1") => ([:id, :delay, :rf, :gx, :gy, :gz, :adc, :ext], fill(Int, 8)),
("DELAYS", v"1.3.1") => ([:id, :delay], [Int, Int]),
("RF", v"1.3.1") => ([:id, :amp, :mag_id, :phase_id, :delay, :freq, :phase], [Int, Float64, Int, Int, Int, Float64, Float64]),
"RF" => ([:id, :amp, :mag_id, :phase_id, :time_id, :delay, :freq, :phase], [Int, Float64, Int, Int, Int, Int, Float64, Float64]),
("GRADIENTS", v"1.3.1") => ([:id, :amp, :shape_id, :delay], [Int, Float64, Int, Int]),
"GRADIENTS" => ([:id, :amp, :shape_id, :time_id, :delay], [Int, Float64, Int, Int, Int]),
"TRAP" => ([:id, :amp, :rise, :flat, :fall, :delay], [Int, Float64, Int, Int, Int, Int]),
"ADC" => ([:id, :num, :dwell, :delay, :freq, :phase], [Int, Int, Float64, Int, Float64, Float64]),
)
function parse_pulseq_section(section::PulseqSection, version=nothing)
if section.title == "VERSION"
props = parse_pulseq_properties(section.content)
result = VersionNumber(
parse(Int, props["major"]),
parse(Int, props["minor"]),
parse(Int, props["revision"]),
)
elseif section.title == "DEFINITIONS"
result = parse_pulseq_properties(section.content)
elseif (section.title, version) in keys(section_headers)
result = parse_pulseq_ordered_dict(section.content, section_headers[(section.title, version)]...)
elseif section.title in keys(section_headers)
result = parse_pulseq_ordered_dict(section.content, section_headers[section.title]...)
elseif section.title == "EXTENSION"
println("Ignoring all extensions in pulseq")
elseif section.title == "SHAPES"
current_id = -1
shapes = CompressedPulseqShape[]
for line in section.content
if length(line) > 8 && lowercase(line[1:8]) == "shape_id"
current_id = parse(Int, line[9:end])
continue
end
for text in ("num_uncompressed", "num_samples")
if startswith(lowercase(line), text)
@assert current_id != -1
push!(shapes, CompressedPulseqShape(current_id, parse(Int, line[length(text)+1:end]), Float64[]))
current_id = -1
break
end
end
if !startswith(lowercase(line), "num")
@assert current_id == -1
push!(shapes[end].samples, parse(Float64, line))
end
end
result = Dict(
[s.id => (s.num, control_points(s)...) for s in shapes]...
)
elseif section.title in ["SIGNATURE"]
# silently ignore these sections
return nothing
_get_amplitude(::Nothing, ::Number, ::Number) = 0.
function _get_amplitude(trap::PulseqTrapezoid, time::Number, raster::Number)
amp = trap.amplitude * 1e-3
edges = _control_times(trap, raster)
if !(edges[1] < time < edges[end])
return 0.
elseif time < edges[2]
return amp * (time - edges[1]) / (1e-3 * trap.rise)
elseif time < edges[3]
return amp
else
error("Unrecognised pulseq section: $(section.title)")
return amp * (edges[end] - time) / (1e-3 * trap.fall)
end
return Symbol(lowercase(section.title)) => result
end
function align_in_time(pairs...)
interps = [linear_interpolation(time, ampl; extrapolation_bc=Flat()) for (time, ampl) in pairs]
all_times = sort(unique(vcat([time for (time, _) in pairs]...)))
return (all_times, [interp.(all_times) for interp in interps]...)
function _get_amplitude(grad::PulseqGradient, time::Number, raster::Number)
amp = grad.amplitude * 1e-3
edges = _control_times(grad, raster)
return amp * linear_interpolation(edges, grad.shape.samples, extrapolation_bc=0.)(time)
end
function build_sequence(; scanner=nothing, B0=3., TR=nothing, definitions, version, blocks, rf=nothing, gradients=nothing, trap=nothing, delays=nothing, shapes=nothing, adc=nothing, name=:Pulseq)
if isnothing(scanner)
scanner = Scanner(B0=B0)
end
if version == v"1.4.0"
# load raster times (converting seconds to milliseconds)
convert = key -> parse(Float64, definitions[key]) * Float64(1e3)
gradient_raster = convert("GradientRasterTime")
rf_raster = convert("RadiofrequencyRasterTime")
adc_raster = convert("AdcRasterTime")
block_duration_raster = convert("BlockDurationRaster")
elseif version == v"1.3.1"
gradient_raster = rf_raster = block_duration_raster = Float64(1e-3) # 1 microsecond as default raster
adc_raster = Float64(1e-6) # ADC dwell time is in ns by default
else
error("Can only load pulseq files with versions v1.3.1 and v1.4.0, not $(version)")
end
full_blocks = BuildingBlock[]
function PulseqSequence(seq::Sequence{S}) where {S}
definitions = (
Name=S,
AdcRasterTime=1e-9,
BlockDurationRaster=1e-9,
RadiofrequencyRasterTime=1e-9,
GradientRasterTime=1e-9,
TotalDuration=duration(seq) * 1e-3,
B0=seq.scanner.B0,
)
blocks = [PulseqBlock(block; definitions...) for (_, block) in iter_blocks(seq)]
return PulseqSequence(
v"1.4.0",
definitions,
blocks
)
end
for block in values(blocks)
block_duration = 0.
events = []
if !iszero(block.rf)
proc = rf[block.rf]
(num, times_shape, amplitudes_shape) = shapes[proc.mag_id]
block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration)
ampl_times = times_shape .* (num * rf_raster)
ampl_size = amplitudes_shape .* (proc.amp * 1e-3)
if iszero(proc.phase_id)
phase_times = ampl_times
phase_size = ampl_times .* 0
else
(num, times_shape, phase_shape) = shapes[proc.phase_id]
block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration)
phase_times = times_shape .* (num * rf_raster)
phase_size = rad2deg.(phase_shape) .+ phase_times .* (proc.freq * 1e-3 * 360) .+ rad2deg(proc.phase)
function PulseqBlock(block::BaseBuildingBlock; BlockDurationRaster, AdcRasterTime, kwargs...)
rf = nothing
adc = nothing
for (delay, event) in events(block)
if event isa RFPulseComponent
if !isnothing(rf)
error("Pulseq does not support a single building block containing multiple RF pulses.")
end
if version != v"1.3.1" && !iszero(proc.time_id)
(num, _, time_shape) = shapes[proc.time_id]
times = time_shape .* rf_raster
ampl = ampl_size
phase = phase_size
@assert length(times) == length(ampl) == length(phase)
else
(times, ampl, phase) = align_in_time((ampl_times, ampl_size), (phase_times, phase_size))
gen = make_generic(event)
rf = PulseqRFPulse(
maximum(gen.amplitude) * 1e3,
PulseqShape(gen.amplitude ./ maximum(gen.amplitude)),
PulseqShape(deg2rad.(gen.phase)),
PulseqShape(gen.time .* 1e-3),
Int(div(delay, 1e-3, RoundNearest)),
0.,
0.
)
elseif event isa ADC
if !isnothing(rf)
error("Pulseq does not support a single building block containing multiple ADC events.")
end
push!(events, (proc.delay * 1e-3, GenericPulse(times, ampl, phase)))
end
if !iszero(block.adc)
proc = adc[block.adc]
push!(events, (proc.delay * 1e-3, ADC(proc.num, proc.dwell * 1e-6, proc.dwell * proc.num * 1e-6 / 2, 1.)))
block_duration = max(proc.delay * 1e-3 + proc.dwell * proc.num * 1e-6, block_duration)
adc = PulseqADC(
nsamples(event),
div(dwell_time(event), AdcRasterTime, RoundNearest),
Int(div(delay, 1e-3, RoundNearest)),
0., 0.
)
else
error("Cannot write $(typeof(event)) to Pulseq.")
end
grad_shapes = []
for symbol_grad in [:gx, :gy, :gz]
grad_id = getfield(block, symbol_grad)
if iszero(grad_id)
push!(grad_shapes, (Float64[], Float64[]))
continue
end
if !isnothing(gradients) && grad_id in keys(gradients)
proc = gradients[grad_id]
start_time = proc.delay * 1e-3
(num, grad_time, grad_shape) = shapes[proc.shape_id]
end
if version != v"1.3.1" && !iszero(proc.time_id)
(num, _, time_shape) = shapes[proc.time_id]
times = time_shape .* gradient_raster .+ start_time
@assert length(times) == length(grad_shape)
else
times = grad_time .* (num * gradient_raster) .+ start_time
end
push!(grad_shapes, (times, grad_shape .* (proc.amp * 1e-9)))
block_duration = max(start_time + num * gradient_raster, block_duration)
elseif !isnothing(trap) && grad_id in keys(trap)
proc = trap[grad_id]
start_time = proc.delay * 1e-3
times = (cumsum([0, proc.rise, proc.flat, proc.fall]) .* 1e-3) .+ start_time
push!(grad_shapes, (times, [0, proc.amp * 1e-9, proc.amp * 1e-9, 0]))
block_duration = max((start_time + proc.rise + proc.flat + proc.fall) * 1e-3, block_duration)
else
error("Gradient ID $grad_id not found in either of [GRADIENTS] or [TRAP] sections")
end
end
if version == v"1.3.1"
if !iszero(block.delay)
actual_duration = max(block_duration, delays[block.delay].delay * 1e-3)
else
actual_duration = block_duration
end
grads = []
times = [t for (t, _) in waveform(block)]
for dim in 1:3
amplitudes = [g[dim] for (_, g) in waveform(block)]
if iszero(maximum(abs.(amplitudes); init=0.))
push!(grads, nothing)
else
actual_duration = block.duration * block_duration_raster
@assert actual_duration >= block_duration
push!(grads, PulseqGradient(
maximum(amplitudes) * 1e3,
PulseqShape(amplitudes ./ maximum(amplitudes)),
PulseqShape(times .* 1e-3),
0.,
))
end
function extend_grad!(pair)
(times, amplitudes) = pair
if iszero(length(times)) || !iszero(times[1])
pushfirst!(times, 0.)
pushfirst!(amplitudes, 0.)
end
if times[end] != actual_duration
push!(times, actual_duration)
push!(amplitudes, 0.)
end
end
extend_grad!.(grad_shapes)
arrs = align_in_time(grad_shapes...)
waveform = [(t, SVector{3, Float64}(gx, gy, gz)) for (t, gx, gy, gz) in zip(arrs...)]
push!(full_blocks, BuildingBlock(waveform, events))
end
full_block_duration = sum(duration.(full_blocks))
if isnothing(TR)
total_duration = parse(Float64, get(definitions, "TotalDuration", "0")) * 1000
TR = iszero(total_duration) ? full_block_duration : total_duration
end
if TR > full_block_duration
push!(full_blocks, Wait(TR - full_block_duration))
elseif TR < full_block_duration
@warn "Given TR ($TR) is shorter than the total duration of all the building blocks ($full_block_duration). Ignoring TR..."
end
Sequence(full_blocks; scanner=scanner, name=name)
return PulseqBlock(
Int(div(1e-3 * duration(block), BlockDurationRaster, RoundNearest)),
rf,
grads...,
adc,
Tuple{PulseqExtension, Int}[]
)
end
end
\ No newline at end of file