Reconstruction plugins#

A reconstruction plugin is a file, <name>.py in a --plugins directory of the reconstruction proxy, holding a ReconPlugin subclass and a module-level PLUGIN instance. This one is bart pics on ESPIRiT maps:

# recon/gre.py
import torch
import bartorch.tools as bt
from bartorch import apps, priors
from pulserver import recon

class Pics(recon.ReconPlugin):
    def recon(self, context, branch, data):
        kspace = torch.from_numpy(data.data.kspace)  # (coils, y, x)
        maps = bt.ecalib(kspace, maps=1)
        image = apps.pics(kspace, maps, regularizers=priors.Wavelet((-1, -2), 0.005))
        return recon.ReconResult(image.abs().numpy())

PLUGIN = Pics()

data is the ReconData of one reconstruction unit, and data.data its k-space, (coils, ..., readout), with the axes data.data.axes names. Readouts are placed by their echo along the readout and by their encoding counters along the encoded axes (Placement); a sequence sets the counters with pypulseqpp.sequences.Labels or pp.make_label as in PyPulseq, and a readout placed over another is warned about. The proxy runs the plugin in a worker process, one per series, over an MRD stream enriched from the sequence’s design (Raw-data enrichment and routing). The client names the plugin of a series in its config text (Reconstruction clients), independently of the scanner-sequence plugin the series was played from.

Reconstruction units#

A reconstruction unit is the set of readouts reconstructed together. It holds the readouts of one branch and one encoding space that share their slice, contrast, cardiac phase, repetition, set and average counters; the segment and user counters do not separate units. A counter named in axes becomes an axis of the unit’s k-space instead, so that the echoes or averages of an image are one unit.

A unit closes when the flag triggers names for its branch has arrived at every position along its axes, and recon is called with it once. A unit is released when recon returns: no reference to it remains, so a series holds only the units open at once, and a plugin that needs the data of an earlier unit keeps the parts it needs. A unit still open at the last readout of the measurement (LAST_IN_MEASUREMENT) or at the end of the stream is reconstructed then, under its own branch, in the order the units opened.

Argument	Sets
`gadgets`	`Gadget` steps run over each readout before it is placed
`triggers`	`{branch: flag}`, the `AcquisitionFlag` closing a unit of the branch, several combined with `\|`. The first branch is that of every readout `branch_for` does not route elsewhere. The default is `{"imaging": LAST_IN_MEASUREMENT}`
`axes`	the counters that are axes of a unit instead of separating units
`merge`	the counters that neither separate units nor are axes: a unit waits for its closing flag at each of their values, and their readouts fill one k-space, as the place in a train (`contrast`, `ECO`) of an echo train or an inversion shot does
`require_flags`, `reject_flags`	the flags an acquisition must all carry, and those any one of which excludes it
`buffered`	whether readouts are placed in `data.data` and `data.ref`; with `False` a plugin reads `data.acquisitions` and a unit closes at its first flag, whatever its `axes`

recon receives the unit as a ReconData:

Attribute	Holds
`data`	the k-space of the imaging readouts, a `ReconBuffer`, or `None` when the unit placed none
`ref`	the k-space of the parallel-imaging calibration readouts, laid out as `data` but holding only the lines of the phase-encode and partition axes its readouts cover, which start at `ref.origin` on the grid of `data`; `None` when the unit has none
`counters`	the unit’s image counters not in `axes`, by MRD name
`waveforms`	the waveforms received since the previous unit closed
`acquisitions`	every readout the unit received, as it arrived

A readout flagged IS_PARALLEL_CALIBRATION is placed in ref only, one flagged IS_PARALLEL_CALIBRATION_AND_IMAGING in both, and one flagged IS_PHASECORR_DATA in neither unless it is also a calibration readout.

Hooks#

The runtime calls three hooks over one stream:

Hook	Called
`startup()`	once, before any acquisition
`recon()`	for each unit, as it closes
`finish()`	once, after the last unit

Only recon has to be written. The runtime passes every acquisition to receive(), which is not overridden: it applies require_flags and reject_flags, runs the gadgets over the readout, adds it to the unit of the branch branch_for() names, and calls recon for each unit that closes. branch_for returns no branch for a noise measurement, nor for a navigator readout unless triggers declares a "navigator" branch; every other readout belongs to the first branch triggers declares.

