.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/07-epi-susceptibility-distortion.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_07-epi-susceptibility-distortion.py: ================================ Susceptibility distortion in EPI ================================ In this tour you correct the susceptibility distortion of a measured EPI pair. In an echo-planar image, the phase-encoding direction is sampled once per echo spacing rather than once per dwell time, so its bandwidth per pixel is only a few tens of hertz. A spin off resonance by :math:`\Delta f` is displaced along phase encoding by :math:`\Delta f` divided by that bandwidth, which near the frontal sinus and the petrous bone comes to several millimetres at 3 T. The orbitofrontal cortex and the temporal poles are compressed or stretched, and signal piles up where neighbouring voxels are displaced onto the same place. The displacement changes sign with the direction in which you traverse k-space. Two acquisitions with opposite phase-encoding polarity, *blip-up* and *blip-down*, are distorted in opposite directions, and the displacement field that brings them into register is the correction [#andersson]_. You correct such a pair, measured at 3 T and published on OpenNeuro [#ds001600]_, with :func:`bartorch.tools.correct_susceptibility`, which runs PyHySCO [#pyhysco]_, and you check the estimated displacement against a gradient-echo field map of the same subject. Prerequisites: lesson 1, :doc:`../01-basics/01-tensors-and-commands`. **Learning objectives** - Read a BIDS EPI series and its sidecar, and compute the bandwidth per pixel along phase encoding and the displacement per hertz of off-resonance. - Tell the direction of the displacement in anatomical terms from the phase-encoding direction and the image orientation. - Recognise the compression, stretching and pile-up of susceptibility distortion, and their reversal between the two polarities. - Estimate the displacement field from a reversed phase-encoding pair and correct both images, intensity included. - Check the estimate against an independent field map, and the correction by the agreement of the two corrected images. .. GENERATED FROM PYTHON SOURCE LINES 43-129 .. GENERATED FROM PYTHON SOURCE LINES 130-142 The data -------- The dataset holds gradient-echo EPI series of one subject acquired on a Siemens Prisma at 3 T: a 64 x 64 matrix, 44 axial slices, 3.75 x 3.75 x 4 mm, with the same protocol and shim and opposite phase-encoding polarity, and a dual-echo gradient-echo field map on the same grid. Each EPI series has five volumes, and you average them to raise the signal-to-noise ratio. :func:`~bartorch.io.read_nifti` returns the volumes as ``(volumes, z, y, x)`` together with the affine from voxel indices ``(x, y, z)`` to RAS millimetres. The JSON sidecars hold the timings of the readout. The download into ``CACHE`` is hidden. .. GENERATED FROM PYTHON SOURCE LINES 142-172 .. code-block:: Python import json import torch import bartorch.tools as bt from bartorch.io import read_nifti from nibabel.orientations import aff2axcodes blip_up_file = CACHE / "sub-1/fmap/sub-1_dir-PA_epi" blip_down_file = CACHE / "sub-1/func/sub-1_task-rest_acq-AP_bold" blip_up = read_nifti(f"{blip_up_file}.nii.gz") blip_down = read_nifti(f"{blip_down_file}.nii.gz") sidecar_up = json.loads(Path(f"{blip_up_file}.json").read_text()) sidecar_down = json.loads(Path(f"{blip_down_file}.json").read_text()) affine = blip_up.affine spacing = affine[:3, :3].norm(dim=0).flip(0) voxel_mm = tuple(float(v) for v in spacing) # (z, y, x) up = blip_up.image.mean(0) down = blip_down.image.mean(0) .. rst-class:: sphx-glr-script-out .. code-block:: none (5, 44, 64, 64) volumes, voxel (z, y, x) (4.0, 3.75, 3.75) mm voxel axes x, y, z point to LAS blip-up BOLD_3mm_matched_0.55esp_topup: j blip-down BOLD_3mm_matched_0.55esp: j- .. GENERATED FROM PYTHON SOURCE LINES 173-198 The voxel axis ``x`` points to the subject's left and ``y`` to anterior, so ``PhaseEncodingDirection`` ``j`` is posterior to anterior and ``j-`` anterior to posterior. Phase encoding is along the anterior-posterior axis of the head, and the readout is left-right. The figures show the slices with anterior at the top and the subject's left on the right. Displacement along phase encoding --------------------------------- You acquire :math:`N` phase-encoding lines at the effective echo spacing :math:`\Delta t_{esp}`, the echo spacing divided by any parallel-imaging acceleration. Off resonance, the phase you accrue is linear in :math:`k_y`, and a linear phase in k-space is a shift in the image: .. math:: d = \frac{\Delta f}{\mathrm{BW}_{PE}}\ \text{voxels}, \qquad \mathrm{BW}_{PE} = \frac{1}{N\, \Delta t_{esp}}, with :math:`\mathrm{BW}_{PE}` the bandwidth per pixel along phase encoding. BIDS writes :math:`(N - 1)\,\Delta t_{esp}` as ``TotalReadoutTime``. For a positive :math:`\Delta f` the displacement is along ``PhaseEncodingDirection``, and reversing the blips reverses it. The readout, at a bandwidth per pixel of ``PixelBandwidth``, is displaced by a fraction of a voxel, and you neglect it. .. GENERATED FROM PYTHON SOURCE LINES 198-216 .. code-block:: Python lines = sidecar_up["AcquisitionMatrixPE"] echo_spacing = sidecar_up["EffectiveEchoSpacing"] bandwidth_pe = 1 / (lines * echo_spacing) .. rst-class:: sphx-glr-script-out .. code-block:: none 64 lines at an effective echo spacing of 0.55 ms bandwidth per pixel: 28.4 Hz along phase encoding (sidecar 28.409 Hz), 2055 Hz along the readout total readout time 34.65 ms 100 Hz off resonance: 3.5 voxels = 13.2 mm along phase encoding .. GENERATED FROM PYTHON SOURCE LINES 217-232 The bandwidth you compute agrees with the one the sidecar states. Along the readout the bandwidth per pixel is about seventy times larger, which is why only phase encoding is distorted. The reference field map ----------------------- To check the correction you need an independent measurement of the field. The gradient-echo field map is the phase difference :math:`\Delta\phi` between two echoes, stored by the scanner as integers from 0 to 4095 for :math:`-\pi` to :math:`\pi`, and the off-resonance is :math:`\Delta f = \Delta\phi / (2\pi\, \Delta\mathrm{TE})`. You smooth the phase over one voxel (as a complex exponential, hidden below). The displacement it predicts for the blip-up image is :math:`\Delta f / \mathrm{BW}_{PE}` voxels towards anterior. .. GENERATED FROM PYTHON SOURCE LINES 232-262 .. code-block:: Python phase_file = CACHE / "sub-1/fmap/sub-1_acq-v4_phasediff" magnitude_file = CACHE / "sub-1/fmap/sub-1_acq-v4_magnitude1" sidecar_phase = json.loads(Path(f"{phase_file}.json").read_text()) stored = read_nifti(f"{phase_file}.nii.gz").image[0] phase = (stored - 2048) / 2048 * torch.pi magnitude = read_nifti(f"{magnitude_file}.nii.gz").image[0] delta_te = sidecar_phase["EchoTime2"] - sidecar_phase["EchoTime1"] off_resonance = phase / (2 * torch.pi * delta_te) predicted_mm = off_resonance / bandwidth_pe * voxel_mm[1] .. rst-class:: sphx-glr-script-out .. code-block:: none echo times 4.00 and 6.46 ms: unambiguous within +/-203 Hz .. GENERATED FROM PYTHON SOURCE LINES 263-280 Reversed-gradient correction ---------------------------- PyHySCO looks for the displacement field :math:`b` for which the blip-up image sampled at :math:`y + b` and the blip-down image sampled at :math:`y - b`, each multiplied by the Jacobian determinant of its transformation, :math:`1 \pm \partial b / \partial y`, agree. The Jacobian factor restores the intensity of voxels compressed into a pile-up or stretched over several voxels. A smoothness penalty on :math:`b`, and a constraint that keeps both transformations invertible, regularise the problem. The estimation is three-dimensional: you pass the volumes with their voxel size and the phase-encoding axis ``y``. You get :math:`b` back in millimetres along the voxel axis ``y``, that is towards anterior, on the faces between voxels, and average neighbouring faces to put it at the voxel centres. .. GENERATED FROM PYTHON SOURCE LINES 280-303 .. code-block:: Python result = bt.correct_susceptibility( up, down, voxel_size=voxel_mm, phase_encoding_axis=1, ) faces = result.field_map estimated_mm = (0.5 * (faces[:, 1:] + faces[:, :-1])).float() corrected_up = result.blip_up.float() corrected_down = result.blip_down.float() .. GENERATED FROM PYTHON SOURCE LINES 304-311 The corrected pair ------------------ If the correction is right, the two corrected images are the same image. You measure their agreement over the brain, where the corrected EPI and the field map's magnitude both have signal, as the correlation coefficient and as the root-mean-square difference relative to the mean image. .. GENERATED FROM PYTHON SOURCE LINES 311-326 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_001.png :alt: blip-up, P to A, blip-down, A to P, corrected :srcset: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_002.png :alt: blip-up, P to A, blip-down, A to P, corrected :srcset: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none blip-up vs blip-down, acquired : correlation 0.871, RMS difference 21 % of the mean image blip-up vs blip-down, corrected: correlation 0.989, RMS difference 6 % of the mean image .. GENERATED FROM PYTHON SOURCE LINES 327-340 Slice 21 passes through the orbitofrontal cortex and slice 14 through the temporal poles and the cerebellum. The lower row of each figure enlarges the anterior part of the slice, in the box, with its own grey scale. Above the frontal sinus the off-resonance is positive. The blip-up image, phase-encoded from posterior to anterior, displaces the orbitofrontal cortex anteriorly and compresses its anterior edge into a bright rim; the blip-down image displaces it posteriorly and flattens it. At the temporal poles, above the petrous bone and the mastoid air cells, the off-resonance is negative and the directions swap: the blip-down image stretches the poles anteriorly into streaks, and the blip-up image compresses them into a few dark voxels. The corrected image, the mean of the two corrected images, puts the cortex between the two and restores its intensity. .. GENERATED FROM PYTHON SOURCE LINES 340-353 .. image-sg:: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_003.png :alt: acquired, slice 21, corrected, slice 21, acquired, slice 14, corrected, slice 14 :srcset: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 354-373 In these overlays the blip-up image is magenta and the blip-down image green, so tissue where the two agree is grey. Before correction the frontal, temporal and occipital edges carry a magenta fringe on one side and a green fringe on the other; after correction they are grey. What is left is green at the temporal poles, where the blip-up image has lost the signal that the blip-down image still has: signal lost within a voxel in one acquisition cannot be brought back by a displacement. Comparison with the field map ----------------------------- Now you compare the estimate with the displacement the gradient-echo field map predicts. You do it over the whole brain, and over the brain voxels where the field map predicts a displacement of more than one voxel, :math:`|\Delta f| > \mathrm{BW}_{PE}`: that is where there is a distortion to measure. Over the rest of the brain the field is close to zero, and a correlation there measures mostly noise. The field map was acquired with its own shim, so the two fields may also differ by a smooth field of first and second order. .. GENERATED FROM PYTHON SOURCE LINES 373-434 .. code-block:: Python distorted = brain & (off_resonance.abs() > bandwidth_pe) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_004.png :alt: PyHySCO estimate, slice 21, GRE field map, slice 21, PyHySCO estimate, slice 14, GRE field map, slice 14 :srcset: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_005.png :alt: brain, displaced > 1 voxel :srcset: /auto_examples/07-tours/images/sphx_glr_07-epi-susceptibility-distortion_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none shim, EPI [1746, -3889, -287, 575, 111, -87, -1, -28] shim, field map [1747, -3883, -285, 583, 125, -23, -45, -40] field over the brain: -45 to +47 Hz (1st to 99th percentile) displaced by more than one voxel: 7 % of the brain voxels estimate vs field map over the brain: correlation 0.38, slope 0.31 estimate vs field map over the displaced voxels: correlation 0.61, slope 0.30 .. GENERATED FROM PYTHON SOURCE LINES 435-469 The estimate follows the pattern of the field-map prediction: anterior displacement above the frontal sinus and in the cerebellum, posterior displacement at the temporal poles. It is smaller in amplitude: the fitted slope is about 0.3, and in the joint histograms the voxels lie flatter than the diagonal. Two likely causes, neither established here, are the smoothness penalty on the displacement field, which at 3.75 mm voxels flattens a displacement that changes within a few voxels, and the difference between the shims of the two acquisitions. Neither EPI image has signal in the voids above the sinus, where the field map predicts the largest displacement, so the estimate there is an extrapolation by the smoothness penalty. In practice you acquire the pair as two short series of a few volumes, estimate the displacement field once, and apply it to every volume of the functional or diffusion series acquired with one of the two polarities. As a spec --------- What this tour built, stated the way you would ask an agent for it: .. code-block:: text Correct the susceptibility distortion of a reversed phase-encoding EPI pair (BIDS, 3 T, phase encoding along y) with bartorch. Read both series with bartorch.io.read_nifti and average their volumes; compute the bandwidth per pixel along phase encoding as 1 / (AcquisitionMatrixPE * EffectiveEchoSpacing) from the sidecar. Run bt.correct_susceptibility(up, down, voxel_size=(z, y, x) in mm, phase_encoding_axis=1) and move the returned displacement (mm, on the faces between voxels) to the voxel centres. Report the correlation and RMS difference of the two images over the brain before and after, and compare the displacement with the one a dual-echo GRE field map predicts (correlation and slope). .. GENERATED FROM PYTHON SOURCE LINES 471-485 References ---------- .. [#andersson] Andersson JLR, Skare S, Ashburner J. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. *NeuroImage* 20(2):870-888 (2003). https://doi.org/10.1016/S1053-8119(03)00336-7 .. [#pyhysco] Julian A, Ruthotto L. PyHySCO: GPU-enabled susceptibility artifact distortion correction in seconds. *Front Neurosci* (2024). .. [#ds001600] Cieslak M, Elliott M, Satterthwaite T. Example Fieldmaps. OpenNeuro, accession ds001600. https://openneuro.org/datasets/ds001600. Licensed under CC-BY-SA. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.709 seconds) .. _sphx_glr_download_auto_examples_07-tours_07-epi-susceptibility-distortion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 07-epi-susceptibility-distortion.ipynb <07-epi-susceptibility-distortion.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 07-epi-susceptibility-distortion.py <07-epi-susceptibility-distortion.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 07-epi-susceptibility-distortion.zip <07-epi-susceptibility-distortion.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_