Segmentation of a sequence for the scanner

This Tour converts shipped pypulseqpp sequences into the scanner representation and reads what the conversion reduces them to: the repetition of each subsequence, the virtual segments it is cut into, and how those counts scale with the length of the scan.

Prerequisites: lessons 1 and 3 of the course.

A scanner interpreter prepares each virtual segment once and plays the scan as an execution stream of segment instances, so the number of virtual segments, not the number of blocks, sets what it prepares. The model is described in Scanner representation (PulSeg).

Repetition of a spin echo

The sequence is pypulseqpp’s shipped 2D spin echo at its shortest TR, written to a file and converted under the scanner limits. convert() writes the IR cache beside the sequence file; summary() reports the segmentation.

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pypulseqpp as pp
from figure_style import FAINT, INK, PAGE_WIDTH, SERIES
from pypulseqpp.sequences import se2D_sequence

from pulserver import ir

system = pp.Opts(max_grad=40, grad_unit="mT/m", max_slew=150, slew_unit="T/m/s")
work = Path(tempfile.mkdtemp())

seq = se2D_sequence(n_x=64, n_y=64, tr=None)
seq.write(work / "se2d.seq")
cache = ir.convert(work / "se2d.seq", system)
report = ir.summary(work / "se2d.seq", system)

unit = report["subsequences"][0]
print(f"{seq.num_blocks} blocks in the file; cache {cache.name}")
print(
    f"repetition: {unit['tr_size']} blocks, TR {unit['tr_duration_us'] / 1e3:.2f} ms, "
    f"played {unit['num_trs']} times"
)
for index, segment in enumerate(report["segments"]):
    first = segment["start_block"] + 1
    last = segment["start_block"] + segment["num_blocks"]
    kind = "pure delay" if segment["pure_delay"] else "events"
    print(
        f"virtual segment {index}: blocks {first}-{last}, "
        f"{segment['duration_us'] / 1e3:.2f} ms, {kind}"
    )

576 blocks in the file; cache se2d.pseg
repetition: 9 blocks, TR 18.02 ms, played 64 times
virtual segment 0: blocks 1-2, 3.56 ms, events
virtual segment 1: blocks 3-3, 0.02 ms, pure delay
virtual segment 2: blocks 4-8, 12.04 ms, events

The repetition is cut only at block boundaries where every gradient is at rest, so blocks joined by a gradient that does not rest stay in one virtual segment. A pure delay has no events: every pure delay of the repetition plays the one pure-delay virtual segment, with its own duration, as the virtual segment of each block, read from the execution stream, shows.

owner = ir.play(work / "se2d.seq")["segment"][: unit["tr_size"]]
print("virtual segment of each block of the repetition:", owner.tolist())

virtual segment of each block of the repetition: [0, 0, 1, 2, 2, 2, 2, 2, 1]

Scan length

The same conversion over the shipped 2D gradient echo, with a growing phase-encoding matrix and slice count. The number of blocks and the number of repetitions grow with the scan; the number of virtual segments does not, because every repetition plays the same base blocks with a different phase-encoding amplitude.

from pypulseqpp.sequences import gre2D_sequence

rows = []
for n_y, n_slices in ((64, 1), (128, 1), (256, 1), (256, 4)):
    path = work / f"gre_{n_y}_{n_slices}.seq"
    gre = gre2D_sequence(n_x=64, n_y=n_y, n_slices=n_slices, tr=None)
    gre.write(path)
    result = ir.summary(path, system)
    rows.append((n_y, n_slices, gre.num_blocks, result))

 ny  slices   blocks  repetitions  virtual segments
 64       1      480           80                 2
128       1      864          144                 2
256       1     1632          272                 2
256       4     6528         1088                 2

Subsequences

pypulseqpp’s echo planar sequence returns a reference prescan, one volume with the phase-encoding direction reversed, and the imaging sequence, as a chain. pypulseqpp.sequences.write() writes them as separate files linked by the NextSequence definition, and the conversion reads the chain as the subsequences of one scan.

from pypulseqpp import sequences
from pypulseqpp.sequences.sequence.epi2D_sequence import epi2d

files = sequences.write(work / "epi.seq", epi2d(system, n_x=64, n_y=64), offline=False)
print([Path(f).name for f in ir.chain(files[0])])

report = ir.summary(files[0], system)
for index, subsequence in enumerate(report["subsequences"]):
    print(
        f"subsequence {index}: {subsequence['tr_size']} blocks per repetition, "
        f"{subsequence['num_trs']} repetitions"
    )
print(f"virtual segments over the chain: {report['num_segments']}")

['epi.seq', 'epi_epi_2d.seq']
subsequence 0: 72 blocks per repetition, 3 repetitions
subsequence 1: 72 blocks per repetition, 3 repetitions
virtual segments over the chain: 2

Virtual segments are deduplicated across subsequences as well as within one. The reference prescan reverses the sign of the phase-encoding gradients, which is an amplitude of each segment instance rather than part of a base block, so both subsequences play the virtual segments printed above.