>>> import numpy as np
>>> import pulserver.mrd as mrd
>>> import pulserver.recon as recon
>>> class DropNoise(recon.Gadget):
...     def __call__(self, acquisition, data):
...         noise = mrd.has_acquisition_flag(acquisition, "ACQ_IS_NOISE_MEASUREMENT")
...         return None if noise else data
>>> class RootSumOfSquares(recon.ReconPlugin):
...     def __init__(self):
...         super().__init__(
...             gadgets=[DropNoise()],
...             triggers={"imaging": mrd.AcquisitionFlag.LAST_IN_SLICE},
...         )
...     def recon(self, context, branch, data):
...         image = np.fft.fftshift(np.fft.ifft2(data.data.kspace))
...         return recon.ReconResult(np.sqrt(np.sum(np.abs(image) ** 2, axis=0)))
>>> PLUGIN = RootSumOfSquares()

A ReconResult is packaged as an MRD image, so a plugin builds no image header. The field of view is that of the unit’s encoding space, the position, orientation and physiology time stamps are those of the unit’s reference acquisition, data.data.reference, and the time stamp is the earliest of the unit’s acquisitions. reference=<acquisition> takes the geometry from another acquisition. The MRD image holds the array as returned, float32 or complex64, with no scale applied.

dicom=True sends the result as DICOM instead, and the integer pixels are made then. The first floating-point image of a series fixes the rescale of its DICOM pixels, and a result states its own in its attributes:

return recon.ReconResult(
    image,
    dicom=True,
    attributes={"RescaleSlope": 1e-4, "RescaleIntercept": 0.0},
)

A value outside the range of the stored integers is clipped and counted; see Raw-data enrichment and routing.

Each series runs on its own copy of PLUGIN (spawn()). context.exam, an ExamCache, is shared by the series of one exam, such as for a coil calibration; under the proxy a stored value reaches the next series pickled, as a copy without its cleanup.

context.device is the GPU the proxy gave the series, such as "cuda:0", or None. A child process started with the spawn method runs functions from importable modules, not from the plugin file.

Readout placement#

A readout of a Cartesian encoding space is placed so that its echo, center_sample, lies at sample N // 2 of the readout axis, the samples from discard_pre to number_of_samples - discard_post being the ones copied. The counter of a line is placed at counter - center + extent // 2 along the encoded axis, center being the counter of the k-space centre in the header’s encoding limits, so the counters of a partial-Fourier or undersampled scan need not start at 0. A readout or a line that falls outside the buffer raises a ValueError. data.data.readout is the first and last sample placed.

A plugin brings its readouts to the reconstruction matrix with two gadgets, run in this order. AsymmetricEcho zero-fills a partial echo to the full echo that is symmetric about its center_sample, and declares the zeros as discarded so that they are not placed. RemoveReadoutOversampling crops a full echo to the readout field of view of the reconstruction, the ratio of the header’s encoded to reconstructed field of view along the readout, and needs bartorch (the coils extra). A gadget states the echo of the readout it returns by assigning center_sample, discard_pre and discard_post of the acquisition it is given, which is a copy of the one the stream delivered:

>>> import ismrmrd
>>> acquisition = ismrmrd.Acquisition()
>>> acquisition.resize(5, 2)
>>> acquisition.center_sample = 1
>>> full = recon.AsymmetricEcho()(acquisition, np.ones((2, 5), dtype=np.complex64))
>>> full.shape[-1], acquisition.discard_pre, acquisition.center_sample
(8, 3, 4)

Non-Cartesian data#

A buffer holds the trajectory of its samples beside the k-space, as the acquisitions carry it: k along the sequence’s x, y and z axes, in 1/m. grid_trajectory() returns it in the grid units and layout bartorch.linop.NUFFT takes, with the coils of the buffer’s kspace as a batch axis of the image:

import torch
from bartorch.linop import NUFFT

buffer = data.data
nufft = NUFFT(
    torch.from_numpy(buffer.grid_trajectory()),
    image_shape=(buffer.coils, *buffer.image_shape),
)
coil_images = nufft.adjoint(torch.from_numpy(buffer.kspace))

The partitions of a stack of spokes or spirals are placed by their kspace_encode_step_2 counter and carry no kz: the partition axis is transformed with an FFT before the in-plane NUFFT.

Coil sensitivities, noise and compression#

