medical-image-reconstruction

How it is built

Architecture overview

A medical scanner does not record an image. An MRI scanner records samples of the Fourier transform of the patient (k-space) through several receiver coils, a CT scanner records line integrals of X-ray attenuation (projections) at many angles, and a PET scanner records pairs of coincident photons counted along lines. Image reconstruction is the system that turns this raw data into a picture a radiologist can read, within the time the patient is on the table or in the next room. Gadgetron, an open-source framework for the task, states the goal in its title: a framework for medical image reconstruction, with a server that receives raw data over the network and returns images (source: Gadgetron documentation).

The one idea that makes it work is that every modality reduces to a forward model y = A x + noise, where x is the image, A is a linear operator (Fourier transform, line integrals, coil weighting) and y is the measured data. A fast direct inverse exists only for the ideal case: an inverse FFT for fully sampled Cartesian MRI, filtered back-projection for dense CT. Everything else, from accelerated MRI to low-dose CT and PET, is an iterative solve for x using fast implementations of A and its adjoint, plus a prior (sparsity, smoothness, a learned network) that replaces the data that was not measured. The engineering is in making A and its adjoint fast on a GPU, in streaming data through a pipeline while the scan is still running, and in shipping the result as DICOM images to the archive.

Parts:

• Raw data and pipeline host — receives the sample stream, chains processing steps, returns images — Convolutional networks (learned steps), CUDA and GPU languages
• Transform engine — fast Fourier transforms on regular and irregular sample positions, wavelets — Fast Fourier transform, non-uniform-fft (no page yet), wavelet-transform (no page yet)
• Direct reconstruction — closed-form inverse for sampling that allows it — filtered-back-projection (no page yet), voxel-ray-traversal-dda (no page yet)
• Iterative solver — least squares, regularised and statistical problems — Conjugate gradient, LU and Cholesky factorization, primal-dual-hybrid-gradient (no page yet), ml-em-tomographic-reconstruction (no page yet)
• Image enhancement — denoising, registration, calibration — image-denoising (no page yet), image-alignment-registration (no page yet), radiometric-calibration (no page yet)
• Fixed-latency hardware path — front-end signal processing near the receiver — fpga-dataflow-pipelines (no page yet), instruction-scheduling (no page yet), Static timing analysis, critical-path-analysis (no page yet), dma-transfers (no page yet), Fixed-rate control loop, PID controller
• Output and archive — images as DICOM objects sent to the archive — dicom-protocol (no page yet)

 scanner ──► receiver chain ──► raw data (header + samples per readout)
 (RF coils)   (FPGA, filters)            │  stream, TCP or pipe
                                          ▼
                          ┌──────────── pipeline host ─────────────┐
                          │ buffer ► noise prewhitening ► regrid   │
                          │   ► coil sensitivity ► solver ► combine│
                          │   ► denoise ► scale, orient           │
                          └───────┬─────────────────┬──────────────┘
                                  │ A, Aᴴ           │ images
                                  ▼                 ▼
                       GPU / SIMD kernels      DICOM writer ► archive / console
                       (FFT, gridding,         (C-STORE)
                        projectors)

How the parts fit together

Raw data interface

Reconstruction code is only shareable when the data is. The ISMRM raw data format (ISMRMRD) was proposed as a standard for MRI raw datasets, and the repository holds both the format description and a C/C++ library (source: ISMRMRD README; Inati et al., Magnetic Resonance in Medicine, 2017). Each acquired readout carries the samples and a header describing where in k-space it was taken, which coil channels it holds and which repetition it belongs to, so the pipeline can place samples without scanner-specific knowledge. In CT and PET the equivalent record is a projection or a list of detected events with detector positions. The contract with the rest of the system: samples arrive in acquisition order, the header is sufficient to build the forward operator A, and the record has a checksum so truncated transfers are detected (dicom-protocol is the output side, see below; CRC-32C checksums is the integrity block).

Pipeline host

The host is a service that accepts a stream of raw records, runs them through a chain of processing nodes and returns images. In Gadgetron the nodes are called gadgets; a gadget is a node in the chain, written as a pure function or a channel-based gadget, and the server can run in two modes: a server mode that listens on a TCP port for data over the network, and a stream mode that reads serialised data from a local file or an anonymous pipe (source: Gadgetron documentation). The steps form a directed acyclic graph; ordering them is a Topological sort. Each scan type needs its own chain. Gadgetron distributes Docker images and documents running it in Kubernetes, so the same host runs on a scanner-side computer or in a cluster (source: Gadgetron documentation).

Transform engine

