Scanning a phantom without a scanner

Scanning a phantom without a scanner#

The virtual scanner plays a stored design’s IR cache, acquires an analytic phantom along the trajectory the cache plays or simulates the cache on the phantom’s isochromats, and sends the series to a reconstruction proxy as the scanner’s reconstruction client does (The virtual scanner). A scanner-sequence plugin and a reconstruction plugin are exercised together, through the production design calls and proxy.

A browser build of the virtual scanner’s console runs without an installation at https://pulserver.github.io/MaRGE/.

Generate a design#

Generate the design as the interpreter host process does, with the prescribed field-of-view offset in the protocol block (Running the services):

printf '[Protocol]\nTE: 5000\nnx: 64\nny: 64\nfov_offset_x: 20.0\n[Protocol End]\n' \
  | pulserver design generate --plugins sequences --plugin gre2d --limits limits.txt --store designs

The reply is GENERATED <id>.

Start a proxy#

Start a reconstruction proxy on the same store (Running the services):

pulserver proxy --store designs --port 9002 --plugins recon

Acquire and send the series#

from pulserver import virtual

design = "..."  # the identifier the generation replied
sequence = f"designs/{design}/sequence.seq"
phantom = virtual.Phantom(
    [virtual.Ellipse((0.02, 0.0, 0.0), (0.08, 0.06))], coils=4
)
readouts = virtual.acquire(sequence, phantom)
received = virtual.send(
    ("127.0.0.1", 9002),
    design,
    readouts,
    position_mm=(20.0, 0.0, 0.0),
    config="pics",
)

config names the reconstruction plugin, here the shipped pics, as the config text of a reconstruction client does (Reconstruction clients); the proxy refuses a series whose config names none.

The phantom lies in the physical frame, whose axes are the logical ones under a prescription without a rotation; the object above is centred on the prescribed field of view, so it appears at the centre of the image. position_mm is written as each acquisition’s position; send writes no fov_offset_mm user parameter, so the proxy corrects the samples for the offset the design was converted at and for no further displacement. received holds the images, DICOM datasets and texts the reconstruction returned; a text beginning pulserver: reports a refused or failed series.

Scan water and fat#

An ellipse of fat carries its chemical shift, and the acquisition the magnet’s field strength, at which the shift is resolved; off_resonance_hz adds a frequency offset common to every spin:

water = virtual.Ellipse((0.02, 0.0, 0.0), (0.06, 0.05))
fat = virtual.Ellipse((0.02, 0.07, 0.0), (0.04, 0.01), shift_ppm=-3.45)
readouts = virtual.acquire(
    sequence, virtual.Phantom([water, fat], coils=4), field_t=3.0
)

Fat is acquired at its chemical shift, and a fat saturation the design plays leaves it the magnetization its pulse leaves at that frequency. Generate the design under limits whose B0 is the same field: the host resolves the ppm offsets of the design’s RF pulses at it when it builds the IR (Running the services).

Simulate relaxation and the RF pulses#

The analytic acquisition acts on the phantom’s density alone. For relaxation, flip angles and slice profiles, sample the phantom as isochromats and play the cache on them (The virtual scanner); the readouts are sent as before:

tissue = virtual.Phantom(
    [
        virtual.Ellipse((0.02, 0.0, 0.0), (0.08, 0.06), t1=0.9, t2=0.07),
        virtual.Ellipse((0.02, 0.02, 0.0), (0.02, 0.015), t1=1.8, t2=0.25),
    ],
    coils=4,
)
readouts = virtual.simulate(sequence, tissue.isochromats(1e-3))

The isochromats approximate the phantom in k-space below \(1/(2\Delta)\) for a grid spacing \(\Delta\), here 1 mm, so choose a spacing several times finer than the pixel; their number, and the time the simulation takes, grow as \(1/\Delta^2\), or as \(1/\Delta^3\) for a phantom that fills a volume. Sample such a phantom only where the design’s excitation pulses excite it, virtual.BrainWeb().isochromats(1e-3, field_t=3.0, region=virtual.excited(sequence)), giving excited the prescription’s rotation where there is one. A phantom with a chemical shift is sampled at the magnet’s field, tissue.isochromats(1e-3, field_t=3.0). The simulation starts from the magnetization the isochromats hold, so a second scan of the same isochromats continues the first: sample them again, or call their reset, to start from equilibrium.