coil_maps() returns the coil sensitivities of a unit as a complex torch tensor (coils, [z,] y, x) on context.device. They are estimated by the function the plugin passes, from the unit’s calibration k-space data.ref, or taken from the maps the stream or the exam holds (Calibration, noise and coil compression states the order and the conditions for reuse). Prewhiten is a gadget and CoilCompression a step called from recon. All three need bartorch, which the coils extra installs.

# recon/gre.py
import torch
from bartorch import apps, priors
from pulserver import mrd, recon

class Pics(recon.ReconPlugin):
    def __init__(self):
        super().__init__(
            gadgets=[
                recon.Prewhiten(),
                recon.AsymmetricEcho(),
                recon.RemoveReadoutOversampling(),
            ],
            triggers={"imaging": mrd.AcquisitionFlag.LAST_IN_SLICE},
        )
        self.compression = recon.CoilCompression(8)

    def recon(self, context, branch, data):
        data = self.compression(context, data)
        maps = recon.coil_maps(context, data, estimate=apps.nlinv_maps)
        if data.data is None:  # a unit of calibration readouts only
            return None
        kspace = torch.from_numpy(data.data.kspace).to(context.device)
        image = apps.pics(kspace, maps, regularizers=priors.Wavelet((-1, -2), 0.005))
        return recon.ReconResult(image.abs().cpu().numpy())

PLUGIN = Pics()

Prewhiten consumes the noise readouts of the stream and whitens every other readout with them, or with the noise covariance a noise series left in the exam. Without either it passes readouts unchanged and logs a warning once; Prewhiten(required=True) raises instead. CoilCompression(n) keeps n virtual channels, at most the channels of the first unit, with the basis of that unit’s calibration k-space (its imaging k-space where it has none), and projects every later unit onto it. A unit is compressed before its maps are requested, so that they are estimated in the basis of the data they are used with.

A calibration unit, whose data.data is None, stores its maps for its slice, and the imaging units of that slice take them from the stream. A series that is only a calibration leaves its maps to the exam by assigning context.coil_sensitivities, and a noise series leaves its whitening with publish():

class Calibration(recon.ReconPlugin):
    def recon(self, context, branch, data):
        recon.coil_maps(context, data, estimate=apps.nlinv_maps)
        context.coil_sensitivities = context.coil_maps[data.counters["slice"]]

class Noise(recon.ReconPlugin):
    def __init__(self):
        super().__init__(gadgets=[recon.Prewhiten()])

    def recon(self, context, branch, data):
        return None

    def finish(self, context):
        self.gadget(recon.Prewhiten).publish()

The exam holds one set of maps. A later series is given them only where the coil labels, the whitening, the compression and the geometry equal its own, and otherwise MissingCalibration is raised, naming the first field that differs with both values. The failure text the client receives carries the message. required=False returns None instead of raising.

DICOM pixel values#

A DICOM dataset stores integer pixels, and the integers are made when a result with dicom=True is converted, from the image the plugin returned. The dataset’s rescale maps them to the values of the image, value = stored × RescaleSlope + RescaleIntercept, with one mapping for each DICOM series, a series being the images that share image_series_index.

An integer image is stored as it is. Its dataset carries a rescale only when the attributes of the result state one.
The first floating-point image of a series fixes the mapping of the series. The intercept is 0, so that zero is stored as zero, and the slope stores the peak magnitude of the image at half the stored range, which leaves headroom for later images of the series up to twice as large. The stored type is 16-bit, signed when that image has a negative value and unsigned when it has none. An image that states ArrayMinimum and ArrayMaximum is mapped as the array of images it belongs to, so the mapping of a series that begins with a volume does not depend on which partition arrives first. A complex image is stored as its magnitude. An image with no nonzero finite value fixes no mapping.
Every later image of the series is stored under the mapping of the series, so the ratios between the values of its images are the ratios between their stored integers.
A mapping stated in the attributes of a result, RescaleSlope or RescaleIntercept, is the mapping of that image: stored = round((value - intercept) / slope), a slope left out being 1 and an intercept left out 0. It is also the mapping of the series when the image is the first floating-point image of the series, and the images that follow it and state none are stored under it. The slope is nonzero and both are finite.
A value outside the range of the stored type is clipped to the nearest end of the range, never wrapped, and NaN is stored as 0. The clipped pixels are counted for each series, and the first clipping of a series is logged with its count. An image whose peak exceeds twice the peak its series was mapped from clips under the mapping of the series, unless the plugin states a mapping.

