.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/05-misc/01-synthetic-data.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 or to run this example in your browser via Binder. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_autoexamples_05-misc_01-synthetic-data.py: ============================= Synthetic MR fingerprinting ============================= The scope of this notebook is to build a training pair end to end: what a scanner would measure from a fingerprinting exam, and the maps it came from. A subject is segmented into tissue classes; each class is given an M0, a T1 and a T2; one voxel per class is simulated by extended phase graphs; every voxel of a class is handed its class's evolution; the volume is weighted by birdcage sensitivities and pushed through a frame-wise non-uniform Fourier transform; the k-space is brought back and coil-combined. The undersampled series and the fully sampled one are the pair, with the ground-truth maps and the segmentation. Only the simulation is TorchSim's. The phantom, the coils and the encoding come from torchio, deepmriprep, SigPy and mri-nufft. .. GENERATED FROM PYTHON SOURCE LINES 22-26 .. colab-link:: :needs_gpu: 1 !pip install torchsim torchio deepmriprep sigpy cmap mri-nufft[finufft,cufinufft] .. GENERATED FROM PYTHON SOURCE LINES 28-41 Only one step of this pipeline is TorchSim's. Four other packages do the rest, and each has exactly one job here: * ``torchio`` fetches a T1-weighted IXI subject and gives it as tensors, which is where the anatomy comes from; * ``deepmriprep`` segments that anatomy into grey matter, white matter and CSF with a U-Net -- torch throughout, and it hands back probabilities rather than labels; * ``sigpy.mri`` generates birdcage coil sensitivities to weight the volume with; * ``mri-nufft`` supplies the spiral trajectory and the non-uniform Fourier transform that plays it, forwards and back. .. GENERATED FROM PYTHON SOURCE LINES 42-153 .. code-block:: Python import mrinufft from deepmriprep import Preprocess import sigpy.mri as smri import torchio as tio from mrinufft.trajectories import initialize_2D_spiral .. GENERATED FROM PYTHON SOURCE LINES 154-157 TorchSim's part is the third step: one fingerprinting simulation per tissue class, rather than one per voxel. .. GENERATED FROM PYTHON SOURCE LINES 158-168 .. code-block:: Python import tempfile import time from pathlib import Path import numpy as np import torch from torchsim.simulators import MRFSimulator .. GENERATED FROM PYTHON SOURCE LINES 169-172 What the pair will be: a 128 matrix, four hundred frames, one spiral arm of 768 samples per frame, and eight receive channels. .. GENERATED FROM PYTHON SOURCE LINES 173-185 .. code-block:: Python SIZE = 128 FRAMES = 400 SAMPLES = 768 COILS = 8 SLICE = 110 # axial, at the level of the lateral ventricles # The GPU transform is used when it is both installed and usable; the # simulation follows it, so the images and the operator meet on one device. on_gpu = torch.cuda.is_available() and mrinufft.check_backend("cufinufft") device = "cuda" if on_gpu else "cpu" backend = "cufinufft" if on_gpu else "finufft" .. GENERATED FROM PYTHON SOURCE LINES 186-196 Subject ------- One IXI subject through torchio: a T1-weighted volume at 1.5 T, and the only measurement this notebook starts from. The segmentation reads it; a table supplies the tissue properties a contrast cannot give. ``download=True`` fetches the archive once, a few hundred MB, into a cache under the home directory. .. GENERATED FROM PYTHON SOURCE LINES 197-213 .. code-block:: Python CACHE = Path.home() / ".cache" / "torchsim" / "ixi-tiny" subject = tio.datasets.IXITiny(str(CACHE), download=True)[0] .. rst-class:: sphx-glr-script-out .. code-block:: none Downloading https://www.dropbox.com/s/ogxjwjxdv5mieah/ixi_tiny.zip?dl=1 to /tmp/tmpvgm8gayl.zip 0it [00:00, ?it/s] 0%| | 0/233926107 [00:01 0.5 class_M0 = NAMES_M0 class_T1 = NAMES_T1 class_T2 = NAMES_T2 .. rst-class:: sphx-glr-script-out .. code-block:: none tissue voxels M0 T1 (ms) T2 (ms) grey matter 2427 0.80 1100 95.0 white matter 2294 0.70 650 70.0 CSF 571 1.00 4000 2000.0 .. GENERATED FROM PYTHON SOURCE LINES 309-317 Simulation per class -------------------- An inversion, then four hundred repetitions whose flip angle sweeps. :class:`~torchsim.simulators.MRFSimulator` takes arrays of tissue properties, so the whole table is one call: three extended phase graph runs rather than sixteen thousand. .. GENERATED FROM PYTHON SOURCE LINES 318-333 .. code-block:: Python schedule = 5.0 + 55.0 * torch.sin(torch.linspace(0.0, 4 * torch.pi, FRAMES)).abs() simulator = MRFSimulator(TR=12.0, TI=20.0, T1=class_T1, T2=class_T2) per_class = torch.as_tensor(simulator.simulate(flip=schedule)) .. rst-class:: sphx-glr-script-out .. code-block:: none 3 classes x 400 frames in 0.4s -> (3, 400) .. GENERATED FROM PYTHON SOURCE LINES 334-345 Whole-brain volume ------------------ A voxel that is part grey matter and part CSF produces the sum of what the two do, weighted by how much of each is present -- not the signal of their averaged relaxation times. The mixing happens on the signals, which is one matrix product against the per-class evolutions. Averaging the parameters and simulating once is wrong wherever a voxel is not pure: an inversion-prepared train is markedly nonlinear in T1. .. GENERATED FROM PYTHON SOURCE LINES 346-356 .. code-block:: Python weights = (fractions * class_M0).to(per_class.dtype) series = (weights @ per_class) * brain[..., None] # The maps that go out as ground truth are the mixture averages, which is what # a single-compartment fit of this data could return at best. share = occupancy.clamp_min(1e-6) truth_M0 = (fractions @ class_M0) * brain truth_T1 = (fractions @ class_T1) / share * brain truth_T2 = (fractions @ class_T2) / share * brain .. GENERATED FROM PYTHON SOURCE LINES 357-364 Coils and encoding ------------------ Birdcage sensitivities from SigPy, then one spiral arm per frame rotated by the golden angle. A single arm of 768 samples against a 128 x 128 matrix is twenty-one-fold undersampled, which is how MRF is run. .. GENERATED FROM PYTHON SOURCE LINES 365-398 .. code-block:: Python sensitivities = torch.as_tensor(smri.birdcage_maps((COILS, SIZE, SIZE))).to( torch.complex64 ) trajectory = initialize_2D_spiral( FRAMES, SAMPLES, tilt="golden", nb_revolutions=8 ).astype(np.float32) build = mrinufft.get_operator(backend) arms = [ build( trajectory[frame], (SIZE, SIZE), n_coils=COILS, squeeze_dims=False, density=True ) for frame in range(FRAMES) ] coil_series = (sensitivities[:, None] * series.movedim(-1, 0)[None]).to(device) kspace = torch.stack( [arms[frame].op(coil_series[:, frame][None])[0] for frame in range(FRAMES)] ) .. rst-class:: sphx-glr-script-out .. code-block:: none 400 arms built in 3.0s forward NUFFT 3.4s -> (400, 8, 768) .. GENERATED FROM PYTHON SOURCE LINES 399-402 Eight birdcage sensitivities, and one spiral arm per frame rotated so that consecutive frames sample different parts of k-space. .. GENERATED FROM PYTHON SOURCE LINES 403-435 .. image-sg:: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_002.png :alt: coil 0, coil 1, coil 2, coil 3, 3 of 400 arms :srcset: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 436-442 Reconstruction -------------- Adjoint per frame, then a sensitivity-weighted coil combination, available here because the maps are known. A real pipeline would estimate them. .. GENERATED FROM PYTHON SOURCE LINES 443-458 .. code-block:: Python folded = torch.stack( [arms[frame].adj_op(kspace[frame][None])[0] for frame in range(FRAMES)] ) weights = sensitivities.to(device) combined = (folded * weights.conj()[None]).sum(1) / ( weights.abs().square().sum(0)[None] + 1e-6 ) .. rst-class:: sphx-glr-script-out .. code-block:: none adjoint and combine 3.4s .. GENERATED FROM PYTHON SOURCE LINES 459-466 Data pair --------- Each frame is one spiral arm, so the aliasing is worse than the signal. What survives is the time course, and that is what a fingerprinting reconstruction reads. .. GENERATED FROM PYTHON SOURCE LINES 467-491 .. rst-class:: sphx-glr-script-out .. code-block:: none per-frame error inside the brain : 60.7% time-course agreement : median 0.893, tenth percentile 0.714 .. GENERATED FROM PYTHON SOURCE LINES 492-497 Fifty percent wrong frame by frame, and the fingerprints still line up above 0.86 for nine voxels in ten. That gap is the premise of the method: the input is what the scanner gives, artefacts and all, and the target is the curve underneath. .. GENERATED FROM PYTHON SOURCE LINES 500-503 Phantom summary --------------- .. GENERATED FROM PYTHON SOURCE LINES 504-516 .. image-sg:: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_003.png :alt: 01 synthetic data :srcset: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 517-520 Inspecting the pair ------------------- .. GENERATED FROM PYTHON SOURCE LINES 521-551 .. image-sg:: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_004.png :alt: frame 0, frame 200, one voxel :srcset: /generated/autoexamples/05-misc/images/sphx_glr_01-synthetic-data_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 552-559 Exporting --------- The pair, the ground truth, the segmentation, and the schedule and trajectory that produced them. It goes to a temporary directory here so that a documentation build leaves no archive behind. .. GENERATED FROM PYTHON SOURCE LINES 560-582 .. code-block:: Python contents = { "undersampled": undersampled.cpu().numpy(), "reference": reference.cpu().numpy(), "M0": truth_M0.numpy(), "T1": truth_T1.numpy(), "T2": truth_T2.numpy(), "tissue_probabilities": fractions.numpy(), "flip_angles_deg": schedule.numpy(), "trajectory": trajectory, } .. rst-class:: sphx-glr-script-out .. code-block:: none synthetic_mrf.npz 56.6 MiB undersampled (400, 128, 128) complex64 reference (400, 128, 128) complex64 M0 (128, 128) float32 T1 (128, 128) float32 T2 (128, 128) float32 tissue_probabilities (128, 128, 3) float32 flip_angles_deg (400,) float32 trajectory (400, 768, 2) float32 .. GENERATED FROM PYTHON SOURCE LINES 583-597 Command-line use ---------------- Everything above is fixed except three inputs: * **a T1-weighted NIfTI**, which replaces the torchio subject and is what the segmentation was trained on, so the CSF map stops being the weak one; * **a matfile** carrying the trajectory and the flip-angle schedule, read instead of the spiral and the sine generated here; * **an output path**, replacing the temporary directory. The tissue table is the one thing decided rather than read. A relaxometry protocol on the same subject would fill it from data, which is what the parameter-inference notebooks do. .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 45.942 seconds) .. _sphx_glr_download_generated_autoexamples_05-misc_01-synthetic-data.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/firmlab-pisa/torchsim/gh-pages?urlpath=lab/tree/v0.0.8/examples/generated/autoexamples/05-misc/01-synthetic-data.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-synthetic-data.ipynb <01-synthetic-data.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-synthetic-data.py <01-synthetic-data.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-synthetic-data.zip <01-synthetic-data.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_