.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/02-parallel-imaging/02-nonlinear-inversion.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_02-parallel-imaging_02-nonlinear-inversion.py: =================== Nonlinear inversion =================== This lesson reconstructs an image and the coil sensitivities together from undersampled data whose fully sampled central region is too small for a separate calibration, and then writes the same reconstruction out as a nonlinear operator and a Gauss-Newton solver. ESPIRiT [#espirit]_ estimates the sensitivities from the autocalibration (ACS) region at the centre of k-space, and a linear SENSE reconstruction then treats them as known. When the ACS region is small, or absent, as in many real-time, non-Cartesian and highly accelerated protocols, the sensitivities are unknowns like the image, and the forward model .. math:: y_c = P F (S_c \cdot x) becomes bilinear: it is a product of two unknowns. Nonlinear inversion (NLINV) [#nlinv]_ solves it by the iteratively regularized Gauss-Newton method (IRGNM) [#bakushinsky]_, with the smoothness of the sensitivities, which resolves the ambiguity of the factorization, built into the model as a weighting of their k-space coefficients rather than added as a penalty. 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** - Reconstruct an image and its sensitivities jointly with :func:`bartorch.tools.nlinv` from an ACS region too small for ESPIRiT. - State the ambiguity of the bilinear factorization and the role of the Sobolev weighting of the sensitivities. - Write the same reconstruction as :class:`bartorch.nlop.NonlinearSense` under :class:`bartorch.nlop.IRGNM`. It follows :doc:`01-coil-calibration`, which used ``nlinv`` as a calibration step. The next lesson, :doc:`03-noise-prewhitening`, turns to the noise model of the receive channels. .. GENERATED FROM PYTHON SOURCE LINES 46-246 .. code-block:: Python import csv 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 nlop SIZE = 128 COILS = 8 ACCELERATION = 3 CALIBRATION = 6 # lines at the centre, far fewer than ESPIRiT needs pattern = lines.reshape(SIZE, 1).to(torch.complex64) measured = kspace[:, None] * pattern print(f"{float(lines.mean()):.0%} of the phase encodes, {CALIBRATION} of them at the centre") .. rst-class:: sphx-glr-script-out .. code-block:: none 32% of the phase encodes, 6 of them at the centre .. GENERATED FROM PYTHON SOURCE LINES 247-262 Six central lines locate the k-space centre but do not calibrate anything on their own. With ESPIRiT's default kernel of six points, a six-line ACS region leaves a single kernel position along the phase-encoding axis, too few rows for the calibration matrix of :func:`bartorch.tools.ecalib`. Joint reconstruction -------------------- :func:`bartorch.tools.nlinv` takes the k-space and returns the image and, when asked, the sensitivities it estimated along the way. Its iteration count is a number of Gauss-Newton steps rather than of linear iterations, and it acts as a regularization parameter rather than a convergence threshold: the regularization weight is halved after every step, so stopping early leaves a smoother image and running longer eventually lets the noise in. Eight steps is BART's default; twelve are used here. .. GENERATED FROM PYTHON SOURCE LINES 263-272 .. code-block:: Python STEPS = 12 reconstruction, estimated = bt.nlinv(measured, maxiter=STEPS, return_sensitivities=True) zero_filled = bartorch.rss(bartorch.ifft(measured[:, 0], axes=(-2, -1), unitary=True), axes=(0,)) print(f"NRMSE, zero-filled {bt.nrmse(image.abs(), zero_filled.abs(), scaled=True):.3f}") print(f"NRMSE, nlinv {bt.nrmse(image.abs(), reconstruction.abs(), scaled=True):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none NRMSE, zero-filled 0.209 NRMSE, nlinv 0.078 .. GENERATED FROM PYTHON SOURCE LINES 273-329 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_001.png :alt: reference, zero-filled, nlinv, |error|, nlinv :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_002.png :alt: simulated, channel 2, simulated, channel 4, simulated, channel 6, nlinv, channel 2, nlinv, channel 4, nlinv, channel 6 :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_003.png :alt: simulated, channel 2, simulated, channel 4, simulated, channel 6, nlinv, channel 2, nlinv, channel 4, nlinv, channel 6 :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_02-nonlinear-inversion_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 330-372 From a third of the phase encodes and six central lines, nonlinear inversion removes the aliasing of the zero-filled image; its error is concentrated at tissue boundaries and in the noise-like residue of the random sampling. The estimated sensitivities are smooth by construction rather than by agreement with the data: the coil unknown is not the sensitivity map but its k-space representation :math:`\hat{s}`, and the map follows as :math:`S = \mathcal{F}^{-1}[(1 + a|k|^2)^{-b/2} \hat{s}]`, a Sobolev-norm weighting that suppresses high spatial frequencies. A Gauss-Newton step in the unknown is therefore a smooth change of the map, and the joint problem needs no separate penalty on the coils. The pair is determined only up to a common factor: multiplying every map by a nonzero function :math:`\gamma(r)` and dividing the image by it leaves the data unchanged (:doc:`../../explanation/nonlinear`). The weighting restricts :math:`\gamma` to smooth functions, so the estimated maps match the simulated ones up to a smooth common magnitude and phase. The figures above therefore show the estimated maps inside the head, divided by their root sum of squares and with the phase of channel 0 subtracted, which removes that factor; so normalized, they reproduce the simulated maps. An ``nlinv`` image is reported after multiplication by the root sum of squares of the maps. Outside the object neither factor is determined at all -- their product is zero for any pair -- so the maps there follow from the initialization and the weighting. The model and the solver ------------------------ :class:`bartorch.nlop.NonlinearSense` is that forward model as a nonlinear operator with two inputs, the image and the coil coefficients, and :class:`bartorch.nlop.IRGNM` is the Gauss-Newton loop over it. Each step linearizes the model at the current estimate :math:`x_k` and solves .. math:: \min_x \, \| DF_{x_k} (x - x_k) - (y - F(x_k)) \|^2 + \alpha_k \| x - x_{\mathrm{ref}} \|^2, with :math:`DF_{x_k}` the derivative of the forward model, :math:`x_{\mathrm{ref}}` zero unless one is given, and :math:`\alpha_k` halved after every step, so the first steps are heavily regularized and the later ones are not. .. GENERATED FROM PYTHON SOURCE LINES 373-379 .. code-block:: Python model = nlop.NonlinearSense( (COILS, 1, SIZE, SIZE), pattern=lines.reshape(1, SIZE, 1).to(torch.complex64) ) print(f"inputs {model.ishapes} -> output {model.oshapes}") .. rst-class:: sphx-glr-script-out .. code-block:: none inputs ((1, 1, 128, 128), (8, 1, 128, 128)) -> output ((8, 1, 128, 128),) .. GENERATED FROM PYTHON SOURCE LINES 380-384 ``nlinv`` scales the data by ``100 / ||y||`` before it starts, which fixes the meaning of :math:`\alpha`, and runs the conjugate gradients of each step to a hundred iterations or a relative tolerance of a tenth. Given the same three settings, the loop written here is the application. .. GENERATED FROM PYTHON SOURCE LINES 385-398 .. code-block:: Python data = model.prepare(measured * (100.0 / float(torch.linalg.vector_norm(measured)))) fitted, coefficients = nlop.IRGNM(iterations=STEPS, cg_maxiter=100, cg_tol=0.1)(data, model) maps = model.coils(coefficients) combined = fitted.squeeze() * bartorch.rss(maps[:, 0], axes=(0,)) difference = float( (fitted.squeeze() - bt.nlinv(measured, maxiter=STEPS, normalize=False)).abs().max() ) print(f"largest difference from nlinv: {difference / float(fitted.abs().max()):.1e}") print(f"NRMSE {bt.nrmse(image.abs(), combined.abs(), scaled=True):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none largest difference from nlinv: 1.0e-03 NRMSE 0.078 .. GENERATED FROM PYTHON SOURCE LINES 399-414 The two agree to single-precision round-off rather than to the last bit, because the data scaling is computed here and inside the application by different expressions. What the operator form adds is access to everything around the step. The linearized problem can go to a solver from :mod:`bartorch.optim` instead of the conjugate gradients inside the library (``inner=optim.CG()`` is the same method written out, and a regularized solver makes the step a regularized one), the loop can be unrolled as :class:`bartorch.nlop.IRGNMBlock`, and a Gauss-Newton step is differentiable with respect to the data, the iterate, the regularization centre and :math:`\alpha` (:doc:`../../explanation/differentiation`). Reconstructing parameter maps rather than an image, by putting a signal model in front of the same encoding, is :doc:`../05-model-based/02-quantitative-models`. .. GENERATED FROM PYTHON SOURCE LINES 417-433 References ---------- .. [#espirit] Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, Vasanawala SS, Lustig M. ESPIRiT -- an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. *Magn Reson Med* 71(3):990-1001 (2014). https://doi.org/10.1002/mrm.24751 .. [#nlinv] Uecker M, Hohage T, Block KT, Frahm J. Image reconstruction by regularized nonlinear inversion -- joint estimation of coil sensitivities and image content. *Magn Reson Med* 60(3):674-682 (2008). https://doi.org/10.1002/mrm.21691 .. [#bakushinsky] Bakushinsky AB, Kokurin MY. *Iterative Methods for Approximate Solution of Inverse Problems.* Springer (2004). https://doi.org/10.1007/978-1-4020-3122-9 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.714 seconds) .. _sphx_glr_download_auto_examples_02-parallel-imaging_02-nonlinear-inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-nonlinear-inversion.ipynb <02-nonlinear-inversion.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-nonlinear-inversion.py <02-nonlinear-inversion.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-nonlinear-inversion.zip <02-nonlinear-inversion.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_