PyHySCO: GPU-Enabled Susceptibility Artifact Distortion Correction in Seconds

TL;DR

PyTorch Hyperelastic Susceptibility Correction (PyHySCO), a user-friendly EPI distortion correction tool implemented in PyTorch that enables multi-threading and efficient use of graphics processing units (GPUs) that is reliable without training is introduced.

Abstract

Over the past decade, reversed Gradient Polarity (RGP) methods have become a popular approach for correcting susceptibility artifacts in Echo-Planar Imaging (EPI). Although several post-processing tools for RGP are available, their implementations do not fully leverage recent hardware, algorithmic, and computational advances, leading to correction times of several minutes per image volume. To enable 3D RGP correction in seconds, we introduce PyHySCO, a user-friendly EPI distortion correction tool implemented in PyTorch that enables multi-threading and efficient use of graphics processing units (GPUs). PyHySCO uses a time-tested physical distortion model and mathematical formulation and is, therefore, reliable without training. An algorithmic improvement in PyHySCO is its novel initialization scheme that uses 1D optimal transport. PyHySCO is published under the GNU public license and can be used from the command line or its Python interface. Our extensive numerical validation using 3T and 7T data from the Human Connectome Project suggests that PyHySCO achieves accuracy comparable to that of leading RGP tools at a fraction of the cost. We also validate the new initialization scheme, compare different optimization algorithms, and test the algorithm on different hardware and arithmetic precision.

Figures (9)

• Figure 1: The Reverse Gradient Polarity correction paradigm. Two images are acquired with opposite phase encoding directions, $+v$ and $-v$. These two images are used to estimate the field map $b$, and the distortion correction model ChangFitzpatrick1992 is applied to obtain a corrected image $\mathcal{I}$.
• Figure 2: Example of one-dimensional optimal transport maps. The top left shows an example of one-dimensional measures. The green signal, $i_{+v}$, corresponds to an intensity pileup in $\mathcal{I}_{+v}$, while the purple signal $i_{-v}$ corresponds to an intensity dispersion in $\mathcal{I}_{-v}$. The red signal corresponds to the intensity of the true image. The top right shows the cumulative distributions for the measures $i_{+v}$ and $i_{-v}$. Bottom left shows the pseudoinverses for $i_{+v}$ and $i_{-v}$ along with the pseudoinverse $C_{i_{\text{half}}}^{-1}$ used in calculating the transport maps $T_{+}=C^{-1}_{i_{\text{half}}} \circ C_{i_{+v}}$ and $T_{-}=C^{-1}_{i_{\text{half}}} \circ C_{i_{-v}}$, shown bottom right.
• Figure 3: The maps $T_{+}$ and $T_{-}$ are calculated using the closed-form one-dimensional optimal transport solution, parallelized in the distortion dimension peyreOT. Note the inverted coloring between $T_{+}$ and $T_{-}$ as the map $T_{-}$ corrects a distortion in the opposite direction as $T_{+}$.
• Figure 4: UML diagram of PyHySCO showing the classes and relationships for the (\ref{['fig:UML_loss']}) loss function and (\ref{['fig:UML_opt']}) optimization. A EPIMRIDistortionCorrection object defining the loss function is an attribute of every EPIOptimize object defining the optimization scheme.
• Figure 5: The usage and workflow of PyHySCO.