Parameter maps from scanner images

In this tour you estimate a \(T_2\) map from the magnitude images a scanner exports, without the raw data: you read a multi-echo spin-echo series from DICOM, fit the decay voxel by voxel, and write the map back as a DICOM series of the same study and as a NIfTI volume. The geometry and the echo times pass from the scanner’s files through the fit to the output unchanged, so the map overlays the images it was computed from.

Besides the pixel values, the fit needs the contrast timings and the geometry, and both are stored in the DICOM headers. Many protocols export each echo of a multi-echo acquisition, or each inversion time of an inversion-recovery experiment, as a DICOM series of its own. bartorch.io.read_dicom() reads several series as one, sorts the images into contrasts by their echo, inversion and repetition times, and returns those times with the images: they are the sampling points of the signal model. The fit is that of lesson 9 without the Fourier encoding: bartorch.apps.mobafit() runs Levenberg-Marquardt on the signal equation alone, one voxel at a time.

The noise of a magnitude image is Rician where the signal is small and has a positive floor, which a decay model fitted to it reads as a longer \(T_2\). With magnitude=True, mobafit() fits the magnitude of the model to the magnitude of the data, which takes the phase out of the problem but not the floor; here the echoes stay above it.

Prerequisites: lesson 9, 9. Nonlinear inversion for T2 mapping.

Learning objectives

• Read a multi-echo series stored as one DICOM series per echo, with its echo times and its voxel-to-world affine.

• Fit bartorch.nlop.MultiEcho to magnitude images with bartorch.apps.mobafit().

• Write the map as a DICOM series of the same study, and as NIfTI, in the geometry of the input.

DICOM export

The input is a scanner export: a CPMG multi-echo spin-echo acquisition of eight echoes, 10 ms apart, stored as one DICOM series per echo, numbered 101 to 108. The slice is a simple head phantom of scalp fat, grey and white matter and two ventricles of cerebrospinal fluid, tilted by 12 degrees about the left-right axis, as a slice aligned to the AC-PC line is. The code that writes the export is hidden; scanner is its folder.

Series import

read_dicom() reads all eight series at once. It returns the images as (contrasts, slices, rows, columns) in order of echo time, the echo times themselves in milliseconds, and the affine from voxel indices to RAS millimetres. Every echo has to cover the same slices.

import torch

from bartorch import apps, io, nlop

echoes = io.read_dicom(scanner, series=range(101, 109))
te = echoes.timings["echo_time"]

print(f"images {tuple(echoes.image.shape)}, TE {te.tolist()} ms")
voxel = echoes.affine[:3, :3].norm(dim=0)
print("voxel size", " x ".join(f"{v:.1f}" for v in voxel), "mm")

images (8, 1, 128, 128), TE [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0] ms
voxel size 1.7 x 1.7 x 4.0 mm

Three of the echoes, on one window. The cerebrospinal fluid in the ventricles, with the longest \(T_2\), keeps its signal across the train, while the scalp fat loses most of it.

Decay fitting

Build bartorch.nlop.MultiEcho on the echo times the files state and on the in-plane shape of the slice. The fit is voxel by voxel, so you can drop the slice axis; for a volume, build the model on (slices, rows, columns) instead. DICOM intensities are in arbitrary units: the fit scales the echoes to unit peak and returns the amplitude in the units of the files. The map is set to zero where the first echo is below 5 % of the peak, outside the head.

images = echoes.image[:, 0].to(torch.complex64)
M = nlop.MultiEcho(te.tolist(), tuple(images.shape[1:]))
maps = apps.mobafit(images, M, magnitude=True, T2=80.0)

head = images[0].abs() > 0.05 * images.abs().max()
t2 = torch.where(head, maps["T2"], torch.zeros(()))

for name, (region, truth, _) in compartments.items():
    median = float(t2[torch.from_numpy(region)].median())
    print(f"{name:>13}: T2 {median:6.1f} ms, simulated {truth:5.1f} ms")

         fat: T2   45.0 ms, simulated  45.0 ms
 grey matter: T2   85.0 ms, simulated  85.0 ms
white matter: T2   70.0 ms, simulated  70.0 ms
         CSF: T2  250.1 ms, simulated 250.0 ms

DICOM output

You write the map back to DICOM as a series of the same study: passing the first echo’s dataset as header lends the new series its patient, study and frame of reference, and you give it a number and a description of its own. ImageType marks it as derived. NIfTI stores the same affine, and that is what a registration or segmentation package reads.

import tempfile

from pathlib import Path

output = Path(tempfile.mkdtemp())
io.write_dicom(
    output / "t2-map",
    t2,
    echoes.affine,
    header=echoes.header,
    SeriesNumber=201,
    SeriesDescription="T2 map",
    ImageType=["DERIVED", "PRIMARY", "T2 MAP"],
)
io.write_nifti(output / "t2-map.nii.gz", t2, echoes.affine)

PosixPath('/tmp/tmpzq1oi1vm/t2-map.nii.gz')

Reading the series back confirms the study and the geometry:

written = io.read_dicom(output / "t2-map")
study = written.header.StudyInstanceUID == echoes.header.StudyInstanceUID
shift = float((written.affine - echoes.affine).abs().max())
print(f"patient {written.header.PatientID}, same study {study}")
print(f"largest affine difference {shift:.1e} mm")

• 
• 

patient P0001, same study True
largest affine difference 2.0e-06 mm

The map you read back from DICOM, beside the \(T_2\) the images were simulated with, and the difference between the two. The error is largest in the cerebrospinal fluid, where eight echoes over 80 ms sample only the start of a 250 ms decay.

As a spec

What this tour built, stated the way you would ask an agent for it:

Estimate a T2 map from a multi-echo spin-echo exported as one DICOM
series per echo (series 101-108), with bartorch. Read all series with
io.read_dicom, taking the echo times from echoes.timings["echo_time"]
and the geometry from echoes.affine. Fit nlop.MultiEcho on those echo
times with apps.mobafit(magnitude=True, T2=80.0), and zero the map
where the first echo is below 5 % of the peak. Write the map with
io.write_dicom as series 201 of the same study (header=echoes.header,
ImageType DERIVED) and with io.write_nifti, both in the input affine.
Read the DICOM back and check the study UID and the affine.