.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/04-non-cartesian/01-trajectories-and-transforms.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_01-trajectories-and-transforms.py: ================================== 7. Basics of non-Cartesian imaging ================================== In lessons 2 to 6 every sample sat on the Cartesian grid, so the FFT and its inverse were all the Fourier transforms you needed. In this lesson you leave the grid. You sample the same brain slice along radial and spiral trajectories, reconstruct it by gridding, and see what each step of that reconstruction does to the image. The mental model: a radial or spiral readout starts at the centre of k-space and traces a curve through it, so its samples do not lie on a grid and the FFT cannot compute them. The non-uniform FFT (NUFFT) takes its place. Every non-Cartesian transform in bartorch is computed by FINUFFT [#finufft]_, which evaluates .. math:: y_j = \frac{1}{\sqrt{N}} \sum_{m} x_m \, \exp\!\Big(-2\pi i \sum_d \frac{k_{j,d}\, m_d}{n_d}\Big) to a requested tolerance. The sum runs over the :math:`N` voxels :math:`m` of an image with :math:`n_d` voxels along dimension :math:`d`, and :math:`k_j` is the position of sample :math:`j` in grid units. :doc:`../../explanation/non-cartesian` states the conventions and the accuracy. **Learning objectives** - Generate radial and golden-angle trajectories in grid units with :func:`bartorch.tools.traj`, and a spiral trajectory by hand. - Apply :func:`bartorch.nufft` and :func:`bartorch.nufft_adjoint`, and check the transform against the sum that defines it. - Explain why a gridding reconstruction needs density compensation, and compute the weights analytically and with :func:`bartorch.estimate_density`. - Apply the normal operator :math:`A^H A` as a convolution, and relate the point spread function (PSF) of an undersampled radial trajectory to its streaks. Previous: :doc:`../03-regularization/02-operators-and-solvers`. Next: :doc:`02-radial-sense`, where you remove the streaks of an undersampled radial scan with the coils and a penalty. .. GENERATED FROM PYTHON SOURCE LINES 47-137 .. GENERATED FROM PYTHON SOURCE LINES 138-159 Radial trajectory ----------------- The object is the BrainWeb slice of lesson 2 [#brainweb]_ on a 128 by 128 grid, acquired with a single receive coil of uniform sensitivity, so that the trajectory is the only encoding. A trajectory in bartorch is a tensor of shape ``(*encoding, shots, samples, 3)`` holding the k-space position of every sample in grid units: one unit is one k-space cell of the image, so a readout of ``SIZE`` samples runs from :math:`-n/2` to :math:`n/2` for an image of :math:`n` voxels. The last axis is :math:`(k_x, k_y, k_z)`; :math:`k_z` is zero throughout here, and that is what makes the transform two-dimensional. You can space the spokes of a radial scan in two ways. Uniformly, at :math:`\pi` over their number, they tile k-space evenly once all of them are acquired. At the golden angle, every next spoke falls into the largest gap left by the previous ones, so *any* run of consecutive spokes tiles k-space nearly evenly [#winkelmann]_. That is what lets you cut a golden-angle scan into frames after you have acquired it, as the Tour :doc:`../08-workflows/01-dynamic-golden-angle` does. .. GENERATED FROM PYTHON SOURCE LINES 159-192 .. code-block:: Python import time import torch import bartorch import bartorch.tools as bt from bartorch import linop SPOKES = 201 # pi/2 * SIZE: enough spokes not to undersample the edge uniform = bt.traj(readout=SIZE, spokes=SPOKES, radial=True) golden = bt.traj(readout=SIZE, spokes=SPOKES, radial=True, golden=True) print(f"trajectory {tuple(golden.shape)}: (spokes, samples, 3)") .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_001.png :alt: first 24 of 201 spokes, uniform, golden angle :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none trajectory (201, 128, 3): (spokes, samples, 3) .. GENERATED FROM PYTHON SOURCE LINES 193-206 You see the first 24 spokes of each scan. The uniform ones cover a narrow wedge, because they are spaced for 201 spokes and you would need all of them to fill the disc. The golden-angle ones already spread over every direction. The transform ------------- :func:`bartorch.nufft` samples an image along a trajectory, and :func:`bartorch.nufft_adjoint` takes samples back onto the grid. Since FINUFFT computes the transform only to a tolerance, check it once against the sum that defines it, evaluated by hand in double precision over the first spoke. That sum is computed without BART or FINUFFT, so it is an independent reference. .. GENERATED FROM PYTHON SOURCE LINES 206-219 .. code-block:: Python samples = bartorch.nufft(image, golden) print(f"samples {tuple(samples.shape)}") spoke = golden[0].real.double() # (samples, 3) m = torch.arange(SIZE, dtype=torch.float64) - SIZE // 2 my, mx = torch.meshgrid(m, m, indexing="ij") km = spoke[:, :1] * mx.reshape(1, -1) + spoke[:, 1:2] * my.reshape(1, -1) fourier = torch.exp(-2j * torch.pi / SIZE * km) explicit = fourier @ image.reshape(-1, 1).to(torch.complex128) / SIZE error = (explicit - samples[0]).abs().max() / explicit.abs().max() print(f"largest relative difference from the sum: {float(error):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none samples (201, 128, 1) largest relative difference from the sum: 1.9e-04 .. GENERATED FROM PYTHON SOURCE LINES 220-237 The difference sits within the tolerance the transform was planned with: a thousandth by default, on a grid a quarter larger than the image. That is meant for reconstruction, where the transform's error stays small beside the noise and the undersampling; :class:`bartorch.linop.NUFFT` takes ``oversampling`` and ``width`` when you need more accuracy. Density compensation -------------------- The adjoint is not the inverse. Every spoke passes through the centre of k-space, so the sample density is highest there and decreases outwards: the sampling density of a radial scan falls as :math:`1/\lvert k \rvert`. The adjoint sums the samples unweighted, so it overweights the low frequencies. To undo that you weight each sample by the inverse of the sampling density [#pipe]_, which for radial sampling is the ramp :math:`\lvert k \rvert`. The clamp keeps a small weight on the samples at the centre, which would otherwise not count at all. .. GENERATED FROM PYTHON SOURCE LINES 237-262 .. code-block:: Python def ramp(trajectory): """The radial density compensation, |k| per sample.""" k = torch.linalg.norm(trajectory.real[..., :2], dim=-1, keepdim=True) return k.clamp(min=0.25).to(torch.complex64) shape = (SIZE, SIZE) plain = bartorch.nufft_adjoint(samples, golden, image_shape=shape) compensated = bartorch.nufft_adjoint( samples * ramp(golden), golden, image_shape=shape, ) .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_002.png :alt: gridding reconstruction, 201 golden-angle spokes, the object, adjoint, no weights, ramp weights :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 263-279 Without weights the image is blurred: it is the object convolved with the PSF of the trajectory, which is broad because the sample density is concentrated at the centre of k-space. With the ramp weights the tissue boundaries are resolved. Neither image recovers what the trajectory never reaches: a radial scan covers a disc, so the corners of the Cartesian k-space are missing whatever the weights are. Estimated density ----------------- The ramp is the density of an ideal radial trajectory. For any other trajectory, :func:`bartorch.estimate_density` estimates the weights with the fixed-point iteration of Pipe and Menon [#pipe]_, which runs on the same NUFFTs. It takes the coordinates as ``(..., samples, 2)``, so you flatten the spokes into one list of samples and reshape the weights back. .. GENERATED FROM PYTHON SOURCE LINES 279-289 .. code-block:: Python flat = golden.real[..., :2].reshape(-1, 2) estimated = bartorch.estimate_density(flat, shape) estimated = estimated.reshape(SPOKES, SIZE, 1).to(torch.complex64) for name, density in (("ramp", ramp(golden)), ("Pipe-Menon", estimated)): weighted = samples * density adjoint = bartorch.nufft_adjoint(weighted, golden, image_shape=shape) nrmse = bt.nrmse(image.abs(), adjoint.abs(), scaled=True) print(f"{name:>10} NRMSE {nrmse:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none ramp NRMSE 0.297 Pipe-Menon NRMSE 0.278 .. GENERATED FROM PYTHON SOURCE LINES 290-301 For a radial trajectory the two weightings give similar errors, and both include the corners the disc does not cover. The estimate is needed where no closed form is available, as for a spiral. Spiral trajectory ----------------- :func:`bartorch.tools.traj` does not generate spirals, so you write one yourself: sixteen Archimedean interleaves reaching :math:`\pm n/2`, each with :math:`n/32` turns so that neighbouring turns of the whole set are one grid unit apart (the Nyquist spacing), and 1024 samples per interleaf. .. GENERATED FROM PYTHON SOURCE LINES 301-350 .. code-block:: Python INTERLEAVES, READOUT = 16, 1024 progress = torch.linspace(0.0, 1.0, READOUT) radius = SIZE / 2 * progress turns = SIZE / (2 * INTERLEAVES) arms = [] for arm in range(INTERLEAVES): angle = 2 * torch.pi * (turns * progress + arm / INTERLEAVES) kx, ky = radius * torch.cos(angle), radius * torch.sin(angle) arms.append(torch.stack([kx, ky, torch.zeros(READOUT)], dim=-1)) spiral = torch.stack(arms) # (interleaves, samples, 3) spiral_samples = bartorch.nufft(image, spiral) flat = spiral[..., :2].reshape(-1, 2) spiral_density = bartorch.estimate_density(flat, shape) spiral_density = spiral_density.reshape(INTERLEAVES, READOUT, 1) spiral_density = spiral_density.to(torch.complex64) spiral_plain = bartorch.nufft_adjoint(spiral_samples, spiral, image_shape=shape) spiral_compensated = bartorch.nufft_adjoint( spiral_samples * spiral_density, spiral, image_shape=shape ) spirals = {"no weights": spiral_plain, "Pipe-Menon": spiral_compensated} for name, adjoint in spirals.items(): nrmse = bt.nrmse(image.abs(), adjoint.abs(), scaled=True) print(f"spiral, {name:>10} NRMSE {nrmse:.3f}") .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_003.png :alt: 16 interleaves, the first in blue, Pipe-Menon weight along the first :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_004.png :alt: the object, spiral, no weights, spiral, Pipe-Menon :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none spiral, no weights NRMSE 0.693 spiral, Pipe-Menon NRMSE 0.089 .. GENERATED FROM PYTHON SOURCE LINES 351-370 The weight grows with :math:`\lvert k \rvert`: the interleaves stay a constant distance apart while each sample covers a longer stretch of arc, so the samples are densest at the centre. It ripples with the period of the turns, levels off in the outer part of k-space, and drops over the last samples, where the outermost turn has no neighbour outside it. Without the weights the spiral image is a low-contrast blur, as the radial one was; with them the tissue contrast and the edges are restored, and the NRMSE printed above drops accordingly. The normal operator ------------------- An iterative reconstruction applies :math:`A^H A` once per iteration, and for a NUFFT that costs a forward and an adjoint transform. Alternatively, :math:`A^H A` is a convolution of the image with the PSF, which is computed once on a grid twice the size of the image [#fessler2005]_ (the Toeplitz embedding). :class:`bartorch.linop.NUFFT` is the transform as an operator, and ``A.gram()`` applies its normal operator as that convolution. The timing below compares both forms: .. GENERATED FROM PYTHON SOURCE LINES 370-388 .. code-block:: Python A = linop.NUFFT(golden, image_shape=shape) repeats = 5 start = time.perf_counter() for _ in range(repeats): toeplitz = A.gram()(image) convolution = (time.perf_counter() - start) / repeats start = time.perf_counter() for _ in range(repeats): pair = A.H(A(image)) transforms = (time.perf_counter() - start) / repeats difference = (toeplitz - pair).abs().max() / pair.abs().max() print(f"A^H A as a convolution {1e3 * convolution:6.1f} ms") print(f"A^H A as two transforms {1e3 * transforms:6.1f} ms") print(f"relative difference {float(difference):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none A^H A as a convolution 1.0 ms A^H A as two transforms 0.9 ms relative difference 2.7e-03 .. GENERATED FROM PYTHON SOURCE LINES 389-399 The two agree to a small multiple of the transform's tolerance, and the convolution is the faster of the two. The point spread function ------------------------- :func:`bartorch.tools.psf` computes that PSF on its own: the adjoint reconstruction of a point object sampled along the trajectory. The reconstructed image is the object convolved with it. Compare the full scan with one of an eighth of the spokes: .. GENERATED FROM PYTHON SOURCE LINES 399-420 .. code-block:: Python full = bt.traj(readout=SIZE, spokes=SPOKES, radial=True, golden=True) eighth = bt.traj(readout=SIZE, spokes=SPOKES // 8, radial=True, golden=True) fully_sampled = bt.psf(full) undersampled = bt.psf(eighth) .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_005.png :alt: point spread function, normalized, 201 spokes, 25 spokes :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 421-431 On a logarithmic scale, the PSF of the full scan is a sharp central peak on a low, even skirt. With an eighth of the spokes the skirt breaks up into bright lines radiating from the centre. Convolving the object with this PSF produces the streaks of an undersampled radial image. Exercise: angular undersampling ------------------------------- Repeat the gridding reconstruction with the ramp weights for fewer and more spokes than 201, and compare the error and the images: .. GENERATED FROM PYTHON SOURCE LINES 431-463 .. code-block:: Python counts = (25, 50, 100, 201, 400) scores, images = [], {} for count in counts: trajectory = bt.traj(readout=SIZE, spokes=count, radial=True, golden=True) acquired = bartorch.nufft(image, trajectory) gridded = bartorch.nufft_adjoint( acquired * ramp(trajectory), trajectory, image_shape=shape, ) scores.append(bt.nrmse(image.abs(), gridded.abs(), scaled=True)) images[count] = gridded .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_006.png :alt: 01 trajectories and transforms :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_006.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_007.png :alt: 25 spokes, 50 spokes, 100 spokes :srcset: /auto_examples/04-non-cartesian/images/sphx_glr_01-trajectories-and-transforms_007.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 464-486 Below 201 spokes the error rises with every halving, and the images show why: streaks across the whole field of view, denser and brighter the fewer spokes you keep. Above 201 the error barely moves: undersampling no longer limits it, and what remains includes the corners of k-space the disc never covers. Try the same with uniform spokes (drop ``golden=True``), or with fewer spiral interleaves. As a spec --------- What this lesson built, stated the way you would ask an agent for it: .. code-block:: text With bartorch, sample a 128x128 complex image along a 2D golden-angle radial trajectory of 201 spokes of 128 samples (bt.traj, grid units), with bartorch.nufft, and check one spoke against the explicit DFT in float64. Reconstruct by gridding with bartorch.nufft_adjoint, without weights, with the ramp |k| (clamped at 0.25) and with bartorch.estimate_density, and report each NRMSE. Repeat for a 16-interleaf Archimedean spiral of 1024 samples. Time linop.NUFFT's A.gram() against A.H(A(x)), and compare bt.psf for 201 and 25 spokes. .. GENERATED FROM PYTHON SOURCE LINES 488-514 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 .. [#finufft] Barnett AH, Magland J, af Klinteberg L. A parallel nonuniform fast Fourier transform library based on an "exponential of semicircle" kernel. *SIAM J Sci Comput* 41(5):C479-C504 (2019). https://doi.org/10.1137/18M120885X .. [#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 .. [#pipe] Pipe JG, Menon P. Sampling density compensation in MRI: rationale and an iterative numerical solution. *Magn Reson Med* 41(1):179-186 (1999). https://doi.org/10.1002/(SICI)1522-2594(199901)41:1%3C179::AID-MRM25%3E3.0.CO;2-V .. [#fessler2005] 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 1.624 seconds) .. _sphx_glr_download_auto_examples_04-non-cartesian_01-trajectories-and-transforms.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-trajectories-and-transforms.ipynb <01-trajectories-and-transforms.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-trajectories-and-transforms.py <01-trajectories-and-transforms.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-trajectories-and-transforms.zip <01-trajectories-and-transforms.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_