ENH: add gradunwarp base workflow for f/smriprep.

niworkflows/interfaces/gradunwarp.py

@@ -0,0 +1,83 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2021 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""GradUnwarp interface."""
+from pathlib import Path
+from nipype.interfaces.base import (
+    traits,
+    TraitedSpec,
+    File,
+    SimpleInterface
+)
+
+
+class _GradUnwarpInputSpec(TraitedSpec):
+    infile = File(exists=True, mandatory=True, desc="input image to be corrected")
+    gradfile = File(exists=True, default=None, desc="gradient file")
+    coeffile = File(exists=True, default=None, desc="coefficients file")
+    outfile = File(
+        "gradunwarped.nii.gz",
+        mandatory=True,
+        usedefault=True,
+        desc="output corrected image"
+    )
+    vendor = traits.Enum("siemens", "ge", usedefault=True, desc="scanner vendor")
+    warp = traits.Bool(desc="warp a volume (as opposed to unwarping)")
+    nojac = traits.Bool(desc="Do not perform Jacobian intensity correction")
+
+    fovmin = traits.Float(desc="the minimum extent of harmonics evaluation grid in meters")
+    fovmax = traits.Float(desc="the maximum extent of harmonics evaluation grid in meters")
+    order = traits.Int(
+        min=1, max=4,
+        usedefault=True,
+        desc="the order of interpolation(1..4) where 1 is linear - default")
+
+
+class _GradUnwarpOutputSpec(TraitedSpec):
+    corrected_file = File(desc="input images corrected")
+    warp_file = File(desc="absolute warp file")
+
+
+class GradUnwarp(SimpleInterface):
+    input_spec = _GradUnwarpInputSpec
+    output_spec = _GradUnwarpOutputSpec
+
+    def version():
+        try:
+            import gradunwarp
+        except ImportError:
+            return
+        return gradunwarp.__version__
+
+    def _run_interface(self, runtime):
+
+        if not GradUnwarp.version():
+            raise RuntimeError('missing gradunwarp dependency')
+        from gradunwarp.core.gradient_unwarp import GradientUnwarpRunner
+        gur = GradientUnwarpRunner(self.inputs)
+        gur.run()
+        gur.write()
+
+        cwd = Path(runtime.cwd)
+        self._results["corrected_file"] = str(cwd / self.inputs.outfile)
+        self._results["warp_file"] = str(cwd / "fullWarp_abs.nii.gz")
+        return runtime

niworkflows/testing.py

@@ -6,6 +6,8 @@
 from nipype.interfaces import afni, fsl
 from nipype.interfaces import freesurfer as fs
 
+from .interfaces import gradunwarp
+
 test_data_env = os.getenv('TEST_DATA_HOME', str(Path.home() / '.cache' / 'stanford-crn'))
 test_output_dir = os.getenv('TEST_OUTPUT_DIR')
 test_workdir = os.getenv('TEST_WORK_DIR')
@@ -44,3 +46,4 @@ def wrapper(*args, **kwargs):
 has_fsl = fsl.Info.version() is not None
 has_freesurfer = fs.Info.version() is not None
 has_afni = afni.Info.version() is not None
+has_gradunwarp = gradunwarp.GradUnwarp.version()

niworkflows/tests/data/gradunwarp_coeffs.grad

@@ -0,0 +1,65 @@
+#*[ Script ****************************************************************\
+#
+# Name        : dummy.grad
+#
+# Date        : 0970-01-01
+#
+# Author      : Mario Bros
+#
+# Language    :
+#
+# Description : Defines Legendre coefficients in spherical harmonics for
+#               dummy Gradient Coil 
+#
+#****************************************************************************/
+
+#*] END: */
+dummy, dummy , Gx,y,z = 80/80 mT/m
+win_low = 0, win_high = 0, win_algo = 0, win_dummy = 2;
+0.25 m = R0, lnorm = 4? A(1,0) = B(1,1) = A(1,1) = 0;
+0 = CoSyMode,
+NO.    TYPE                  SPECTRUM                                AXIS
+
+1 A( 3, 0)           -0.0510000                  z
+ 2 A( 5, 0)           -0.11000000                  z
+ 3 A( 7, 0)            0.04600000                  z
+ 4 A( 9, 0)           -0.00600000                  z
+ 5 A(11, 0)           -0.00400000                  z
+ 6 A(13, 0)            0.00400000                  z
+ 7 A(15, 0)           -0.00200000                  z
+ 8 A(17, 0)            0.00100000                  z
+ 9 A(19, 0)           -0.00000000                  z
+101 A( 3, 1)           -0.02500000                  x
+102 A( 3, 3)           -0.00300000                  x
+103 A( 5, 1)           -0.09400000                  x
+104 A( 5, 3)           -0.00300000                  x
+105 A( 5, 5)           -0.00200000                  x
+106 A( 7, 1)            0.01400000                  x
+107 A( 7, 3)            0.00700000                  x
+108 A( 7, 5)           -0.00100000                  x
+109 A( 7, 7)            0.00300000                  x
+110 A( 9, 1)            0.00800000                  x
+111 A( 9, 3)           -0.00400000                  x
+112 A(11, 1)           -0.00600000                  x
+113 A(11, 3)            0.00200000                  x
+114 A(13, 1)            0.00200000                  x
+115 A(13, 3)           -0.00100000                  x
+116 A(15, 1)           -0.00000000                  x
+117 A(15, 3)            0.00000000                  x
+201 B( 3, 1)           -0.03000000                  y
+202 B( 3, 3)            0.00700000                  y
+203 B( 5, 1)           -0.08800000                  y
+204 B( 5, 3)           -0.00400000                  y
+205 B( 5, 5)           -0.00100000                  y
+206 B( 7, 1)            0.01900000                  y
+207 B( 7, 3)           -0.00300000                  y
+208 B( 7, 5)           -0.00100000                  y
+209 B( 7, 7)            0.00100000                  y
+210 B( 9, 1)            0.00200000                  y
+211 B( 9, 3)            0.00300000                  y
+212 B( 9, 5)            0.00000000                  y
+213 B(11, 1)           -0.00300000                  y
+214 B(11, 3)           -0.00100000                  y
+215 B(13, 1)            0.00100000                  y
+216 B(13, 3)            0.00000000                  y
+217 B(15, 1)           -0.00000000                  y

niworkflows/workflows/gradunwarp.py

@@ -0,0 +1,101 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2021 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Workflow for the unwarping ."""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+
+from ..engine.workflows import LiterateWorkflow as Workflow
+
+
+def init_gradunwarp_wf(
+    name="gradunwarp_wf",
+):
+    from nipype.interfaces.fsl import ConvertWarp
+    from ..interfaces.gradunwarp import GradUnwarp
+
+    wf = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=["input_file", "coeff_file", "grad_file"]), name="inputnode"
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                "warp_file",
+                "corrected_file"
+            ]
+        ),
+        name="outputnode",
+    )
+
+    gradient_unwarp = pe.MapNode(
+        GradUnwarp(),
+        name="grad_unwarp",
+        iterfield=["infile"]
+    )
+
+    convert_warp = pe.MapNode(
+        ConvertWarp(
+            abswarp=True,
+            out_relwarp=True,
+        ),
+        name='convert_warp_abs2rel',
+        iterfield=["reference", "warp1"]
+    )
+
+    def _warp_fsl2itk(in_warp):
+        import os
+        import nitransforms.io
+        from nipype.utils.filemanip import fname_presuffix
+        fsl_warp = nitransforms.io.fsl.FSLDisplacementsField.from_filename(in_warp)
+        out_fname = fname_presuffix(in_warp, suffix="_itk", newpath=os.getcwd())
+        nitransforms.io.itk.ITKDisplacementsField.to_filename(fsl_warp, out_fname)
+        return out_fname
+
+    warp_fsl2itk = pe.MapNode(
+        niu.Function(
+            function=_warp_fsl2itk,
+            output_names=["itk_warp"],
+            input_names=["in_warp"]
+        ),
+        iterfield=['in_warp'],
+        name="warp_fsl2itk"
+    )
+
+    # fmt:off
+    wf.connect([
+        (inputnode, gradient_unwarp, [
+            ("input_file", "infile"),
+            ("coeff_file", "coeffile"),
+            ("grad_file", "gradfile"),
+        ]),
+        (inputnode, convert_warp, [("input_file", "reference")]),
+        (gradient_unwarp, convert_warp, [("warp_file", "warp1")]),
+        (gradient_unwarp, outputnode, [("corrected_file", "corrected_file")]),
+        (convert_warp, warp_fsl2itk, [("out_file", "in_warp")]),
+        (warp_fsl2itk, outputnode, [("itk_warp", "warp_file")])
+    ])
+    # fmt:on
+
+    return wf

