diff --git a/CHANGELOG.md b/CHANGELOG.md index 832fb18e14346886f769e1ee668a46e51d7b4324..02239c641fa6ea6837ff3e9d21b89f56466dbb4b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Finite RF pulses (`RFPulse`) to more realistically model the effect of these pulses. - Hierchical properties: MRI relaxation parameters and collision parameters (e.g., permeability) can now set at both the global level and for each obstruction. -- `run_tests.sh` function which will run the tests and produce a coverage that can be visualised in VS Code by the Coverage Gutters plugin +- Sequences can now be read from the pulseq format (http://pulseq.github.io/) using `read_pulseq`. +- `run_tests.sh` function which will run the tests and produce a coverage that can be visualised in VS Code by the Coverage Gutters plugin. ### Changed - Made units consistent with angles (i.e., phases and flip angles) in degrees, off-resonance fields in kHz, and gradients in kHz/um. - `time` function has been renamed `get_time`, so as not to conflict with `Base.time` diff --git a/Project.toml b/Project.toml index e59e3e72084660b4a9e5b42ef2b16a5d4f97b1bb..7f4456b0df32f7682e6fc242bd501a56a0855ef6 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.4.0" Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/src/sequence.md b/docs/src/sequence.md index f29795f7add334c0be63cd5f12939d2d846074af..eaf80525e9e751910ece3991bf9e51568abcad12 100644 --- a/docs/src/sequence.md +++ b/docs/src/sequence.md @@ -1,4 +1,8 @@ # [MR sequences](@id sequence) +## Pulseseq sequences +The recommended way to define sequences is to use pulseq sequence format (http://pulseq.github.io/). +Pulseq sequences can be generated using matlab (http://pulseq.github.io/) or python (https://pypulseq.readthedocs.io/en/master/). +They can be loaded into MCMRSimulator using [`read_pulseq`](@ref). ## Built-in MR sequences ### Diffusion-weighted MRI A pulsed-gradient spin-echo can be created using [`dwi`](@ref) diff --git a/src/MCMRSimulator.jl b/src/MCMRSimulator.jl index bd582d8317df73cef446a09ec42e03e20b5a6718..b62d53e05d7b98f9da251cdd008e55017f716eff 100644 --- a/src/MCMRSimulator.jl +++ b/src/MCMRSimulator.jl @@ -25,6 +25,7 @@ import Accessors: @set import Roots import PlyIO import Optim +import DataStructures: OrderedDict include("constants.jl") @@ -40,7 +41,7 @@ include("plot/plot.jl") export Sequence, InstantRFPulse, Readout, InstantGradient export MRGradients, gradient, rotate_bvec -export dwi +export dwi, read_pulseseq export Scanner, B0, max_gradient, max_slew_rate export Siemens_Connectom, Siemens_Prisma, Siemens_Terra export longitudinal, transverse, phase, vector, Spin, Snapshot, SpinOrientation, isinside, off_resonance diff --git a/src/sequence/gradients.jl b/src/sequence/gradients.jl index 91e92606a6aa84f91d59f4dcb45798a2740c7851..96ac028dc210782bd8b706a8eb7c417a64caccb0 100644 --- a/src/sequence/gradients.jl +++ b/src/sequence/gradients.jl @@ -49,11 +49,11 @@ function MRGradients(times::AbstractVector{<:Number}, amplitudes::AbstractVector end add_TR(g::MRGradients, TR::Number) = MRGradients(add_TR(g.Gx, Tr),add_TR(g.Gy, TR), add_TR(g.Gz, TR), g.origin) -control_points(g::MRGradients) = [ +control_points(g::MRGradients) = sort(unique([ control_points(g.Gx)..., control_points(g.Gy)..., control_points(g.Gz)..., -] +])) """ diff --git a/src/sequence/pulseseq.jl b/src/sequence/pulseseq.jl new file mode 100644 index 0000000000000000000000000000000000000000..c389f64ee1d05f43bca57f039e05aa94106a1d4c --- /dev/null +++ b/src/sequence/pulseseq.jl @@ -0,0 +1,273 @@ +struct PulseSeqSection + title :: String + content :: Vector{String} +end + +struct CompressedPulseSeqShape + id :: Int + num :: Int + samples :: Vector{Float} +end + + +""" + Shape(shape::CompressedPulseSeqShape, step_size) + +Create a `Shape` object based on the pulseseq Shape format presuming each sample is `step_size` milliseconds apart. +""" +function Shape(pulseseq::CompressedPulseSeqShape) + compressed = length(pulseseq.samples) != pulseseq.num + if !compressed + times = range(0, 1, length=pulseseq.num) + return Shape(times, amplitudes) + end + times = [zero(Float)] + amplitudes = [pulseseq.samples[1]] + repeating = false + time_norm = 1 / (pulseseq.num - 1) + prev_sample = pulseseq.samples[1] + prev_applied = true + for sample in pulseseq.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 + else + if !prev_applied + push!(times, times[end] + time_norm) + push!(amplitudes, amplitudes[end] + prev_sample) + end + prev_sample = sample + prev_applied = false + end + end + if !prev_applied + push!(times, times[end] + time_norm) + push!(amplitudes, amplitudes[end] + prev_sample) + end + return Shape(times, amplitudes) +end + +""" + read_pulseseq(filename; 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_pulseseq(filename; kwargs...) + keywords = open(read_pulseseq_sections, filename) + build_sequence(; kwargs..., keywords...) +end + +function read_pulseseq_sections(io::IO) + sections = PulseSeqSection[] + 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, PulseSeqSection(title, String[])) + elseif length(sections) > 0 + push!(sections[end].content, line) + else + error("Content found in pulseseq file before first section") + end + end + Dict(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 +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 + 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" => ([:id, :amp, :mag_id, :phase_id, :delay, :freq, :phase], [Int, Float, Int, Int, Int, Float, Float]), + "Gradients" => ([:id, :amp, :shape_id, :delay], [Int, Float, Int, Int]), + "TRAP" => ([:id, :amp, :rise, :flat, :fall, :delay], [Int, Float, Int, Int, Int, Int]), + "ADC" => ([:id, :num, :dwell, :delay, :freq, :phase], [Int, Int, Float, Int, Float, Float]), +) + +function parse_pulseq_section(section::PulseSeqSection, 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 pulseseq") + elseif section.title == "SHAPES" + current_id = -1 + shapes = CompressedPulseSeqShape[] + 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, CompressedPulseSeqShape(current_id, parse(Int, line[length(text)+1:end]), Float[])) + current_id = -1 + break + end + end + if !startswith(lowercase(line), "num") + @assert current_id == -1 + push!(shapes[end].samples, parse(Float, line)) + end + end + result = Dict( + [s.id => (s.num, Shape(s)) for s in shapes]... + ) + elseif section.title in ["SIGNATURE"] + # silently ignore these sections + else + error("Unrecognised pulseq section: $(section.title)") + end + return Symbol(lowercase(section.title)) => result +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) + if isnothing(scanner) + scanner = Scanner(B0=B0) + end + if version == v"1.4.0" + # load raster times (converting seconds to milliseconds) + convert = key -> parse(Float, definitions[key]) * Float(1e3) + gradient_raster = convert("GradientRasterTime") + rf_raster = convert("RadiofrequencyRasterTime") + adc_raster = convert("AdcRasterTime") + block_duration_raster = convert("BlockDurationRasterTime") + elseif version == v"1.3.1" + gradient_raster = rf_raster = block_duration_raster = Float(1e-3) # 1 microsecond as default raster + adc_raster = Float(1e-6) # ADC dwell time is in ns by default + else + error("Can only load pulseseq files with versions v1.3.1 and v1.4.0, not $(version)") + end + + gradient_control_points = ( + # times and amplitudes of gradient control points + ([zero(Float)], [zero(Float)]), + ([zero(Float)], [zero(Float)]), + ([zero(Float)], [zero(Float)]), + ) + + block_start_time = zero(Float) + rf_pulses = RFPulse[] + readouts = Readout[] + for block in values(blocks) + block_duration = zero(Float) + if !iszero(block.rf) + proc = rf[block.rf] + start_time = block_start_time + proc.delay * 1e-3 + @assert version == v"1.3.1" || iszero(proc.time_id) + (num, mag_shape) = shapes[proc.mag_id] + block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration) + mag = ConcreteShape(start_time, start_time + num * rf_raster, proc.amp * 1e-3, 0, mag_shape) + if iszero(proc.phase_id) + phase = ConcreteShape([start_time, start_time + num * rf_raster], [0, 0]) + else + (num, phase_shape) = shapes[proc.phase_id] + block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration) + phase = ConcreteShape(start_time, start_time + num * rf_raster, proc.freq * 1e-3 * 360, rad2deg(proc.phase), phase_shape) + end + push!(rf_pulses, RFPulse(mag, phase, sqrt(proc.amp^2 + proc.freq^2))) + end + if !iszero(block.adc) + proc = adc[block.adc] + mid_readout = block_start_time + proc.delay * 1e-3 + proc.dwell * proc.num * adc_raster / 2 + push!(readouts, Readout(mid_readout)) + block_duration = max(proc.delay * 1e-3 + proc.dwell * proc.num * adc_raster, block_duration) + end + for ((times_grad, amps_grad), symbol_grad) in zip(gradient_control_points, [:gx, :gy, :gz]) + grad_id = getfield(block, symbol_grad) + if iszero(grad_id) + continue + end + if grad_id in keys(gradients) + proc = gradients[grad_id] + @assert version == v"1.3.1" || iszero(proc.time_id) + start_time = block_start_time + proc.delay * 1e-3 + (num, grad_shape) = shapes[proc.shape_id] + append!(times_grad, grad_shape.times .* (num * gradient_raster) .+ start_time) + append!(amps_grad, grad_shape.amplitudes .* (proc.amp * 1e3)) + block_duration = max(proc.delay * 1e-3 + num * gradient_raster, block_duration) + elseif grad_id in keys(trap) + proc = trap[grad_id] + @assert version == v"1.3.1" || iszero(proc.time_id) + start_time = block_start_time + proc.delay * 1e-3 + times = (cumsum([0, proc.rise, proc.flat, proc.fall]) .* 1e-3) .+ start_time + amplitudes = [0, proc.amp * 1e3, proc.amp * 1e3, 0] + append!(times_grad, times) + append!(amps_grad, amplitudes) + block_duration = max((proc.delay + 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) + duration = max(block_duration, delays[block.delay].delay * 1e-3) + else + duration = block_duration + end + else + duration = block.duration * block_duration_raster + @assert duration >= block_duration + end + block_start_time += duration + end + if isnothing(TR) + TR = get(definitions, "TotalDuration", block_start_time) + end + for (times, amplitudes) in gradient_control_points + push!(times, TR) + push!(amplitudes, 0) + end + gradients = MRGradients([ConcreteShape(t, a) for (t, a) in gradient_control_points]..., zero(PosVector)) + Sequence(; scanner=scanner, TR=TR, pulses=[rf_pulses..., readouts...], gradients=gradients) +end \ No newline at end of file diff --git a/src/sequence/radio_frequency.jl b/src/sequence/radio_frequency.jl index 3cd879cc0ca43feef8f48e6eda696d16b0fb5f7f..064561134ce7fa7c3c58425915d8fa672c25c005 100644 --- a/src/sequence/radio_frequency.jl +++ b/src/sequence/radio_frequency.jl @@ -31,7 +31,7 @@ function RFPulse(times::AbstractVector, amplitudes::AbstractVector, phases=nothi end start_time(pulse::RFPulse) = min(start_time(pulse.amplitude), start_time(pulse.phase)) -end_time(pulse::RFPulse) = min(end_time(pulse.amplitude), end_time(pulse.phase)) +end_time(pulse::RFPulse) = max(end_time(pulse.amplitude), end_time(pulse.phase)) """ amplitude(rf_pulse, t1[, t2]) diff --git a/src/sequence/sequence.jl b/src/sequence/sequence.jl index 6fc7a715a16b71ea4f44f9e917c99eff934f8d3a..13a1849ddf60fd2f336d5e3ffeb651fb5d8dc853 100644 --- a/src/sequence/sequence.jl +++ b/src/sequence/sequence.jl @@ -4,7 +4,7 @@ include("gradients.jl") include("radio_frequency.jl") """ - Sequence(;TR, gradients=nothing, pulses=nothing, B0=3., interplate_gradients=:step) + Sequence(;TR, gradients=nothing, pulses=nothing, scanner=Scanner(B0), B0=3., interplate_gradients=:step) An MR sequence represented by a series of pulses repeated with a given repetition time (`TR`). @@ -135,4 +135,5 @@ gradient(seq::Sequence, time::Number) = gradient(seq.gradient, mod(time, seq.TR) gradient(seq::Sequence, time1::Number, time2::Number) = gradient(seq.gradient, mod(time1, seq.TR), iszero(mod(time2, seq.TR)) ? seq.TR : mod(time2, seq.TR)) include("timestep.jl") -include("diffusion.jl") \ No newline at end of file +include("diffusion.jl") +include("pulseseq.jl") \ No newline at end of file diff --git a/src/sequence/shape.jl b/src/sequence/shape.jl index bbd77702776f32151db759f8ecd32b7fb2107dfe..15003cd593430a6e55080081a38a42fa3ae94cac 100644 --- a/src/sequence/shape.jl +++ b/src/sequence/shape.jl @@ -135,8 +135,9 @@ struct ConcreteShape{N} t0::Float t1::Float max_amplitude::Float + shift::Float shape::Shape{N} - ConcreteShape(t0::Number, t1::Number, max_amplitude::Number, shape::Shape{N}) where {N} = new{N}(Float(t0), Float(t1), Float(max_amplitude), shape) + ConcreteShape(t0::Number, t1::Number, max_amplitude::Number, shift::Number, shape::Shape{N}) where {N} = new{N}(Float(t0), Float(t1), Float(max_amplitude), Float(shift), shape) end """ @@ -156,7 +157,7 @@ function ConcreteShape(times::AbstractVector{<:Number}, amplitudes::AbstractVect t1 = times[end] max_amplitude = maximum(abs.(amplitudes)) ConcreteShape( - t0, t1, max_amplitude, + t0, t1, max_amplitude, 0, Shape( (times .- t0) ./ (t1 - t0), iszero(max_amplitude) ? amplitudes : (amplitudes ./ max_amplitude) @@ -172,14 +173,14 @@ Produces a ConcreteShape with only zero amplitude. ConcreteShape() = ConcreteShape(Float[], Float[]) convert_time(concrete::ConcreteShape, time::Number) = (time - concrete.t0) / (concrete.t1 - concrete.t0) -amplitude(concrete::ConcreteShape, time::Number) = amplitude(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude +amplitude(concrete::ConcreteShape, time::Number) = amplitude(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude + concrete.shift amplitude_derivative(concrete::ConcreteShape, time::Number) = amplitude_derivative(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude / (concrete.t1 - concrete.t0) -amplitude_integral(concrete::ConcreteShape, t0::Number, t1::Number) = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, t1)) * concrete.max_amplitude * (concrete.t1 - concrete.t0) -amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractVector{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time.(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) -amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractRange{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) +amplitude_integral(concrete::ConcreteShape, t0::Number, t1::Number) = (amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, t1)) * concrete.max_amplitude + concrete.shift) * (concrete.t1 - concrete.t0) +amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractVector{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time.(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) .+ (concrete.shift * (concrete.t1 - concrete.t0)) +amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractRange{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) .+ (concrete.shift * (concrete.t1 - concrete.t0)) control_points(concrete::ConcreteShape) = control_points(concrete.shape) .* (concrete.t1 - concrete.t0) .+ concrete.t0 -add_TR(concrete::ConcreteShape, TR::Number) = ConcreteShape(concrete.t0 + TR, concrete.t1 + TR, concrete.max_amplitude, concrete.shape) +add_TR(concrete::ConcreteShape, TR::Number) = ConcreteShape(concrete.t0 + TR, concrete.t1 + TR, concrete.max_amplitude, concrete.shift, concrete.shape) start_time(concrete::ConcreteShape) = concrete.t0 end_time(concrete::ConcreteShape) = concrete.t1 diff --git a/test/example_pulseseq/fiddisp.seq b/test/example_pulseseq/fiddisp.seq new file mode 100644 index 0000000000000000000000000000000000000000..bfab4c79c7acfeccf22ba8d7ee60b7cb15c38e07 --- /dev/null +++ b/test/example_pulseseq/fiddisp.seq @@ -0,0 +1,55 @@ +# Pulseq sequence file +# Created by MATLAB mr toolbox + +[VERSION] +major 1 +minor 3 +revision 1 + +[DEFINITIONS] +Name fid + +# Format of blocks: +# # D RF GX GY GZ ADC EXT +[BLOCKS] +1 0 1 0 0 0 0 0 +2 1 0 0 0 0 0 0 +3 0 0 0 0 0 1 0 + +# Format of RF events: +# id amplitude mag_id phase_id delay freq phase +# .. Hz .... .... us Hz rad +[RF] +1 2500 1 2 100 0 0 + +# Format of ADC events: +# id num dwell delay freq phase +# .. .. ns us Hz rad +[ADC] +1 1024 312500 20 0 0 + +# Format of delays: +# id delay (us) +[DELAYS] +1 5000 + +# Sequence Shapes +[SHAPES] + +shape_id 1 +num_samples 120 +1 +0 +0 +97 +-1 +0 +0 +17 + +shape_id 2 +num_samples 120 +0 +0 +118 + diff --git a/test/runtests.jl b/test/runtests.jl index 5aee1ea9500e89868a88db9dddd3ce881dbf29c5..b3b1f00d037b8d609a055e3770c925124e5b931b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,8 @@ all_tests = [ "transfer", "permeability", "radio_frequency", - "hierarchical_mri" + "hierarchical_mri", + "pulseseq", ] if length(ARGS) == 0 diff --git a/test/test_pulseseq.jl b/test/test_pulseseq.jl new file mode 100644 index 0000000000000000000000000000000000000000..823c8be95dcd6e787797d5e4e78792b275a0188a --- /dev/null +++ b/test/test_pulseseq.jl @@ -0,0 +1,38 @@ +@testset "Test pulseseq input" begin + directory = joinpath(pwd(), "example_pulseseq") + @testset "read fiddisp.seq" begin + keywords = open(mr.read_pulseseq_sections, joinpath(directory, "fiddisp.seq")) + seq = mr.build_sequence(;keywords...); + + # starting with RF pulse + @test length(seq.pulses) == 1 + pulse = seq.pulses[1] + @show pulse + @test mr.start_time(pulse) == 0.1 # start after 100 microseconds delay + @test mr.end_time(pulse) == 0.22 # pulse lasts 120 microseconds + for cshape in (pulse.amplitude, pulse.phase) + @test mr.start_time(cshape) == 0.1 + @test mr.end_time(cshape) == 0.22 + end + for (time, ampl) in [ + (0.09, 0.), + (0.11, 2.5), + (0.19, 2.5), + (0.21, 0.), + (0.23, 0.), + ] + @test iszero(mr.phase(pulse, time)) + @test mr.amplitude(pulse, time) ≈ ampl + end + + # no gradients + @test mr.control_points(seq.gradient) == [0, seq.TR] + + # Single readout + @test length(seq.readout_times) == 1 + rt = seq.readout_times[1] + @assert rt ≈ 0.22 + 5 + 0.02 + (1024 * 0.3125) / 2 + + @assert seq.TR ≈ 0.22 + 5 + 0.02 + (1024 * 0.3125) + end +end \ No newline at end of file