The Fourier transform is the core of MRI. For fully sampled Cartesian data, reconstruction is an inverse Fast Fourier transform per coil followed by combining the coil images. Radial or spiral trajectories sample off the grid, so the transform needs a non-uniform FFT: samples are gridded onto a regular grid with a convolution kernel, then an ordinary FFT is applied, and the kernel's footprint is corrected afterwards (non-uniform-fft, no page yet; Fessler and Sutton, Nonuniform fast Fourier transforms using min-max interpolation, IEEE Trans. Signal Processing, 2003). Wavelet transforms (wavelet-transform, no page yet) serve the sparsity priors below. FIR and IIR filters cover the signal chain: the ramp filter in CT, the decimation filters after the receiver.

Direct reconstruction

For dense CT, filtered-back-projection (no page yet) applies a ramp filter to each projection and smears it back across the image along its ray. For circular cone-beam geometries, the practical algorithm of Feldkamp, Davis and Kress extends this to three dimensions (J. Opt. Soc. Am. A, 1984). It is cheap, one pass over the data, and it is the baseline that iterative methods are judged against. Its weakness is noise: the ramp filter amplifies high frequencies, which is what low-dose protocols suffer from.

Iterative solver

Iterative reconstruction treats the image as the unknown of an optimisation problem. Three solver families cover most systems:

• Linear least squares. min ||A x − y||² with a quadratic penalty gives a symmetric positive definite system, solved by Conjugate gradient using only products with A and Aᴴ; LU and Cholesky factorization are used on small dense problems such as per-pixel calibration. This is the engine for coil-sensitivity-weighted (SENSE-type) reconstruction, where multi-coil data and measured coil sensitivities unfold an undersampled image (Pruessmann et al., SENSE: Sensitivity encoding for fast MRI, Magnetic Resonance in Medicine, 1999). GRAPPA instead fills missing k-space lines by interpolation using calibration lines acquired in the scan (Griswold et al., 2002).
• Sparsity-regularised problems. Compressed sensing MRI undersamples k-space in a way that produces incoherent aliasing and recovers the image by exploiting sparsity in a transform domain (Lustig, Donoho and Pauly, Sparse MRI, Magnetic Resonance in Medicine, 2007). The non-smooth penalty is handled by first-order methods such as primal-dual hybrid gradient (primal-dual-hybrid-gradient, no page yet). The ODL library lists PDHG and Douglas-Rachford splitting among its optimisation routines (source: ODL README).
• Statistical emission tomography. PET and SPECT counts are Poisson distributed; maximum likelihood by expectation maximisation (Shepp and Vardi, IEEE Trans. Medical Imaging, 1982) is accelerated by ordered subsets (Hudson and Larkin, 1994) (ml-em-tomographic-reconstruction, no page yet). STIR describes its emphasis as iterative image reconstruction in PET and SPECT (source: STIR README).

The operator A is either applied by a fast transform (MRI), by projectors on the fly (CT), or stored as a Sparse matrices (CSR) system matrix (small PET geometries). ASTRA exposes GPU-accelerated forward and back projection, callable from MATLAB and Python, as building blocks for new algorithms (source: ASTRA README).

Learned reconstruction

Deep-learning reconstruction replaces the hand-built prior with a network. One structure is the unrolled iteration: a fixed number of solver steps, each with a learned regulariser, trained end to end; the variational network of Hammernik et al. learns such a model for accelerated MRI (Magnetic Resonance in Medicine, 2017). The fastMRI dataset, released with benchmarks, provides raw measurements and clinical images for training and evaluation of such approaches (Zbontar et al., arXiv:1811.08839). GE AIR Recon DL and Siemens Deep Resolve are vendor products of this kind. The network is a Convolutional networks stage inside the same pipeline, and the data-consistency step (the CG or FFT solve) stays classical, which is what keeps the output tied to the measurement.

Image enhancement

After reconstruction: noise reduction (image-denoising, no page yet), registration of repeated frames (image-alignment-registration, no page yet), intensity calibration (radiometric-calibration, no page yet) and resampling (spline-interpolation, no page yet). Robust statistics gives the median and MAD used for noise-level estimation and outlier rejection.

Compute layer

The cost is dominated by A and Aᴴ. They are FFTs, gridding, and ray integrals, all data-parallel, so they run on the GPU: the vendor-neutral launch model is in GPU compute shaders and the CUDA specifics in CUDA and GPU languages. Gadgetron's install check reports CUDA support and the GPU's memory (source: Gadgetron documentation). Where a GPU is not available or the problem is small, SIMD programming serves the CPU. Iterative methods accumulate rounding error over hundreds of iterations and thousands of voxels; Floating-point error lists the countermeasures (compensated summation, accumulating in higher precision), and projector and back-projector pairs must remain exact adjoints of each other or CG-type methods diverge.

Fixed-latency hardware path

Near the receiver, data rates are high and latency must be fixed, so the first stages (down-conversion, decimation filters) run as FPGA pipelines with a constant initiation interval (fpga-dataflow-pipelines, no page yet), checked by static timing analysis (static-timing-analysis, no page yet).

