.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/04-non-cartesian/02-radial-sense.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_04-non-cartesian_02-radial-sense.py: =========================== Radial SENSE reconstruction =========================== This lesson reconstructs an undersampled golden-angle radial acquisition with eight receive coils: the density-compensated gridding reconstruction first, then an iterative SENSE reconstruction with coil sensitivities estimated from the radial data themselves, with and without a total-variation penalty. The aim is to see which of the streak artefacts of radial undersampling the coil encoding removes, which the regularization removes, and what each costs. A radial acquisition that satisfies the Nyquist criterion at the edge of k-space needs :math:`\pi/2` times as many spokes as the matrix has lines; with fewer, the azimuthal gaps between spokes alias into streaks that run across the whole field of view. Radial undersampling is nevertheless benign compared with Cartesian undersampling: every spoke passes through the k-space centre, so the low spatial frequencies stay fully sampled and the aliasing is incoherent rather than a discrete fold-over. Parallel imaging removes it by fitting the image to the non-Cartesian SENSE model [#pruessmann2001]_ .. math:: A = W \, \mathrm{NUFFT} \, S, with :math:`S` the coil sensitivities, the NUFFT evaluated along the trajectory and :math:`W` an optional weighting of the samples. The fit is solved iteratively, since :math:`A^H A` is not diagonal in any basis. The measured data are simulated with the same transform the reconstruction uses, so the comparison isolates the undersampling, the noise and the error of the estimated sensitivities from any mismatch between the forward model and the measurement. 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** - Simulate a multichannel radial acquisition with :class:`bartorch.linop.NoncartesianSense`. - Estimate sensitivities from the radial data with :func:`bartorch.tools.ncalib`. - Compare gridding, unregularized CG-SENSE and total-variation-regularized SENSE, and identify the residual artefact of each. - Reconstruct with :func:`bartorch.apps.pics` and with the operator under :class:`bartorch.optim.ADMM`, and compare the two forms of the normal operator. It follows :doc:`01-trajectories-and-transforms`. The next lesson, :doc:`03-dynamic-golden-angle`, adds a time axis. .. GENERATED FROM PYTHON SOURCE LINES 55-242 .. 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, optim, priors SIZE = 192 COILS = 8 SPOKES = 48 # against pi/2 * SIZE = 302 for a trajectory that is not undersampled TV_WEIGHT = 0.0005 ITERATIONS = 30 .. GENERATED FROM PYTHON SOURCE LINES 243-257 Acquisition ----------- Forty-eight golden-angle spokes of 192 samples each, across a 192 matrix: :math:`\pi/2 \times 192 \approx 302` spokes would sample the edge of k-space at the Nyquist rate, so the acquisition is undersampled by a factor of about 6.3. Successive spokes are rotated by the golden angle, 111.25 degrees, so any contiguous subset of them covers k-space nearly uniformly [#winkelmann]_; :doc:`03-dynamic-golden-angle` relies on that property. :class:`bartorch.linop.NoncartesianSense` maps an image to the samples of every channel along the trajectory. The measurement is that operator applied to the phantom, with complex Gaussian noise of variance :math:`10^{-4}` per sample added, as in the Cartesian lessons. .. GENERATED FROM PYTHON SOURCE LINES 258-267 .. code-block:: Python trajectory = bt.traj(readout=SIZE, spokes=SPOKES, radial=True, golden=True) E = linop.NoncartesianSense(sensitivities, (SIZE, SIZE), traj=trajectory) measured = bt.noise(E(image), n=1e-4, s=7) print(f"{E.ishape} -> {E.oshape}") print(E.plan) .. rst-class:: sphx-glr-script-out .. code-block:: none (192, 192) -> (8, 48, 192) Plan(transform=nufft, image=sensitivities, normal=kernel, coil_batch=1, streamed=coils, executor=slab) .. GENERATED FROM PYTHON SOURCE LINES 268-288 The operator's samples are ``(coils, spokes, samples)``. BART's applications carry non-Cartesian k-space in their own layout, ``(coils, spokes, samples, 1)``, whose trailing axis is the readout dimension of a Cartesian acquisition, so an application is given ``measured[..., None]``. Sensitivity calibration ----------------------- ESPIRiT reads its calibration matrix from a Cartesian neighbourhood of the k-space centre, so on radial data it needs that region gridded first. :func:`bartorch.tools.ncalib` estimates the sensitivities from the samples as they were measured, by nonlinear inversion [#nlinv]_ at low resolution; the densely sampled centre of a radial acquisition acts as its own autocalibration region. ``N=True`` divides the estimated maps by their root sum of squares. A SENSE fit recovers the image :math:`x` for which :math:`Sx` explains the data, so maps whose root sum of squares varies across the field of view leave its reciprocal in the image as a smooth intensity shading. ESPIRiT maps are normalized by construction; nonlinear inversion maps are not. .. GENERATED FROM PYTHON SOURCE LINES 289-303 .. code-block:: Python maps = bt.ncalib(measured[..., None], t=trajectory, N=True) .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_001.png :alt: coil 1, coil 4, coil 7 :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 304-318 The estimated maps reproduce the magnitude and phase of the simulated ones over the head. They are smoother, because nonlinear inversion penalizes the high spatial frequencies of the sensitivities, and outside the head, where there is no signal to calibrate from, they are extrapolated. Gridding -------- The gridding reconstruction is the density-compensated adjoint: each sample is weighted by its distance from the k-space centre (the ramp filter of filtered back-projection), the samples are interpolated onto the grid by the adjoint NUFFT, and the channels are combined by root sum of squares. It uses no model of the coil encoding, so the missing spokes appear in it as the streaks the point spread function of the trajectory predicts. .. GENERATED FROM PYTHON SOURCE LINES 319-326 .. code-block:: Python weights = torch.linalg.norm(trajectory.real[..., :2], dim=-1, keepdim=True) weights = weights.clamp(min=0.25).to(torch.complex64) channels = bartorch.nufft_adjoint(measured[..., None] * weights, trajectory, (SIZE, SIZE)) gridded = bartorch.rss(channels[:, 0], axes=(0,)) .. GENERATED FROM PYTHON SOURCE LINES 327-338 Iterative SENSE --------------- Conjugate gradients on the normal equations :math:`A^H A x = A^H y`, with no penalty, is CG-SENSE [#pruessmann2001]_. The coil encoding separates the aliased signal the streaks consist of, but at this undersampling the problem is ill-conditioned, and each further iteration fits more of the noise; the number of iterations acts as the regularization. A total-variation penalty [#rof]_ adds prior knowledge instead: streaks and noise have a large total variation, the anatomy a small one. ADMM is the algorithm ``pics`` selects for this penalty. .. GENERATED FROM PYTHON SOURCE LINES 339-356 .. code-block:: Python cg_sense = apps.pics(measured, maps, traj=trajectory, maxiter=30) term = priors.TotalVariation(axes=(-1, -2), weight=TV_WEIGHT) start = time.perf_counter() reconstruction = apps.pics( measured, maps, traj=trajectory, regularizers=term, solver="admm", maxiter=ITERATIONS ) print(f"pics: {time.perf_counter() - start:.2f} s") results = {"gridding": gridded, "CG-SENSE": cg_sense, "SENSE + TV": reconstruction} for name, estimate in results.items(): error = bt.nrmse(image.abs(), estimate.abs(), scaled=True) similarity = bt.ssim(image.abs(), scaled(estimate, image)) print(f"{name:>12} NRMSE {error:.3f} SSIM {similarity:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none pics: 0.55 s gridding NRMSE 0.324 SSIM 0.418 CG-SENSE NRMSE 0.087 SSIM 0.566 SENSE + TV NRMSE 0.085 SSIM 0.858 .. GENERATED FROM PYTHON SOURCE LINES 357-384 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_002.png :alt: 48 golden-angle spokes, 8 coils, reference, gridding, CG-SENSE, SENSE + TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_003.png :alt: error magnitude, gridding, CG-SENSE, SENSE + TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_004.png :alt: enlarged, reference, gridding, CG-SENSE, SENSE + TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_02-radial-sense_004.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 385-406 Gridding shows the streaks of radial undersampling over the whole field of view, superimposed on an image that is otherwise sharp: the low spatial frequencies are fully sampled, and the streaks come from the high ones. The error map shows them extending outside the head, where the object has no signal. CG-SENSE removes the streaks, whose aliased signal the coil encoding separates, and leaves amplified noise across the head; its largest errors are at the scalp, whose bright, thin edge has most of its energy at spatial frequencies the sparse outer k-space samples poorly. The total-variation penalty suppresses the noise as well. Its residual error lies along the tissue boundaries, and in the enlarged region the thin cortical folds are flattened, since a boundary between two tissues of similar intensity also has a small total variation. The disc of k-space the trajectory samples limits the resolution of all three. The same solve through the operator ----------------------------------- Besides the iteration, ``pics`` scales the data. Off the Cartesian grid it estimates the scale from the adjoint reconstruction and therefore needs the operator, which :func:`bartorch.optim.data_scaling` takes. The encoding is the operator built above, now over the estimated sensitivities. .. GENERATED FROM PYTHON SOURCE LINES 407-418 .. code-block:: Python A = linop.NoncartesianSense(maps[:, 0], (SIZE, SIZE), traj=trajectory) data = measured / optim.data_scaling(measured[..., None], A=A) start = time.perf_counter() assembled = optim.ADMM(term, maxiter=ITERATIONS)(data, A) print(f"operator and solver: {time.perf_counter() - start:.2f} s") difference = (assembled.squeeze() - reconstruction.squeeze()).abs().max() print(f"relative difference from pics: {float(difference / reconstruction.abs().max()):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none operator and solver: 0.62 s relative difference from pics: 0.0e+00 .. GENERATED FROM PYTHON SOURCE LINES 419-434 The two run the same iteration over the same operator. The NUFFT spreads samples onto the grid over several threads and sums in the order they finish in, so the two can differ at the level of floating-point round-off. The normal operator ------------------- Each iteration applies :math:`A^H A`. For a single coil this is a convolution with the point spread function of the trajectory, so it can be computed exactly by FFTs on a grid of twice the matrix size (the Toeplitz embedding [#fessler]_) instead of by a NUFFT and an adjoint NUFFT; with coils it is that convolution between multiplications by the sensitivities. The operator uses the convolution by default, and ``toeplitz=False`` requests the transform pair. The two differ by the tolerance of the transforms, and the iterations carry that difference into the reconstructions. .. GENERATED FROM PYTHON SOURCE LINES 435-443 .. code-block:: Python start = time.perf_counter() pair = optim.ADMM(term, maxiter=ITERATIONS)( data, linop.NoncartesianSense(maps[:, 0], (SIZE, SIZE), traj=trajectory, toeplitz=False) ) print(f"without the Toeplitz normal: {time.perf_counter() - start:.2f} s") print(f"relative difference {float((pair - assembled).abs().max() / assembled.abs().max()):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none without the Toeplitz normal: 0.36 s relative difference 3.2e-02 .. GENERATED FROM PYTHON SOURCE LINES 444-451 The convolution costs an FFT, a pointwise multiplication and an inverse FFT on the doubled grid per coil, independent of the number of samples; the pair costs two non-uniform transforms, whose spreading and interpolation grow with the number of samples. With 48 spokes there are fewer samples than grid points, and the pair is not the slower of the two; as the number of samples grows, with more spokes or with the frames of a dynamic series sharing one normal operator, the convolution becomes the cheaper. .. GENERATED FROM PYTHON SOURCE LINES 454-479 References ---------- .. [#pruessmann2001] Pruessmann KP, Weiger M, Börnert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. *Magn Reson Med* 46(4):638-651 (2001). https://doi.org/10.1002/mrm.1241 .. [#winkelmann] Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. An optimal radial profile order based on the golden ratio for time-resolved MRI. *IEEE Trans Med Imaging* 26(1):68-76 (2007). https://doi.org/10.1109/TMI.2006.885337 .. [#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 .. [#rof] Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. *Physica D* 60(1-4):259-268 (1992). https://doi.org/10.1016/0167-2789(92)90242-F .. [#fessler] Fessler JA, Lee S, Olafsson VT, Shi HR, Noll DC. Toeplitz-based iterative image reconstruction for MRI with correction for magnetic field inhomogeneity. *IEEE Trans Signal Process* 53(9):3393-3402 (2005). https://doi.org/10.1109/TSP.2005.853152 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.361 seconds) .. _sphx_glr_download_auto_examples_04-non-cartesian_02-radial-sense.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-radial-sense.ipynb <02-radial-sense.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-radial-sense.py <02-radial-sense.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-radial-sense.zip <02-radial-sense.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_