RescaleType is US, unspecified: the units are those of the reconstruction.

Logging#

A plugin logs with the standard logging module, under a logger of its own (logging.getLogger(__name__)). Under the proxy, what a plugin logs, prints or warns is appended to the proxy’s --logfile, each record stamped with its time, process ID, level and logger name.

Running a plugin offline#

run() reconstructs an ISMRMRD HDF5 file in the calling process. A file recorded as the scanner sends it, such as pulserver scan --mrd raw.h5 writes, is enriched from its design store first:

images = PLUGIN.run("raw.h5", store="designs")

Calling the plugin on an AcquisitionBucket runs startup, passes every acquisition to receive, reconstructs the units still open and runs finish, and returns the last output. A header only needs to describe the encoded space and the receiver channels:

>>> from types import SimpleNamespace
>>> matrix = SimpleNamespace(matrixSize=SimpleNamespace(x=8, y=4, z=1))
>>> header = SimpleNamespace(
...     encoding=[SimpleNamespace(encodedSpace=matrix, reconSpace=matrix)],
...     acquisitionSystemInformation=SimpleNamespace(receiverChannels=2),
... )
>>> bucket = mrd.AcquisitionBucket.from_arrays(
...     np.ones((4, 2, 8), dtype=np.complex64),
...     labels={"kspace_encode_step_1": np.arange(4)},
... )
>>> result = PLUGIN(bucket, recon.ReconContext.offline(header))
>>> result.data.shape
(4, 8)

Six shipped plugins are complete reconstructions, returning the values of their transform or solve unscaled. They are searched after every reconstruction plugin directory, so a client can name one, such as nufft, without a file of its own:

Plugin	Reconstructs
`nufft`	The same images as `pics`, for a radial, spiral, PROPELLER or zero-echo-time trajectory or a stack of them, from readouts whitened with the stream’s noise measurement (`Prewhiten`); each solved as `bart pics -t -R W` over `bart nlinv -t` sensitivities of the unit’s samples near the k-space centre. The readouts of every segment of a unit, the blades of a PROPELLER scan or the shells of a zero-echo-time one, enter one solve; a trajectory that encodes the partition direction without a partition axis is solved as one volume. Needs the `coils` extra
`nufft_train`	`nufft`’s images from readouts that number their place in a train of excitations (`contrast`, `ECO`), merged so that the readouts of every place fill one image. Needs the `coils` extra
`pics`	One image per slice, contrast, cardiac phase, set and repetition, averages summed, from readouts placed by their encoding counters, whitened with the stream’s noise measurement, completed to full echoes and cropped to the reconstruction field of view as they arrive (`Prewhiten`, `AsymmetricEcho`, `RemoveReadoutOversampling`); each image solved as `bart pics -R W` over the coil sensitivities of `coil_maps()` (`bart nlinv` of the unit’s calibration readouts, else the maps stored for its slice or the exam’s, else a fit to the unit’s own low-resolution centre), then completed by `bart homodyne` along a partial-Fourier axis. A unit of calibration readouts only makes no image and stores its maps. Needs the `coils` extra
`pics_train`	`pics`’s images from readouts that number their place in an echo train or an inversion shot (`contrast`, `ECO`), merged so that the readouts of every place fill one image. Needs the `coils` extra
`epi`	`pics`’s images from EPI readouts, each whitened, resampled off its ramps onto the matrix and corrected for its odd/even phase against the shot’s navigator as it arrives (`bart` ramp operator, `estimate_epi_phase`). The phase-encode-reversed reference set is an image of its own, and with PyHySCO installed (GPL-3.0, not a dependency) each later image of its slice is corrected for susceptibility distortion against it. Needs the `coils` extra
`pmc`	`pics_train`’s images, and a pose from each three-plane navigator (`NAV` or `RTFEEDBACK` readouts): planes gridded with Pipe-Menon density on one thread, coils combined by root sum of squares, registered against the first navigator and filtered by an extended Kalman filter. The pose is published to the scan when its design sets `EnablePmc`. Needs the `coils` extra

Every shipped plugin whitens its readouts with Prewhiten, which consumes the noise readouts of the stream and logs a warning once where the stream has none and the exam stores no noise covariance. A volume is sent as one image per partition. A plugin file that reexports one under another name, from pulserver.recon.handlers.pics import PLUGIN, is the same reconstruction.

Reconstruction plugins

Contents