.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/03-regularization/01-regularized-reconstruction.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_03-regularization_01-regularized-reconstruction.py: ========================== Regularized reconstruction ========================== This lesson compares three regularization terms on the same undersampled, noisy SENSE acquisition, shows how the choice of the regularization weight trades residual noise and aliasing against loss of detail, and combines two terms in one reconstruction. At high acceleration, or low SNR, the SENSE problem is ill-conditioned: the encoding does not determine the image, and a plain least-squares fit amplifies the noise by the g-factor. A reconstruction therefore adds prior knowledge about the image as a penalty, .. math:: \hat{x} = \arg\min_x \tfrac12 \| A x - y \|_2^2 + \lambda\, g(G x), with :math:`A` the SENSE encoding, :math:`g` a convex functional, :math:`G` a linear transform and :math:`\lambda` the regularization weight. Tikhonov regularization favours a small image, a wavelet :math:`\ell_1` penalty an image that is sparse in a wavelet basis, as in compressed sensing, and total variation (TV) an image that is piecewise constant. :mod:`bartorch.priors` provides BART's terms :math:`g(Gx)` as objects, and :func:`bartorch.apps.pics` solves the problem with the iteration each term admits. :doc:`../../explanation/inverse-problems` introduces the formulation and the algorithms. **Learning objectives** - Pass regularization terms from :mod:`bartorch.priors` to :func:`bartorch.apps.pics`, and choose a solver the term admits. - Compare Tikhonov, wavelet :math:`\ell_1` and total-variation regularization on the same data, by their images and error maps. - Select a regularization weight by the error against a reference, and recognize under- and over-regularization in the image. - Combine two terms in one reconstruction. The previous lessons, :doc:`../01-basics/02-from-kspace-to-image` and :doc:`../02-parallel-imaging/01-coil-calibration`, reconstructed with a fixed weight. The next lesson, :doc:`02-operators-and-solvers`, assembles the same reconstruction from an operator, a term and a solver. .. GENERATED FROM PYTHON SOURCE LINES 47-165 .. 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 apps, priors SIZE = 192 COILS = 8 ACCELERATION = 3 CALIBRATION = 24 .. GENERATED FROM PYTHON SOURCE LINES 166-170 The phantom is the BrainWeb [#brainweb]_ slice of :doc:`../01-basics/02-from-kspace-to-image`, with the same eight-channel sensitivities; the cell that builds both is hidden on this page and present in the script this page can be downloaded as. .. GENERATED FROM PYTHON SOURCE LINES 171-242 .. GENERATED FROM PYTHON SOURCE LINES 243-253 Acquisition ----------- A third of the phase encodes (:math:`R = 3`), drawn at random from a variable density around a fully sampled ACS region of 24 lines, as in :doc:`../01-basics/02-from-kspace-to-image`. The noise is three times stronger than in that lesson: complex Gaussian noise of variance :math:`3 \times 10^{-4}` per sample of the unitary transform of an image whose peak is one. At this level noise amplification, and not only aliasing, determines the error of an unregularized reconstruction. .. GENERATED FROM PYTHON SOURCE LINES 254-271 .. code-block:: Python kspace = bt.noise(bartorch.fft(sensitivities * image, axes=(-2, -1), unitary=True), n=3e-4, s=42) encodes = torch.arange(SIZE) - SIZE // 2 centre = (encodes.abs() < CALIBRATION // 2).to(torch.float32) drawn = torch.multinomial( (1.0 + 2.0 * encodes.abs() / SIZE) ** -3.0 * (1.0 - centre), SIZE // ACCELERATION - CALIBRATION, replacement=False, generator=torch.Generator().manual_seed(11), ) lines = centre.clone() lines[drawn] = 1.0 measured = kspace[:, None] * lines.reshape(SIZE, 1).to(torch.complex64) maps = bt.ecalib(measured, maps=1, calib_size=CALIBRATION, crop=0.8) .. GENERATED FROM PYTHON SOURCE LINES 272-290 Three terms ----------- Tikhonov regularization, :math:`g(x) = \tfrac12\|x\|_2^2`, is the ``l2`` argument of ``pics`` and keeps the problem quadratic, so conjugate gradients solve it. The wavelet term :math:`\|\Psi x\|_1` promotes an image whose wavelet coefficients are sparse [#lustig]_; its transform is orthogonal and is applied inside its proximal operator, so FISTA [#beck]_ solves it. Total variation, :math:`\sum_r \|(\nabla x)_r\|_2`, promotes a piecewise-constant image [#block]_; the finite-difference operator :math:`\nabla` has no such closed-form proximal operator, so TV requires a splitting method, here ADMM [#boyd]_. The axes of a term are the tensor axes it acts along, here the two spatial axes. Each term is run over a range of weights. The weights are relative to the data divided by the scaling :func:`bartorch.optim.data_scaling` estimates, which ``pics`` applies, so the same weight means the same thing for data of a different overall scale. .. GENERATED FROM PYTHON SOURCE LINES 291-325 .. code-block:: Python sweeps = { "Tikhonov": [0.01, 0.03, 0.1, 0.3], "wavelet": [0.003, 0.006, 0.012, 0.03], "total variation": [0.002, 0.006, 0.012, 0.04], } def reconstruct(name, weight): if name == "Tikhonov": return apps.pics(measured, maps, l2=weight, maxiter=100) if name == "wavelet": term, solver = priors.Wavelet((-1, -2), weight), "fista" else: term, solver = priors.TotalVariation((-1, -2), weight), "admm" return apps.pics(measured, maps, regularizers=term, solver=solver, maxiter=100) reconstructions = { name: {weight: reconstruct(name, weight) for weight in weights} for name, weights in sweeps.items() } errors_by_weight = { name: {w: bt.nrmse(image.abs(), x.abs(), scaled=True) for w, x in results.items()} for name, results in reconstructions.items() } best = {name: min(values, key=values.get) for name, values in errors_by_weight.items()} for name, weight in best.items(): estimate = reconstructions[name][weight] error = errors_by_weight[name][weight] similarity = bt.ssim(image.abs(), scaled(estimate, image)) print(f"{name:>16} weight {weight:<6} NRMSE {error:.3f} SSIM {similarity:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Tikhonov weight 0.03 NRMSE 0.082 SSIM 0.783 wavelet weight 0.006 NRMSE 0.056 SSIM 0.867 total variation weight 0.006 NRMSE 0.059 SSIM 0.926 .. GENERATED FROM PYTHON SOURCE LINES 326-352 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_001.png :alt: reference, Tikhonov, $\lambda$ = 0.03, wavelet, $\lambda$ = 0.006, total variation, $\lambda$ = 0.006 :srcset: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_002.png :alt: Tikhonov error, wavelet error, total variation error :srcset: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_003.png :alt: reference, enlarged, Tikhonov, wavelet, total variation :srcset: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 353-373 The Tikhonov reconstruction retains amplified noise and the incoherent aliasing of the random sampling across the whole head: a quadratic penalty that is strong enough to suppress them also blurs the image, so its best weight leaves them in. Both sparsity-promoting terms remove most of the noise while keeping the tissue boundaries, which the error maps show as a much darker background inside the brain. They differ in their residual artefacts: in the enlarged region the wavelet penalty leaves a blotchy residual texture in white matter, and total variation renders the gradual intensity variations within white matter as patches of constant intensity (staircasing). The regularization weight ------------------------- A weight that is too small leaves the noise in; one that is too large removes image detail with it, and the error against the reference has a minimum between the two. The sweep above spans a factor of ten or more for each term. It is possible here because the phantom is known; for measured data the weight is chosen by a criterion that does not require the reference, or fixed once for a protocol. .. GENERATED FROM PYTHON SOURCE LINES 374-395 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_004.png :alt: 01 regularized reconstruction :srcset: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_005.png :alt: TV, $\lambda$ = 0.002 too small, TV, $\lambda$ = 0.006 best, TV, $\lambda$ = 0.04 too large :srcset: /auto_examples/03-regularization/images/sphx_glr_01-regularized-reconstruction_005.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 396-412 Each curve has an interior minimum. The weights of the Tikhonov term and of the two :math:`\ell_1` terms are not comparable with each other, because the functionals differ; only the location of each minimum is meaningful. At its minimum each sparsity-promoting term reaches a lower error than Tikhonov regularization at its own minimum, for this image and this noise level. The three TV reconstructions show the two failure modes. At the smallest weight the noise and the incoherent aliasing remain; at the largest the cortex is flattened into patches of constant intensity and thin gyri merge. Combining terms --------------- ``regularizers`` accepts a list, whose terms are summed. ADMM splits the variable once per term with a nontrivial transform, so it accepts any combination; FISTA accepts only terms whose transform is the identity. .. GENERATED FROM PYTHON SOURCE LINES 413-426 .. code-block:: Python combined = apps.pics( measured, maps, regularizers=[ priors.Wavelet((-1, -2), best["wavelet"] / 2), priors.TotalVariation((-1, -2), best["total variation"] / 2), ], solver="admm", maxiter=100, ) print(f"wavelet + TV NRMSE {bt.nrmse(image.abs(), combined.abs(), scaled=True):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none wavelet + TV NRMSE 0.058 .. GENERATED FROM PYTHON SOURCE LINES 427-435 With each weight halved the sum reaches an error comparable to either term alone. Whether a combination improves on its parts depends on the image and the weights, and is established by a comparison such as the one above rather than assumed. :func:`bartorch.apps.pics` builds three objects -- the encoding operator, the terms and the iteration -- and runs BART's solver on them. The next lesson, :doc:`02-operators-and-solvers`, builds them separately. .. GENERATED FROM PYTHON SOURCE LINES 438-461 References ---------- .. [#brainweb] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. *IEEE Trans Med Imaging* 17(3):463-468 (1998). https://doi.org/10.1109/42.712135 .. [#lustig] Lustig M, Donoho D, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. *Magn Reson Med* 58(6):1182-1195 (2007). https://doi.org/10.1002/mrm.21391 .. [#beck] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. *SIAM J Imaging Sci* 2(1):183-202 (2009). https://doi.org/10.1137/080716542 .. [#block] Block KT, Uecker M, Frahm J. Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. *Magn Reson Med* 57(6):1086-1098 (2007). https://doi.org/10.1002/mrm.21236 .. [#boyd] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Found Trends Mach Learn* 3(1):1-122 (2011). https://doi.org/10.1561/2200000016 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.017 seconds) .. _sphx_glr_download_auto_examples_03-regularization_01-regularized-reconstruction.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-regularized-reconstruction.ipynb <01-regularized-reconstruction.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-regularized-reconstruction.py <01-regularized-reconstruction.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-regularized-reconstruction.zip <01-regularized-reconstruction.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_