Stream a scan in real time, with its sound#

A scan played against its clock releases each readout once it is acquired, and the client sends it on; the sound of its gradients can be written as it plays:

import wave

import numpy as np

scan = virtual.Scan(sequence, tissue.isochromats(1e-3))
with wave.open("scan.wav", "wb") as audio:
    audio.setnchannels(2)
    audio.setsampwidth(2)
    audio.setframerate(int(virtual.SAMPLE_RATE))

    def acquired():
        for chunk in scan.chunks(speed=1.0):
            print(f"{chunk.stop:6.2f} s of {scan.duration:.2f} s")
            audio.writeframes(np.round(32767 * chunk.sound.T).astype("<i2").tobytes())
            yield from chunk.readouts

    received = virtual.send(("127.0.0.1", 9002), design, acquired(), config="pics")

The spans are simulated in a thread of their own, ahead of the clock. At speed=1.0 the clock starts once the simulation, at the rate it has run so far, stays ahead of it to the end of the scan, which then lasts as long as it would on a scanner; preparing, when given, receives the time left before the clock starts. A span simulated after its end on the clock holds the clock until it is. Without speed, the spans come as fast as they are simulated.

Scan from the command line#

pulserver scan runs a scan without a script. It imports a sequence file or generates a design from a scanner-sequence plugin, as the design calls do, plays its IR cache on a phantom in the Bloch simulation, and writes the series to an ISMRMRD file, streams it to a reconstruction proxy, or both:

pulserver scan --seq sequence.seq --limits limits.txt \
  --orientation coronal --center 10 -5 3 --mrd raw.h5 --sound scan.wav

pulserver scan --plugins sequences --plugin gre2d --protocol protocol.txt \
  --limits limits.txt --store designs --recon 127.0.0.1:9002 --output images

The design call’s reply, IMPORTED <id> or GENERATED <id>, or its ERROR line, is written to standard output, and the scan clock to standard error.

  • --limits is the [Limits] block of Running the services; the phantom’s chemical shifts are resolved at its B0. --protocol holds a protocol block, as pulserver design generate reads it on standard input; without it, the plugin’s defaults are designed.

  • --orientation names the rotation from logical to physical axes by the physical direction of each logical axis:

    Orientation

    Readout

    Phase encoding

    Slice selection

    axial

    +x

    +y

    +z

    coronal

    +x

    +z

    −y

    sagittal

    +y

    +z

    +x

    --rotation gives the nine elements of another, row by row, and --center the field-of-view centre, in mm along the physical axes. They replace the field-of-view entries of the protocol block.

  • --phantom is a JSON file of ellipses, each an object of the fields of Ellipse, in its units:

    {"ellipses": [
  {"centre": [0.0, 0.0, 0.0], "semi_axes": [0.08, 0.06], "t1": 0.9, "t2": 0.07},
  {"centre": [0.02, 0.02, 0.0], "semi_axes": [0.02, 0.015], "shift_ppm": -3.45}
]}

    --phantom brainweb scans BrainWeb, BrainWeb’s normal brain, downloaded on first use by the brainweb extra (pip install 'pulserver[brainweb]'); its --spacing is a whole number of millimetres. Without --phantom, the phantom is seven vials of water around one of fat, with T1 from 0.3 s to 2.0 s and T2 from 0.04 s to 0.3 s. Its isochromats are those in the slabs the design’s excitation pulses excite (excited()). --spacing, in mm, and --coils set its isochromats and its receive coils; --coil scans it with one of the scanner’s coils instead, body, body/head48 or head8/head32, named transmit/receive, whose sensitivities bartorch samples from BART’s coil models, which the coils extra installs (pip install 'pulserver[coils]'). The sensitivities at every isochromat are written to a temporary file mapped into memory; where the temporary directory is a memory file system, set TMPDIR to one on disk.

  • --spins spreads each voxel over that many isochromats, at the quantiles of the Lorentzian line of its tissue’s T2′, at its centre or, with --voxel box, over it; --diffusion lets BrainWeb’s tissue classes diffuse (Bloch engine). A phantom file gives each ellipse its t2_prime, in s, and diffusion, in m²/s. Isochromats that diffuse are played block by block, and a voxel’s diffusion attenuation needs many spins to be resolved.

  • --nod DEGREES PERIOD turns the subject about the physical x axis through the isocentre by DEGREES times the sine of 2π t / PERIOD, and --drift X Y Z translates it along the physical axes, in mm/min, t being the time from the start of the scan (RigidMotion). Isochromats that move are played block by block.

  • --recon streams the series to a reconstruction proxy, which looks the design up in its store: give that store as --store, or the proxy’s design intake as --push. The images go to images.h5 in --output, the DICOM datasets to the files they are named by, and a text beginning pulserver: ends the command with status 1. --reconstruction names the reconstruction plugin the proxy reconstructs the series with, independently of the sequence the design was generated from. Without it, a shipped --plugin is reconstructed with the shipped reconstruction paired with it (Shipped sequences), and a design of another plugin or imported with --seq needs it.

  • --speed plays the scan that many times as fast as a scanner, once the simulation is far enough ahead, writing the time left before it is to standard error; without it, the scan runs as fast as the simulation.

  • --device cuda reads the ADC windows and carries the runs of repetitions on the CUDA device. A window outside a run under a changing gradient, such as a spiral’s, whose isochromats lie on a lattice along the axes its k moves along, as a phantom’s do, is summed onto the lattice by a Triton kernel and transformed by cuFINUFFT; every other one, under a held gradient or with its isochromats off any lattice, as a moving subject’s are, is summed sample by sample by a second Triton kernel. A run of repetitions, Cartesian, radial or of ZTE spokes, is carried and spread onto its windows’ grids by more Triton kernels, a tile of repetitions at a time. The gpu extra installs all of them (pip install 'pulserver[gpu]'). The rest of the scan is simulated on the CPU.

