.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/04-non-cartesian/03-dynamic-golden-angle.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_03-dynamic-golden-angle.py: =============================== Dynamic golden-angle radial MRI =============================== This lesson reconstructs a dynamic contrast-enhanced series from one continuous golden-angle radial acquisition, cut into frames of thirteen spokes each. Each frame on its own is undersampled fifteenfold and cannot be reconstructed; the series can, because consecutive frames are strongly correlated, and a total-variation penalty along the time axis states that correlation. The lesson compares frame-by-frame gridding with this joint reconstruction on the images and on the time-intensity curve a perfusion analysis would use. In a golden-angle acquisition [#winkelmann]_ each spoke is rotated from the previous one by :math:`180^\circ / \phi \approx 111.25^\circ`, with :math:`\phi` the golden ratio, so that any block of consecutive spokes covers k-space approximately uniformly, whatever its length and wherever it starts. The acquisition therefore runs without interruption, and the temporal resolution is chosen at reconstruction: fewer spokes per frame give a finer temporal resolution and stronger streak artefacts. Combined with parallel imaging and a sparsity penalty along time, this is GRASP [#feng]_. The encoding is that of :doc:`02-radial-sense` with a frame axis added: the image is ``(frames, y, x)``, the trajectory indexes frames as well as spokes, and the coil sensitivities are shared by all frames. 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** - Divide a continuous golden-angle acquisition into frames after the fact. - Build an encoding whose image and trajectory carry a frame axis. - Regularize along time with a total-variation term over the frame axis. - Compare frame-by-frame gridding with the joint reconstruction in the images, in an x-t profile and in the time-intensity curve of a region. It follows :doc:`02-radial-sense`. The next section begins with :doc:`../05-model-based/01-subspace-t1-mapping`, which constrains the time axis by a signal model. .. GENERATED FROM PYTHON SOURCE LINES 45-160 .. 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 linop, optim, priors SIZE = 128 COILS = 8 FRAMES = 16 SPOKES = 13 # per frame .. GENERATED FROM PYTHON SOURCE LINES 161-171 A contrast-enhanced series -------------------------- The phantom is the BrainWeb slice of the previous lessons with a contrast agent bolus passing through it. A gamma-variate curve describes the first-pass concentration over time, and each tissue enhances in proportion to its blood volume: strongly in grey matter, weakly in white matter, and not at all in cerebrospinal fluid. The series is therefore smooth in time, with the same anatomy in every frame, which is the structure the temporal penalty exploits. .. GENERATED FROM PYTHON SOURCE LINES 172-260 .. code-block:: Python CLASSES = {"grey matter": ("GM", 0.8), "white matter": ("WM", 0.25)} .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_001.png :alt: phantom: before, during and after the first pass, frame 0, frame 4, frame 15 :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 261-272 Acquisition ----------- ``FRAMES * SPOKES`` spokes are generated as one golden-angle trajectory and reshaped, so that the first axis indexes frames and the second the spokes within a frame. With :math:`\pi/2 \times 128 \approx 201` spokes needed for a fully sampled frame, thirteen spokes undersample each frame by a factor of about fifteen. The image varies along the frame axis, so the operator holds one NUFFT per frame, all applied inside the same loop over the coils. The noise is that of the Cartesian lessons, of variance :math:`10^{-4}` per sample. .. GENERATED FROM PYTHON SOURCE LINES 273-290 .. code-block:: Python trajectory = bt.traj(readout=SIZE, spokes=FRAMES * SPOKES, radial=True, golden=True) trajectory = trajectory.reshape(FRAMES, SPOKES, SIZE, 3) A = linop.NoncartesianSense(sensitivities, (FRAMES, SIZE, SIZE), traj=trajectory) measured = bt.noise(A(series), n=1e-4, s=3) print(f"{A.ishape} -> {A.oshape}") print(A.plan) .. rst-class:: sphx-glr-script-out .. code-block:: none (16, 128, 128) -> (8, 16, 13, 128) Plan(transform=nufft, items=16, image=sensitivities, normal=kernel, coil_batch=1, streamed=coils, executor=slab) .. GENERATED FROM PYTHON SOURCE LINES 291-306 ``plan.items`` is the number of frames the encoding carries, each with its own transform. Reconstruction -------------- Two reconstructions of the same data. The first treats the frames as independent: the adjoint of the encoding applied to density-compensated samples, which is the gridding reconstruction of thirteen spokes per frame, with the coils combined by the sensitivities. The second solves for the whole series at once, with a total-variation penalty along the frame axis: the solution is the series that explains all the data and changes least from frame to frame. The streak pattern of each frame is different, because each frame has different spokes, so it has a large temporal total variation and is suppressed, while the anatomy, which is the same in every frame, is not. .. GENERATED FROM PYTHON SOURCE LINES 307-317 .. code-block:: Python weights = torch.linalg.norm(trajectory.real[..., :2], dim=-1).clamp(min=0.25) gridded = A.H(measured * weights.to(torch.complex64)) data = measured / optim.data_scaling(measured[..., None], A=A) temporal = optim.ADMM(priors.TotalVariation(axes=(-3,), weight=0.02), maxiter=30)(data, A) for name, volume in (("gridding", gridded), ("temporal TV", temporal)): print(f"{name:>12} NRMSE over all frames {bt.nrmse(series.abs(), volume, scaled=True):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none gridding NRMSE over all frames 0.449 temporal TV NRMSE over all frames 0.126 .. GENERATED FROM PYTHON SOURCE LINES 318-322 The frame axis is ``-3``, the axis in front of the two spatial ones. A term given ``(-1, -2)`` would penalize the spatial gradient instead, and one given all three would penalize both; the axes a term acts on are the whole difference between a spatial and a temporal regularizer. .. GENERATED FROM PYTHON SOURCE LINES 323-343 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_002.png :alt: frame 4, peak of the first pass, 13 spokes, reference, gridding, temporal TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_003.png :alt: error magnitude, frame 4, gridding, temporal TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 344-353 In the gridding reconstruction the streaks of thirteen spokes dominate the frame, and only the ventricles and the outline of the head are recognizable. The joint reconstruction recovers the anatomy and the enhanced cortex; its residual error is small and concentrated at the tissue boundaries, which carry the high spatial frequencies each frame samples most sparsely. An x-t profile, one line of the image plotted against time, shows the time axis directly. The line below runs left to right through the ventricles and the grey matter on either side. .. GENERATED FROM PYTHON SOURCE LINES 354-374 .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_004.png :alt: x-t profile, row 58, reference, gridding, temporal TV :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 375-388 In the reference profile the grey matter brightens and fades over a few frames, while cerebrospinal fluid and the scalp stay constant. Gridding shows the same enhancement under a different streak pattern in every frame, so the profile changes from one row to the next even where the object does not; temporal total variation removes that variation and keeps the time course. Time-intensity curve -------------------- A perfusion study reports the signal in a region as a function of time, so the reconstructions are compared on that curve too. The region is the grey matter, where the enhancement is strongest. .. GENERATED FROM PYTHON SOURCE LINES 389-416 .. code-block:: Python region = memberships[CLASS["GM"]] > 0.6 curves = { "reference": series.abs(), "gridding": scaled(gridded, series), "temporal TV": scaled(temporal, series), } truth = curves["reference"][:, region].mean(-1) for name, volume in curves.items(): if name == "reference": continue enhancement = volume[:, region].mean(-1) print(f"{name:>12} curve NRMSE {float((enhancement - truth).norm() / truth.norm()):.3f}") .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_005.png :alt: grey matter :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_03-dynamic-golden-angle_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none gridding curve NRMSE 0.153 temporal TV curve NRMSE 0.039 .. GENERATED FROM PYTHON SOURCE LINES 417-425 Averaged over the grey matter, the streaks largely cancel, and gridding recovers the shape of the curve but not its level: part of the signal of each frame is spread into streaks across the field of view, outside the region. The joint reconstruction follows the reference curve closely, with the peak slightly attenuated and the baseline slightly raised: a temporal total-variation penalty flattens a signal change that lasts only a few frames more than any other feature, and a larger weight trades more of the peak for less noise. .. GENERATED FROM PYTHON SOURCE LINES 428-441 References ---------- .. [#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 .. [#feng] Feng L, Grimm R, Block KT, Chandarana H, Kim S, Xu J, Axel L, Sodickson DK, Otazo R. Golden-angle radial sparse parallel MRI: combination of compressed sensing, parallel imaging, and golden-angle radial sampling for fast and flexible dynamic volumetric MRI. *Magn Reson Med* 72(3):707-717 (2014). https://doi.org/10.1002/mrm.24980 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.086 seconds) .. _sphx_glr_download_auto_examples_04-non-cartesian_03-dynamic-golden-angle.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-dynamic-golden-angle.ipynb <03-dynamic-golden-angle.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-dynamic-golden-angle.py <03-dynamic-golden-angle.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-dynamic-golden-angle.zip <03-dynamic-golden-angle.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_