.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/01-readout-oversampling-and-apodization.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_07-tours_01-readout-oversampling-and-apodization.py: ==================================== Readout oversampling and apodization ==================================== Two operations are applied to Cartesian k-space before the image is reconstructed. Readout oversampling doubles the field of view along the frequency-encoding direction, and is removed so that the image has the prescribed matrix. Truncation of k-space at the edge of the acquired matrix convolves the image with a sinc, whose side lobes appear as Gibbs ringing parallel to every sharp edge; an apodization window reduces the ringing at the cost of spatial resolution. This example removes twofold readout oversampling from the k-space of a Shepp-Logan phantom and compares the result with an acquisition without oversampling, then reconstructs a 64 x 64 acquisition with no window, a Fermi window and a Hann window, and measures the ringing and the resolution of each. The k-space is evaluated analytically at the sample positions, so its truncation and its oversampling are those of a continuous object. **Learning objectives** * Remove readout oversampling in the image domain with :func:`bartorch.remove_readout_oversampling`, and distinguish it from discarding the outer readout samples. * Relate the Gibbs ringing of a truncated acquisition to the side lobes of its point spread function. * Apodize k-space with :func:`bartorch.apodize`, and quantify the trade-off between ringing amplitude and the full width at half maximum. * Choose between the radial and the separable extension of a window over k-space. .. GENERATED FROM PYTHON SOURCE LINES 35-59 .. code-block:: Python import math import torch import bartorch import bartorch.tools as bt .. GENERATED FROM PYTHON SOURCE LINES 60-72 Readout oversampling -------------------- The receiver digitises the echo at twice the sampling rate the prescribed matrix requires, so that its anti-aliasing filter can roll off outside the field of view. Twofold oversampling halves the sample spacing :math:`\Delta k_x`, which doubles the field of view along the readout; the extent of k-space, and so the spatial resolution, is unchanged. The acquisition has 128 phase-encoding lines and 256 readout samples at half the phase-encoding spacing, in the grid units of a 128 matrix. The first trajectory component is :math:`k_x`, along the last image axis. .. GENERATED FROM PYTHON SOURCE LINES 73-82 .. code-block:: Python MATRIX = 128 phase_encodes = (torch.arange(MATRIX) - MATRIX // 2).float() readout = (torch.arange(2 * MATRIX) - MATRIX).float() / 2 ky, kx = torch.meshgrid(phase_encodes, readout, indexing="ij") trajectory = torch.stack([kx, ky, torch.zeros_like(kx)], dim=-1) oversampled = bt.phantom(traj=trajectory).reshape(MATRIX, 2 * MATRIX) .. GENERATED FROM PYTHON SOURCE LINES 83-91 :func:`bartorch.remove_readout_oversampling` transforms the readout to the image domain, crops the field of view to the prescribed matrix and transforms back. Discarding the outer half of the readout samples instead keeps the doubled field of view and halves the resolution along the readout. The crop is compared with an acquisition without oversampling; the crop uses the unitary transform, whose normalization depends on the number of samples, so the two differ by the factor :math:`\sqrt{2}` between the transform lengths. .. GENERATED FROM PYTHON SOURCE LINES 92-108 .. code-block:: Python cropped = bartorch.remove_readout_oversampling(oversampled, MATRIX, axis=-1) / math.sqrt(2) truncated = bartorch.resize(oversampled, (MATRIX, MATRIX)) / math.sqrt(2) reference = bt.phantom(MATRIX, kspace=True) def image_of(kspace): return bartorch.fft(kspace, axes=(-2, -1), inverse=True).abs() def nrmse(estimate, target): return float((estimate - target).norm() / target.norm()) print(f"crop against no oversampling: NRMSE {nrmse(image_of(cropped), image_of(reference)):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none crop against no oversampling: NRMSE 4.8e-03 .. GENERATED FROM PYTHON SOURCE LINES 109-126 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_001.png :alt: 256 readout samples :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_002.png :alt: no oversampling, image-domain crop, crop - no oversampling, k-space truncation :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_002.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 127-147 The image-domain crop reproduces the acquisition without oversampling to the part of the truncation ringing that extends past the prescribed field of view along the readout and is cropped with it. Discarding k-space samples instead widens the pixel to twice the prescribed size along :math:`x`: the object occupies half of the matrix. Apodization ----------- An acquisition truncated at the edge of k-space is the object's spectrum multiplied by a rectangle, and its image the object convolved with a sinc. The first side lobe of the sinc is 22 % of its peak, which appears as an overshoot of about 9 % at a step edge and as ringing that decays over a few pixels. An apodization window rolls the data off towards the edge of k-space, which lowers the side lobes and widens the main lobe of the point spread function [#bernstein]_. The phantom is acquired on a 64 matrix and reconstructed on a 128 grid by zero-filling, which interpolates the image and makes the ringing visible between the pixels of the acquired grid. .. GENERATED FROM PYTHON SOURCE LINES 148-158 .. code-block:: Python GRID, ACQUIRED = 128, 64 image = bt.phantom(GRID).abs() measured = bt.phantom(ACQUIRED, kspace=True) def reconstruct(kspace): return bartorch.fft(bartorch.resize(kspace, (GRID, GRID)), axes=(-2, -1), inverse=True).abs() .. GENERATED FROM PYTHON SOURCE LINES 159-166 :func:`bartorch.fermi_window` sets the radius of the half height and the width of the transition separately, and keeps a wide passband. :func:`bartorch.hann_window` tapers from the centre of k-space to zero at the edge. The point spread functions are evaluated on a grid eight times finer than the acquired one: the side lobes of an unwindowed acquisition have their zeros at the pixels of the acquired grid, so a point spread function read off that grid shows none. .. GENERATED FROM PYTHON SOURCE LINES 167-187 .. code-block:: Python fermi = bartorch.fermi_window((ACQUIRED, ACQUIRED), radius=0.8, width=0.08) hann = bartorch.hann_window((ACQUIRED, ACQUIRED)) windows = {"no window": torch.ones(ACQUIRED, ACQUIRED), "Fermi": fermi, "Hann": hann} UPSAMPLE = 8 offset = (torch.arange(UPSAMPLE * ACQUIRED) - UPSAMPLE * ACQUIRED // 2) / UPSAMPLE def psf(window): fine = bartorch.resize(window.to(torch.complex64), (UPSAMPLE * ACQUIRED,) * 2) profile = bartorch.fft(fine, axes=(-2, -1), inverse=True).abs()[UPSAMPLE * ACQUIRED // 2] return profile / profile.max() def fwhm(profile): above = torch.nonzero(profile > 0.5).flatten() return float(offset[above[-1]] - offset[above[0]]) .. GENERATED FROM PYTHON SOURCE LINES 188-193 :func:`bartorch.apodize` multiplies k-space by either window over the axes it names. The ringing is measured as the standard deviation of the image over the parenchyma within six pixels of the skull, where the object is uniform, and the resolution as the full width at half maximum of the point spread function, in pixels of the acquired grid. .. GENERATED FROM PYTHON SOURCE LINES 194-217 .. code-block:: Python reconstructions = { "no window": reconstruct(measured), "Fermi": reconstruct(bartorch.apodize(measured, kind="fermi", radius=0.8, width=0.08)), "Hann": reconstruct(bartorch.apodize(measured, kind="hann")), } skull = (image > 0.9).float()[None, None] near_skull = torch.nn.functional.max_pool2d(skull, 13, stride=1, padding=6)[0, 0] > 0 flat = torch.nn.functional.max_pool2d( (image - 0.2).abs().gt(1e-3).float()[None, None], 5, stride=1, padding=2 )[0, 0].eq(0) ringing_region = near_skull & flat print(f"{'window':>9} {'ringing':>7} {'side lobe':>9} {'FWHM':>7}") metrics = {} for name, window in windows.items(): profile = psf(window) ringing = float(reconstructions[name][ringing_region].std()) lobe = float(profile[offset.abs() > 3].max()) metrics[name] = (ringing, fwhm(profile)) print(f"{name:>9} {ringing:7.4f} {lobe:9.1e} {fwhm(profile):5.2f} px") .. rst-class:: sphx-glr-script-out .. code-block:: none window ringing side lobe FWHM no window 0.0404 9.1e-02 1.00 px Fermi 0.0325 2.6e-02 1.50 px Hann 0.0100 5.6e-03 2.00 px .. GENERATED FROM PYTHON SOURCE LINES 218-248 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_003.png :alt: object, no window, FWHM 1.00 px, Fermi, FWHM 1.50 px, Hann, FWHM 2.00 px :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_004.png :alt: |PSF|, row 64, across the skull :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_004.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 249-271 The zoomed panels show the left edge of the skull, displayed from zero to half the skull intensity. Without a window the ringing is visible as bands parallel to the skull across the adjacent parenchyma. The Fermi window lowers the far side lobes by a factor of three and the ringing next to the skull by about a fifth, and widens the point spread function by half a pixel; the Hann window removes the ringing almost entirely and doubles the full width at half maximum, which blurs the skull into the parenchyma. The Fermi window's wide passband is the usual compromise for anatomical imaging. Radial and separable windows ---------------------------- The one-dimensional kernel is extended over k-space either on the Euclidean norm of the normalized coordinates (``geometry="radial"``, an ellipse) or as a product along each axis (``geometry="separable"``), which retains more of the corners of k-space. Bernstein et al. [#bernstein]_ give the ratio of the two at the diagonal Nyquist point as 52.4 % in two dimensions and 50.7 % in three, for a Fermi window of transition width 10/128. The radial window is 0.5 there; the separable one is the product of the one-dimensional kernel at :math:`1/\sqrt{d}` along each of :math:`d` axes. The one-dimensional kernel at :math:`u` is the window at the centre of a grid with its radius moved to :math:`1 - u`. .. GENERATED FROM PYTHON SOURCE LINES 272-284 .. code-block:: Python TRANSITION = 10.0 / 128 def kernel(u): return float(bartorch.fermi_window((4, 4), radius=1.0 - u, width=TRANSITION)[2, 2]) for dimensions in (2, 3): separable = kernel(1 / math.sqrt(dimensions)) ** dimensions print(f"{dimensions}D: radial / separable at the diagonal {100 * 0.5 / separable:.1f} %") .. rst-class:: sphx-glr-script-out .. code-block:: none 2D: radial / separable at the diagonal 52.4 % 3D: radial / separable at the diagonal 50.7 % .. GENERATED FROM PYTHON SOURCE LINES 285-310 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_005.png :alt: radial, separable :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_006.png :alt: window along the diagonal :srcset: /auto_examples/07-tours/images/sphx_glr_01-readout-oversampling-and-apodization_006.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 311-323 The radial window has the more isotropic point spread function and the higher signal-to-noise ratio; the separable one the better resolution along the diagonals. Both normalize each axis to its own Nyquist edge, so a rectangular matrix receives an ellipse matched to its grid. References ---------- .. [#bernstein] Bernstein MA, Fain SB, Riederer SJ. Effect of windowing and zero-filled reconstruction of MRI data on spatial resolution and acquisition strategy. *J Magn Reson Imaging* 14(3):270-280 (2001). https://doi.org/10.1002/jmri.1183 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.814 seconds) .. _sphx_glr_download_auto_examples_07-tours_01-readout-oversampling-and-apodization.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-readout-oversampling-and-apodization.ipynb <01-readout-oversampling-and-apodization.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-readout-oversampling-and-apodization.py <01-readout-oversampling-and-apodization.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-readout-oversampling-and-apodization.zip <01-readout-oversampling-and-apodization.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_