Prescribe an orientation#

Add the nine fov_rotation_ij entries to the protocol block the design is generated from: element (i, j) of the rotation \(R\) from logical to physical axes, the identity’s where absent (Protocol resolution). The design is checked in the physical frame \(R\) gives, where logical axes played together add on one physical axis, so derate the limits it is designed under for \(R\) in the design_max_grad and design_max_slew lines of the limits file (Running the services). The acquisition and the client are given the same \(R\), and the client sends \(R\) times the offset as the field-of-view centre:

import numpy as np

rotation = np.array([[0.866025, -0.5, 0.0], [0.5, 0.866025, 0.0], [0.0, 0.0, 1.0]])
centre = rotation @ (0.02, 0.0, 0.0)
phantom = virtual.Phantom(
    [virtual.Ellipse((0.0, 0.0, 0.0), (0.08, 0.06))],
    coils=4,
    rotation=rotation,
    position=centre,
)
readouts = virtual.acquire(sequence, phantom, rotation=rotation)  # or virtual.simulate
received = virtual.send(
    ("127.0.0.1", 9002),
    design,
    readouts,
    position_mm=1e3 * centre,
    rotation=rotation,
    config="pics",
)

The phantom, posed at the prescribed centre with its axes turned by \(R\), appears at the centre of the image as it would at the isocentre without a rotation, and the images carry the columns of \(R\) as their read, phase and slice directions.

Serve a scanner console#

pulserver console answers a scanner console over a WebSocket: the design calls on the interpreter’s text blocks, an exam on a subject, and scans of the designs it generates, played on the subject’s phantom:

pulserver console --plugins sequences --limits limits.txt --store designs \
  --recon 127.0.0.1:9002 --port 8765

--plugins and --recon-plugins may be repeated, the first directory holding a plugin supplying it. With --recon-plugins DIR in place of --recon, each scan is reconstructed in the console’s own process by the plugins of DIR, checked and enriched as the proxy checks and enriches a series (LocalReconstruction), and no proxy runs. --origin names an origin whose browser pages the console serves, and may be repeated; without it, pages from every origin are served, and a client that sends no Origin, which a browser always sends, is served either way. --fields DIR takes the coils from the field maps in DIR, <coil>.npz for each of body, head8, head32 and head48, and <coil>_vops.npz for the two that transmit, as mariepy writes them for BrainWeb’s head (The virtual scanner). Every exam is then on BrainWeb, and every design is made under the VOP limits of the exam’s transmit coil; maps solved at another frequency than the Larmor frequency of the limits’ B0 are refused.

A scan is simulated on the phantom’s isochromats in the slabs its excitation pulses excite (excited()), sampled --spacing apart, 1 mm by default. Where they number more than --max-isochromats, two million by default, as over a whole head, the spacing is coarsened in steps of 1 mm until they do not, counting --spins isochromats per voxel. --device, --spins, --voxel, --diffusion, --nod and --drift act as they do for pulserver scan; a subject’s motion starts anew with each scan.

