.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/04-gradient-nonlinearity.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_04-gradient-nonlinearity.py: ===================== Gradient nonlinearity ===================== Spatial encoding assumes that every gradient field changes linearly with position. A real gradient coil departs from that more and more with the distance from isocentre, so each spin is encoded at a position displaced from its true one. The image is warped, by a few millimetres at the edge of a head-sized field of view and by centimetres at the edge of a body-sized one, and the voxel volume changes with the warp. The displacement is a property of the coil: the manufacturer states it as the coefficients of a spherical-harmonic expansion of each gradient field [#janke]_. To correct it you evaluate that expansion at every voxel and resample the image at the positions where the voxels were actually encoded. In this tour you warp a grid phantom in a coronal slice over a 450 mm field of view with a coil described by third-order harmonics, and correct it with :class:`bartorch.tools.Gradunwarp`. The coil is a generic one, defined in the code; no manufacturer's table is used. Prerequisites: lesson 1, :doc:`../01-basics/01-tensors-and-commands`. **Learning objectives** - Describe a gradient coil's nonlinearity by its spherical-harmonic coefficients with :class:`~bartorch.tools.GradientCoefficients`. - Relate the displacement to the distance from isocentre and to the gradient axis. - Correct the geometry and the intensity of an image with :class:`~bartorch.tools.Gradunwarp`, and tell the two apart. - Place an image in scanner coordinates by its orientation and field of view. .. GENERATED FROM PYTHON SOURCE LINES 37-69 .. code-block:: Python import numpy as np import torch import bartorch.tools as bt .. GENERATED FROM PYTHON SOURCE LINES 70-88 The coil -------- :class:`~bartorch.tools.GradientCoefficients` holds the cosine and sine coefficients :math:`\alpha_{nm}` and :math:`\beta_{nm}` of each gradient field's departure from linearity, in the Siemens convention: harmonics normalized over a reference radius :math:`R_0`, and positions in scanner coordinates, with :math:`z` along the bore. An all-zero table is a perfectly linear coil. Each gradient field is odd along its own axis, so its lowest-order nonlinear terms are of third order: :math:`\alpha_{31}` for the :math:`x` gradient, :math:`\beta_{31}` for the :math:`y` gradient and :math:`\alpha_{30}` for the :math:`z` gradient, whose harmonic is proportional to :math:`z\,(2z^2 - 3x^2 - 3y^2)`. The signs and sizes of the terms decide where the image is compressed and where it is stretched. On a short whole-body coil the :math:`z` gradient is usually the least linear, so you give it the largest term: .. GENERATED FROM PYTHON SOURCE LINES 88-98 .. code-block:: Python ORDER = 3 alpha = np.zeros((3, ORDER + 1, ORDER + 1)) beta = np.zeros_like(alpha) alpha[0, 3, 1] = -0.04 # x gradient beta[1, 3, 1] = -0.04 # y gradient alpha[2, 3, 0] = -0.06 # z gradient coil = bt.GradientCoefficients( basis="normalized", alpha=alpha, beta=beta, reference_radius_mm=250.0 ) .. GENERATED FROM PYTHON SOURCE LINES 99-109 The slice --------- The correction needs to know where each voxel sits in the scanner. Your slice is coronal through isocentre, 256 x 256 over a 450 mm field of view (1.8 mm voxels): rows run from superior to inferior along :math:`-z`, and columns from right to left along :math:`x`. The orientation matrix says, for each array axis, along which scanner direction it increases: its columns are the row axis and the column axis, its rows are :math:`x`, :math:`y` and :math:`z`. .. GENERATED FROM PYTHON SOURCE LINES 109-119 .. code-block:: Python SIZE, FOV_MM = 256, 450.0 VOXEL_MM = FOV_MM / SIZE coronal = np.array([[0.0, 1.0], [0.0, 0.0], [-1.0, 0.0]]) unwarp = bt.Gradunwarp( coil, shape=(SIZE, SIZE), fov_mm=(FOV_MM, FOV_MM), orientation=coronal, ) .. GENERATED FROM PYTHON SOURCE LINES 120-124 :class:`~bartorch.tools.Gradunwarp` evaluates the expansion at every voxel of the corrected grid. Its ``source_grid`` is the index into the acquired image at which each voxel was encoded, so the displacement is the difference from the voxel's own index: .. GENERATED FROM PYTHON SOURCE LINES 124-137 .. code-block:: Python grid = np.meshgrid(np.arange(SIZE), np.arange(SIZE), indexing="ij") index = np.stack(grid, axis=-1) displacement_mm = (unwarp.source_grid - index) * VOXEL_MM offset_mm = (index - (SIZE - 1) / 2) * VOXEL_MM .. rst-class:: sphx-glr-script-out .. code-block:: none 100 mm from isocentre: displacement 0.9 mm along z, 0.5 mm along x 150 mm from isocentre: displacement 3.3 mm along z, 1.8 mm along x 200 mm from isocentre: displacement 7.6 mm along z, 4.1 mm along x .. GENERATED FROM PYTHON SOURCE LINES 138-154 The displacement grows with the cube of the distance from isocentre: under a millimetre at 100 mm, several at 200 mm, and about twice as large along :math:`z` as along :math:`x`, as the coefficients say it should be. The acquisition --------------- The object is a grid phantom: a disc 420 mm across carrying lines 30 mm apart, written in closed form so that you can evaluate it at any position. The acquired image at index :math:`p` holds the object at the position :math:`r` encoded there, the solution of :math:`r + d(r) = p`, found by fixed-point iteration. Its intensity is divided by the Jacobian determinant of the mapping, because a voxel whose volume the nonlinearity enlarges collects the signal of a larger region. That code simulates the artefact rather than correcting it, so it is hidden; ``acquired`` is the warped image and ``truth`` the object on the true grid. .. GENERATED FROM PYTHON SOURCE LINES 154-225 .. image-sg:: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_001.png :alt: acquired, displacement, Jacobian determinant :srcset: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 226-242 The centre of the field of view is untouched, and the displacement rises towards the edge, fastest at the top and bottom, along :math:`z`. There the Jacobian determinant is below one: the image is compressed towards isocentre, each voxel collects the signal of a larger region, and the periphery comes out brighter. Along :math:`x` the image is stretched and darker. Off the axes the :math:`z` displacement follows its harmonic and changes sign where :math:`2z^2 = 3x^2`, which is the cross you see in the Jacobian determinant. The correction -------------- Calling the :class:`~bartorch.tools.Gradunwarp` on an image resamples it at the source grid by cubic B-spline interpolation and multiplies it by the Jacobian determinant, which restores the intensity. With ``jacobian=False`` it corrects the geometry only, which lets you see the two parts separately: .. GENERATED FROM PYTHON SOURCE LINES 242-317 .. code-block:: Python image = torch.as_tensor(acquired, dtype=torch.float32) corrected = unwarp(image).numpy() geometry = bt.Gradunwarp( coil, shape=(SIZE, SIZE), fov_mm=(FOV_MM, FOV_MM), orientation=coronal, jacobian=False, ) geometry_only = geometry(image).numpy() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_002.png :alt: orange: grid of the object, acquired, corrected :srcset: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_003.png :alt: superior right corner (blue box), object, acquired, geometry only, corrected :srcset: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_004.png :alt: geometry only - object, corrected - object :srcset: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_005.png :alt: along z, along x :srcset: /auto_examples/07-tours/images/sphx_glr_04-gradient-nonlinearity_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none acquired NRMSE 0.427 geometry only NRMSE 0.118 corrected NRMSE 0.025 .. GENERATED FROM PYTHON SOURCE LINES 318-345 In the acquired image the grid lines bend away from the grid of the object (orange) towards the edge; after correction they lie on it. In the enlarged corner the geometry-only correction already puts the lines back in place, but leaves the intensity error of the Jacobian: its difference from the object is positive (red) at the top and bottom and negative (blue) at the sides. The full correction removes both, and the NRMSE against the object drops from the acquired image's to a few per cent. What is left is the error of interpolating the acquired image, largest on the thin lines where compression has undersampled them. The profiles show the same thing on the two axes: the acquired lines move inwards and brighten along :math:`z`, move outwards along :math:`x`, and the corrected lines sit on the object's. For a real coil, :meth:`~bartorch.tools.GradientCoefficients.from_file` reads the manufacturer's table, and :meth:`~bartorch.tools.Gradunwarp.from_mrd` and :meth:`~bartorch.tools.Gradunwarp.from_affine` take the geometry of the image from an MRD header or an affine. Scanners apply a two-dimensional correction in the plane of the slice by default; the through-plane displacement moves signal between slices, and correcting it needs the whole volume. Exercise: coil nonlinearity --------------------------- The displacement is linear in the coefficients. Scale the :math:`z` gradient's term and see how far a voxel 200 mm above isocentre moves: .. GENERATED FROM PYTHON SOURCE LINES 345-363 .. code-block:: Python for term in (-0.03, -0.06, -0.12): stronger = alpha.copy() stronger[2, 3, 0] = term trial = bt.Gradunwarp( bt.GradientCoefficients( basis="normalized", alpha=stronger, beta=beta, reference_radius_mm=250.0, ), shape=(SIZE, SIZE), fov_mm=(FOV_MM, FOV_MM), orientation=coronal, ) moved_mm = (trial.source_grid - index) * VOXEL_MM along_z, _ = displacement_at(moved_mm, offset_mm, 200) print(f"alpha_30 = {term:+.2f}: {along_z:4.1f} mm at z = 200 mm") .. rst-class:: sphx-glr-script-out .. code-block:: none alpha_30 = -0.03: 3.8 mm at z = 200 mm alpha_30 = -0.06: 7.6 mm at z = 200 mm alpha_30 = -0.12: 15.2 mm at z = 200 mm .. GENERATED FROM PYTHON SOURCE LINES 364-382 Doubling the coefficient doubles the displacement. Try a positive sign as well, and look at which way the grid lines move. As a spec --------- What this tour built, stated the way you would ask an agent for it: .. code-block:: text Correct the gradient-nonlinearity distortion of a 2D coronal slice (256 x 256, 450 mm field of view, rows superior to inferior along -z, columns along +x) with bt.Gradunwarp. Describe the coil with bt.GradientCoefficients (normalized basis, R0 250 mm, third-order alpha_31 for x, beta_31 for y, alpha_30 for z). Apply the full correction and a geometry-only one (jacobian=False), report the NRMSE of each against the object, and the displacement in mm at 100, 150 and 200 mm from isocentre along z and x. .. GENERATED FROM PYTHON SOURCE LINES 384-391 References ---------- .. [#janke] Janke A, Zhao H, Cowin GJ, Galloway GJ, Doddrell DM. Use of spherical harmonic deconvolution methods to compensate for nonlinear gradient effects on MRI images. *Magn Reson Med* 52(1):115-122 (2004). https://doi.org/10.1002/mrm.20122 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.533 seconds) .. _sphx_glr_download_auto_examples_07-tours_04-gradient-nonlinearity.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04-gradient-nonlinearity.ipynb <04-gradient-nonlinearity.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04-gradient-nonlinearity.py <04-gradient-nonlinearity.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04-gradient-nonlinearity.zip <04-gradient-nonlinearity.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_