.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/08-workflows/03-maps-from-scanner-images.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_08-workflows_03-maps-from-scanner-images.py: ================================== Parameter maps from scanner images ================================== In this tour you estimate a :math:`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. :func:`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: :func:`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 :math:`T_2`. With ``magnitude=True``, :func:`~bartorch.apps.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, :doc:`../05-model-based/01-quantitative-models`. **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 :class:`bartorch.nlop.MultiEcho` to magnitude images with :func:`bartorch.apps.mobafit`. - Write the map as a DICOM series of the same study, and as NIfTI, in the geometry of the input. .. GENERATED FROM PYTHON SOURCE LINES 42-110 .. GENERATED FROM PYTHON SOURCE LINES 111-128 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 ------------- :func:`~bartorch.io.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. .. GENERATED FROM PYTHON SOURCE LINES 128-150 .. code-block:: Python 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") .. image-sg:: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_001.png :alt: TE 10 ms, TE 40 ms, TE 80 ms :srcset: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 151-165 Three of the echoes, on one window. The cerebrospinal fluid in the ventricles, with the longest :math:`T_2`, keeps its signal across the train, while the scalp fat loses most of it. Decay fitting ------------- Build :class:`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. .. GENERATED FROM PYTHON SOURCE LINES 165-176 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 177-185 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. .. GENERATED FROM PYTHON SOURCE LINES 185-201 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none PosixPath('/tmp/tmpzq1oi1vm/t2-map.nii.gz') .. GENERATED FROM PYTHON SOURCE LINES 202-203 Reading the series back confirms the study and the geometry: .. GENERATED FROM PYTHON SOURCE LINES 203-227 .. code-block:: Python 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") .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_002.png :alt: simulated, fitted :srcset: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_003.png :alt: |fitted - simulated| :srcset: /auto_examples/08-workflows/images/sphx_glr_03-maps-from-scanner-images_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none patient P0001, same study True largest affine difference 2.0e-06 mm .. GENERATED FROM PYTHON SOURCE LINES 228-249 The map you read back from DICOM, beside the :math:`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: .. code-block:: text 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. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.079 seconds) .. _sphx_glr_download_auto_examples_08-workflows_03-maps-from-scanner-images.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-maps-from-scanner-images.ipynb <03-maps-from-scanner-images.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-maps-from-scanner-images.py <03-maps-from-scanner-images.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-maps-from-scanner-images.zip <03-maps-from-scanner-images.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_