pulserver’s image runs such a console by default, with BrainWeb’s normal brain and the field maps solved in it, --fields=/console/fields, so that every exam is on BrainWeb and downloads nothing, and with bartorch for the non-Cartesian reconstructions. It serves the pages of https://pulserver.github.io, http://localhost:8000 and http://127.0.0.1:8000 on port 8765 of this computer, and plays each scan at the scanner’s speed, --speed=1:

docker run -d --restart unless-stopped --name pulserver \
  -p 127.0.0.1:8765:8765 ghcr.io/pulserver/pulserver

Options after the image’s name are added to its console’s: an option of one value takes the value given last, and --plugins, --recon-plugins and --origin add to the image’s. A console whose scans spread each voxel over four isochromats, in tissue that diffuses, of a subject that nods through ±2° with a period of 8 s:

docker run -d --restart unless-stopped --name pulserver \
  -p 127.0.0.1:8765:8765 ghcr.io/pulserver/pulserver \
  --spins 4 --diffusion --nod 2 8

Its sequences are the ones listed under Shipped sequences, reconstructed by the shipped plugins of Reconstruction plugins unless a scan names another, and each protocol starts at values the image’s limits play.

Directories mounted at /console/user/plugins and /console/user/recon add sequences and reconstructions to those pulserver ships. They are searched first, so a file there takes the place of the shipped file of that name. A link must resolve inside the container, as a relative link within the mounted directory does. A file added or changed is used from the next call on:

docker run -d --restart unless-stopped --name pulserver \
  -p 127.0.0.1:8765:8765 \
  -v "$PWD/sequences:/console/user/plugins:ro" \
  -v "$PWD/recon:/console/user/recon:ro" \
  ghcr.io/pulserver/pulserver

Each request is a JSON object carrying call and an id that every reply to it repeats:

call

Fields

Replies

plugins

plugins: the scanner-sequence plugin names

recons

recons: the reconstruction plugin names, the shipped ones and those of --recon-plugins

coils

coils: each coil’s name and its transmit and receive channels

list, validate, generate, import

plugin, block

status and reply, as pulserver design answers; design for a generated or imported design

exam

subject, coil

localizer: the axial, coronal and sagittal images of the subject’s phantom, as base64 DICOM files

scan

design, rotation (nine elements), centre_mm, sound, recon

at a speed, preparing about twice a second until the clock starts, the wall-clock time left in s or null before there is an estimate; clock and duration after each span played, with the span’s sound when asked, as base64 of 16-bit little-endian stereo samples at rate Hz; dicom and name for each image the reconstruction returns, text, then done with the status

cancel

stops the scan in progress

The recon field of a scan names the reconstruction plugin of the scan, independently of the plugin the design was generated from. Without it, a design of a shipped sequence is reconstructed with the shipped reconstruction paired with it (Shipped sequences), and any other design, such as an imported one, is refused with an error naming the missing reconstruction. A console that reconstructs through a proxy lists the shipped reconstructions only, and uses a name the proxy holds all the same.

Console.answer answers one request in process with the same replies, without the id, for a console that reaches pulserver by another route than a WebSocket, such as the messages of a Web Worker running pulserver beside a browser page.

A subject named brainweb is BrainWeb’s normal brain, and so is every subject of a console with field maps; any other is the vials. An exam is scanned in the coil it names, or in the one it had, --coil at first; a head coil of BART’s models needs the coils extra. Scans of an exam that excite the same slabs play on one set of isochromats, each from equilibrium, so that a repeated scan neither builds them nor computes again the pulses an earlier one played. A console loads BrainWeb once, and computes the field its head adds while the first exam’s localizer is viewed. The localizer is drawn from the phantom’s proton density, so an exam can be planned before any scan. MaRGE is such a console when MARGE_PULSERVER holds the address, ws://127.0.0.1:8765 here: its sequences are then the plugins, its subject names the phantom, its RF coil names the coil, and it opens each exam on the localizer. The browser build of MaRGE in pulserver/MaRGE runs it in a browser tab, opened with ?console=ws://127.0.0.1:8765, and plays each scan’s sound as it streams; it is published at https://pulserver.github.io/MaRGE/.

See also#