niworkflows/workflows/tests/test_gradunwarp.py

@@ -0,0 +1,43 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2021 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Check the refmap module."""
+import unittest
+from ..gradunwarp import init_gradunwarp_wf
+from ...testing import has_gradunwarp
+from niworkflows.tests.data import load_test_data
+
+
+@unittest.skipUnless(has_gradunwarp, "Needs Gradunwarp")
+def test_gradunwarp(tmpdir, ds000030_dir, workdir, outdir):
+    """Exercise the Gradunwarp workflow."""
+    tmpdir.chdir()
+
+    wf = init_gradunwarp_wf()
+    if workdir:
+        wf.base_dir = str(workdir)
+
+    print(ds000030_dir)
+    wf.inputs.inputnode.input_file = str(next(ds000030_dir.glob("sub-10228/anat/*_T1w.nii.gz")))
+    wf.inputs.inputnode.grad_file = str(load_test_data('gradunwarp_coeffs.grad'))
+
+    wf.run()

pyproject.toml

@@ -26,13 +26,14 @@ dependencies = [
     "acres",
     "attrs >=20.1",
     "importlib_resources >= 5.7; python_version < '3.11'",
+    "gradunwarp@https://github.com/nipreps/gradunwarp/archive/refs/heads/master.zip",
     "jinja2 >=3",
     "looseversion",
     "matplotlib >= 3.5",
     "nibabel >= 3.0",
     "nilearn >= 0.8",
     "nipype >= 1.8.5",
-    "nitransforms >= 22.0.0",
+    "nitransforms >= 24.0.0",
     "numpy >= 1.20",
     "packaging",
     "pandas >= 1.2",