The acquisition hardware also moves. In CT the gantry must rotate at a steady speed and the table advance at a set pitch, and in MRI the gradient amplifiers must deliver the commanded current, otherwise the sample positions that the raw records claim are not where the data was really taken. A Fixed-rate control loop runs each of these regulators at an absolute period, and a PID controller closes the loop on measured speed or current. Their contract with reconstruction is the header: the achieved angle or gradient waveform is recorded with each readout, so the reconstruction uses the trajectory that actually happened, not the planned one.

Output and archive

Images are written as DICOM objects and sent to the picture archive. The DICOM standard defines the service classes in Part 4 (including the C-FIND query operation), and Part 5 defines transfer syntaxes as encoding rules that are negotiated when an association is established between two application entities (source: DICOM standard, Parts 4 and 5). The reconstruction output needs patient and series identifiers copied from the raw data header, a unique image identifier per slice, and pixel scaling attributes, because the display depends on them. dicom-protocol (no page yet) is the block that would explain the message exchange.

Data and control flow

One accelerated MRI exam, scanner to archive:

1. The scanner sends a header describing the protocol, then one record per readout: samples for every coil channel plus the position in k-space.
2. The host buffers readouts per slice or volume, applies noise decorrelation across coil channels from a noise scan, and removes known corruption by robust statistics.
3. A calibration step estimates coil sensitivities from the fully sampled centre of k-space or from a separate scan.
4. The solver runs: gridding or FFT for A, CG or an unrolled network for the iterations, each iteration costing a few FFTs per coil (the dominant cost, order of tens of iterations).
5. Coil data are combined into one image per slice; magnitude, scaling and orientation are set from the header.
6. Images go to the DICOM writer, with new series identifiers, and are sent to the archive and the console; the exam is complete when the last series is stored.

One CT slab with iterative reconstruction: projections arrive and are calibrated; an FBP image is computed first as start value and fallback; iterations of forward projection, comparison, back-projection and update refine it with a noise penalty; if the iterative job misses its deadline the FBP result is stored so the study is not delayed.

Scaling

In the order you hit them:

1. Operator cost. Each iteration needs A and Aᴴ for every coil or view. Levers: GPU kernels for the FFT and projectors, fewer iterations through a good start value, and ordered subsets, which accelerate emission tomography (Hudson and Larkin).
2. Memory. 3D multi-coil or cone-beam datasets are large (order of gigabytes), and a GPU holds less than the host. Levers: process in slabs or per slice after a transform along the fully sampled axis, use half precision for storage, or chunk projections.
3. Latency. The radiologist waits for images. Levers: stream readouts into the pipeline during the scan, overlap transfer and compute, keep a cheap direct reconstruction as the first image.
4. Throughput of the site. Several scanners share a recon farm. Levers: more hosts behind a queue, container images with the same pipeline, GPU sharing by job.
5. Numerical drift. Long iterations and large sums lose accuracy. Levers: compensated summation, double-precision accumulators, matched adjoints.

Failure modes

Real implementations

The alternative: use the vendor's reconstruction. In a clinical system the vendor's cleared software is the reference, and open-source tools are for research, for new sequences, and for prototyping that the vendor interface accepts as an add-on. Build your own pipeline when the sampling pattern or the prior is the research question.

Sources

• Gadgetron: documentation; M. S. Hansen, T. S. Sørensen, Gadgetron: An open source framework for medical image reconstruction, Magn. Reson. Med., 2013. DOI
• S. J. Inati et al., ISMRM Raw data format: A proposed standard for MRI raw datasets, Magn. Reson. Med., 2017. DOI; repo
• J. A. Fessler, B. P. Sutton, Nonuniform fast Fourier transforms using min-max interpolation, IEEE Trans. Signal Process., 2003. DOI
• L. A. Feldkamp, L. C. Davis, J. W. Kress, Practical cone-beam algorithm, J. Opt. Soc. Am. A, 1984. DOI
• L. A. Shepp, Y. Vardi, Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imaging, 1982. DOI
• H. M. Hudson, R. S. Larkin, Accelerated image reconstruction using ordered subsets of projection data, IEEE Trans. Med. Imaging, 1994. DOI
• K. P. Pruessmann et al., SENSE: Sensitivity encoding for fast MRI, Magn. Reson. Med., 1999. DOI
• M. A. Griswold et al., Generalized autocalibrating partially parallel acquisitions (GRAPPA), Magn. Reson. Med., 2002. DOI
• M. Lustig, D. Donoho, J. M. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magn. Reson. Med., 2007. DOI
• K. Hammernik et al., Learning a variational network for reconstruction of accelerated MRI data, Magn. Reson. Med., 2018. DOI
• J. Zbontar et al., fastMRI: An Open Dataset and Benchmarks for Accelerated MRI, 2018. arXiv:1811.08839
• DICOM standard: Part 4, Part 5
• Toolboxes: BART, ASTRA, TIGRE, STIR, ODL, fastMRI