.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/05-model-based/02-quantitative-models.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_05-model-based_02-quantitative-models.py: ==================================== Parameter maps straight from k-space ==================================== This lesson estimates a :math:`T_2` map from an undersampled multi-echo spin-echo acquisition in two ways, and compares them: reconstructing an image per echo and fitting the decay voxel by voxel afterwards, and fitting the signal model directly to the k-space data. The aim is to show why the second, model-based reconstruction, tolerates undersampling that ruins the first. In a multi-echo spin-echo (CPMG) acquisition the signal of each voxel decays from echo to echo as :math:`M_0 \exp(-\mathrm{TE}/T_2)`. Undersampling each echo shortens the scan, but a reconstruction of each echo on its own is an ill-posed problem, and its aliasing and noise differ from echo to echo; a voxelwise fit cannot tell them apart from decay, and carries them into the map. The model-based approach [#sumpf]_ [#wang]_ puts the signal model inside the forward operator, .. math:: y_{c,e} = P_e F \, (S_c \cdot M_e(\theta)), where :math:`P_e` is the sampling pattern of echo :math:`e`, :math:`F` the Fourier transform, :math:`S_c` the sensitivity of coil :math:`c`, :math:`M` the signal model and :math:`\theta` the parameter maps, and solves for :math:`\theta` from the k-space data of all echoes at once. The unknowns are then three real maps rather than eight complex images, and every echo constrains all of them. The operator is nonlinear in :math:`\theta`, so the problem is solved by the iteratively regularized Gauss-Newton method of :doc:`../02-parallel-imaging/02-nonlinear-inversion`, over a different model. The model here is :class:`bartorch.nlop.MultiEcho`, a TorchSim simulator as a BART nonlinear operator. The phantom and the coil sensitivities are built as in :doc:`../01-basics/02-from-kspace-to-image`; the cell that does it is hidden on this page and present in the script this page can be downloaded as. **Learning objectives** - Represent a relaxation model as a TorchSim-backed :class:`bartorch.nlop.SignalModel`. - Fit it to reconstructed echo images, and directly to k-space by composing it with the encoding, with :class:`bartorch.nlop.IRGNM`. - Run the same fits through :func:`bartorch.apps.mobafit` and :func:`bartorch.apps.moba`. - Explain, from the echo images and the error maps, why the model-based fit is more accurate at the same undersampling. It follows :doc:`01-subspace-t1-mapping`, and :doc:`03-maps-from-scanner-images` fits the same model to images read from DICOM. .. GENERATED FROM PYTHON SOURCE LINES 55-173 .. code-block:: Python import csv import time from pathlib import Path import brainweb_dl import numpy as np import torch from brainweb_dl import get_mri import bartorch import bartorch.tools as bt from bartorch import apps, linop, nlop, optim SIZE = 64 COILS = 8 ECHOES = 8 ACCELERATION = 4 ECHO_TIMES = torch.tensor([12.5 * (echo + 1) for echo in range(ECHOES)]) # ms .. GENERATED FROM PYTHON SOURCE LINES 174-181 Phantom ------- The :math:`T_2` of each tissue class from the BrainWeb table, combined by membership, and the echo images from the mono-exponential decay :math:`M_0 \exp(-\mathrm{TE}/T_2)` written out here rather than taken from the model that will be fitted. .. GENERATED FROM PYTHON SOURCE LINES 182-254 .. code-block:: Python contrasts = (amplitude[None] * torch.exp(-ECHO_TIMES[:, None, None] / t2[None])).to(torch.complex64) .. GENERATED FROM PYTHON SOURCE LINES 255-266 Acquisition ----------- Eight echoes at an echo spacing of 12.5 ms, each sampled at a quarter of the phase encodes (:math:`R = 4`) around eight fully sampled central lines, with a different random draw per echo, so that the missing phase encodes differ between echoes. The echoes are a batch of the encoding rather than an axis inside it: the sensitivities and the transform are shared, and only the pattern differs, so the operator is the Cartesian SENSE encoding of :doc:`../03-regularization/02-operators-and-solvers` with the pattern of each echo applied to its samples. .. GENERATED FROM PYTHON SOURCE LINES 267-308 .. code-block:: Python encoding = linop.CartesianSense(sensitivities, (ECHOES, SIZE, SIZE), ndim=2) E = linop.Diagonal(lines.to(torch.complex64), encoding.oshape) @ encoding measured = bt.noise(E(contrasts), n=1e-6, s=9) data = measured / float(E.H(measured).abs().max()) print(f"{E.ishape} -> {E.oshape}") .. image-sg:: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_001.png :alt: sampled phase encodes (white) :srcset: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (8, 64, 64) -> (8, 8, 64, 64) .. GENERATED FROM PYTHON SOURCE LINES 309-317 The signal model ---------------- :class:`bartorch.nlop.MultiEcho` maps parameter maps to one image per echo. Its unknowns are :math:`T_2` and a complex amplitude, carried as three real maps in a bounded parameterisation rather than in their own units, so :meth:`~bartorch.nlop.SignalModel.initial` builds a starting point from values and :meth:`~bartorch.nlop.SignalModel.split` reads the fit back. .. GENERATED FROM PYTHON SOURCE LINES 318-324 .. code-block:: Python M = nlop.MultiEcho([float(te) for te in ECHO_TIMES], (SIZE, SIZE)) start = M.initial(T2=80.0) print(f"unknowns {M.names}: {M.ishapes[0]} -> {M.oshapes[0]}") .. rst-class:: sphx-glr-script-out .. code-block:: none unknowns ('T2', 'amplitude.real', 'amplitude.imag'): (3, 64, 64) -> (8, 64, 64) .. GENERATED FROM PYTHON SOURCE LINES 325-335 Two routes ---------- The first reconstructs the echo images by conjugate gradients on the SENSE normal equations and fits the model to them voxel by voxel, which is what :func:`bartorch.apps.mobafit` does given the images. The second composes the model with the encoding and fits the k-space with :class:`bartorch.nlop.IRGNM`. Both are Gauss-Newton loops of the same number of steps and differ only in the forward operator that maps the unknowns to the data. .. GENERATED FROM PYTHON SOURCE LINES 336-348 .. code-block:: Python STEPS = 20 start_time = time.perf_counter() images = optim.CG(maxiter=40)(data, E) two_step = apps.mobafit(images, M, iterations=STEPS, T2=80.0) print(f"reconstruct, then fit: {time.perf_counter() - start_time:5.1f} s") start_time = time.perf_counter() model_based = nlop.IRGNM(iterations=STEPS, cg_maxiter=100, cg_tol=0.1)(data, E @ M, x0=start) print(f"model inside the operator: {time.perf_counter() - start_time:5.1f} s") .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/torch/jit/_script.py:1491: FutureWarning: `torch.jit.script` is deprecated. Please switch to `torch.compile` or `torch.export`. warnings.warn( reconstruct, then fit: 3.8 s model inside the operator: 9.6 s .. GENERATED FROM PYTHON SOURCE LINES 349-355 ``E @ M`` composes a linear operator with a nonlinear one; the derivative of the composition at a point is the encoding applied to the derivative of the model, which is the derivative a Gauss-Newton step requires. The echo images of the first route show what its fit is given. They are compared here with the fully sampled echo images of the phantom. .. GENERATED FROM PYTHON SOURCE LINES 356-371 .. image-sg:: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_002.png :alt: TE = 12.5 ms, TE = 50.0 ms, TE = 100.0 ms :srcset: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 372-383 Each reconstructed echo carries residual aliasing and noise, which differ from echo to echo because each echo has its own sampling pattern. At the later echoes the signal has decayed and the relative error grows, so the late echoes, which determine :math:`T_2` most, are also the least accurate. :func:`bartorch.apps.moba` assembles the model-based composition from the k-space, the model, the sensitivities and the sampling pattern, and returns the maps in their own units. It scales the data by the rule of :func:`bartorch.optim.data_scaling` and regularizes each step towards the starting maps rather than towards zero, so its result is not identical to the fit above. .. GENERATED FROM PYTHON SOURCE LINES 384-402 .. code-block:: Python start_time = time.perf_counter() one_call = apps.moba(measured, M, sensitivities, pattern=lines, iterations=STEPS, T2=80.0) print(f"apps.moba: {time.perf_counter() - start_time:5.1f} s") estimates = { "two-step": two_step["T2"], "model-based": M.split(model_based)["T2"], "apps.moba": one_call["T2"], } for name, estimate in estimates.items(): error = float((estimate[support] - t2[support]).norm() / t2[support].norm()) median = float(estimate[support].median()) print(f"{name:>22} median {median:5.1f} ms relative error {error:.3f}") print(f"{'phantom':>22} median {float(t2[support].median()):5.1f} ms") .. rst-class:: sphx-glr-script-out .. code-block:: none apps.moba: 10.3 s two-step median 86.1 ms relative error 0.477 model-based median 81.5 ms relative error 0.007 apps.moba median 81.3 ms relative error 0.057 phantom median 81.5 ms .. GENERATED FROM PYTHON SOURCE LINES 403-421 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_003.png :alt: reference, two-step, model-based, apps.moba :srcset: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_004.png :alt: $T_2$ error, two-step, model-based, apps.moba :srcset: /auto_examples/05-model-based/images/sphx_glr_02-quantitative-models_004.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 422-444 The two-step :math:`T_2` map is dominated by the errors of the echo images: a voxelwise fit cannot distinguish residual aliasing from decay, and in voxels where a late echo is too bright or too dark the fitted :math:`T_2` is far off. The model-based fits reach a much lower error from the same data, since the model admits only images that decay exponentially from echo to echo, and the aliasing of eight different sampling patterns is not such an image. The fit inside the operator reproduces the phantom almost exactly, because every voxel of the phantom decays with a single :math:`T_2`, as the model assumes; a measured voxel holding two tissues decays with two, and a single exponential cannot represent it. :func:`bartorch.apps.moba` regularizes each Gauss-Newton step towards the starting maps, and its residual error is largest where :math:`T_2` is farthest from the starting value of 80 ms: in cerebrospinal fluid and the scalp. The maps are drawn with the navia colormap [#fuderer]_ in a window that spans white and grey matter. Since the maps are the solver's unknowns, a regularizer passed to the linearized problem -- the ``inner`` solver of :func:`bartorch.apps.moba` -- penalizes the maps rather than the echo images. Without ``sensitivities``, :func:`bartorch.apps.moba` estimates the coils jointly with the maps, as :doc:`../02-parallel-imaging/02-nonlinear-inversion` estimates them jointly with an image. :doc:`../../explanation/nonlinear` explains the model-based approach in more detail. .. GENERATED FROM PYTHON SOURCE LINES 447-462 References ---------- .. [#sumpf] Sumpf TJ, Uecker M, Boretius S, Frahm J. Model-based nonlinear inverse reconstruction for T2 mapping using highly undersampled spin-echo MRI. *J Magn Reson Imaging* 34(2):420-428 (2011). https://doi.org/10.1002/jmri.22634 .. [#wang] Wang X, Tan Z, Scholand N, Roeloffs V, Uecker M. Physics-based reconstruction methods for magnetic resonance imaging. *Phil Trans R Soc A* 379(2200):20200196 (2021). https://doi.org/10.1098/rsta.2020.0196 .. [#fuderer] Fuderer M, Wichtmann B, Crameri F, de Souza NM, Baeßler B, Gulani V, et al. Color-map recommendation for MR relaxometry maps. *Magn Reson Med* 93(2):490-506 (2025). https://doi.org/10.1002/mrm.30290 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 24.581 seconds) .. _sphx_glr_download_auto_examples_05-model-based_02-quantitative-models.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-quantitative-models.ipynb <02-quantitative-models.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-quantitative-models.py <02-quantitative-models.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-quantitative-models.zip <02-quantitative-models.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_