RF: Replace Convert3d with nitransforms in MCFLIRT2ITK (#835)

* TEST: Add regression test for MCFLIRT2ITK

* TEST: Update test to account for minor differences with nitransforms

* RF: Replace Convert3d with nitransforms in MCFLIRT2ITK

* MNT: Remove convert3d from the Dockerfile

Author: effigies

Dockerfile

@@ -152,17 +152,6 @@ ENV FSLDIR="/opt/fsl" \
 ENV MKL_NUM_THREADS=1 \
     OMP_NUM_THREADS=1
 
-# Convert3D (neurodocker build)
-RUN echo "Downloading Convert3D ..." \
-    && mkdir -p /opt/convert3d-1.0.0 \
-    && curl -fsSL --retry 5 https://sourceforge.net/projects/c3d/files/c3d/1.0.0/c3d-1.0.0-Linux-x86_64.tar.gz/download \
-    | tar -xz -C /opt/convert3d-1.0.0 --strip-components 1 \
-    --exclude "c3d-1.0.0-Linux-x86_64/lib" \
-    --exclude "c3d-1.0.0-Linux-x86_64/share" \
-    --exclude "c3d-1.0.0-Linux-x86_64/bin/c3d_gui"
-ENV C3DPATH="/opt/convert3d-1.0.0" \
-    PATH="/opt/convert3d-1.0.0/bin:$PATH"
-
 # Create a shared $HOME directory
 RUN useradd -m -s /bin/bash -G users niworkflows
 WORKDIR /home/niworkflows

niworkflows/interfaces/itk.py

@@ -25,17 +25,22 @@
 from mimetypes import guess_type
 from tempfile import TemporaryDirectory
 
+import nibabel as nb
+import nitransforms as nt
+import numpy as np
 from nipype import logging
-from nipype.utils.filemanip import fname_presuffix
 from nipype.interfaces.base import (
-    traits,
-    TraitedSpec,
     BaseInterfaceInputSpec,
     File,
     InputMultiObject,
     OutputMultiObject,
     SimpleInterface,
+    TraitedSpec,
+    isdefined,
+    traits,
 )
+from nipype.utils.filemanip import fname_presuffix
+
 from .fixes import _FixTraitApplyTransformsInputSpec
 
 LOGGER = logging.getLogger("nipype.interface")
@@ -45,13 +50,9 @@ class _MCFLIRT2ITKInputSpec(BaseInterfaceInputSpec):
     in_files = InputMultiObject(
         File(exists=True), mandatory=True, desc="list of MAT files from MCFLIRT"
     )
-    in_reference = File(
-        exists=True, mandatory=True, desc="input image for spatial reference"
-    )
+    in_reference = File(exists=True, mandatory=True, desc="input image for spatial reference")
     in_source = File(exists=True, mandatory=True, desc="input image for spatial source")
-    num_threads = traits.Int(
-        1, usedefault=True, nohash=True, desc="number of parallel processes"
-    )
+    num_threads = traits.Int(nohash=True, desc="number of parallel processes")
 
 
 class _MCFLIRT2ITKOutputSpec(TraitedSpec):
@@ -65,53 +66,22 @@ class MCFLIRT2ITK(SimpleInterface):
     output_spec = _MCFLIRT2ITKOutputSpec
 
     def _run_interface(self, runtime):
-        num_threads = self.inputs.num_threads
-        if num_threads < 1:
-            num_threads = None
-
-        with TemporaryDirectory(prefix="tmp-", dir=runtime.cwd) as tmp_folder:
-            # Inputs are ready to run in parallel
-            if num_threads is None or num_threads > 1:
-                from concurrent.futures import ThreadPoolExecutor
-
-                with ThreadPoolExecutor(max_workers=num_threads) as pool:
-                    itk_outs = list(
-                        pool.map(
-                            _mat2itk,
-                            [
-                                (
-                                    in_mat,
-                                    self.inputs.in_reference,
-                                    self.inputs.in_source,
-                                    i,
-                                    tmp_folder,
-                                )
-                                for i, in_mat in enumerate(self.inputs.in_files)
-                            ],
-                        )
-                    )
-            else:
-                itk_outs = [
-                    _mat2itk(
-                        (
-                            in_mat,
-                            self.inputs.in_reference,
-                            self.inputs.in_source,
-                            i,
-                            tmp_folder,
-                        )
-                    )
-                    for i, in_mat in enumerate(self.inputs.in_files)
-                ]
-
-        # Compose the collated ITK transform file and write
-        tfms = "#Insight Transform File V1.0\n" + "".join(
-            [el[1] for el in sorted(itk_outs)]
+        if isdefined(self.inputs.num_threads):
+            LOGGER.warning("Multithreading is deprecated. Remove the num_threads input.")
+
+        source = nb.load(self.inputs.in_source)
+        reference = nb.load(self.inputs.in_reference)
+        affines = [
+            nt.linear.load(mat, fmt='fsl', reference=reference, moving=source)
+            for mat in self.inputs.in_files
+        ]
+
+        affarray = nt.io.itk.ITKLinearTransformArray.from_ras(
+            np.stack([a.matrix for a in affines], axis=0),
         )
 
         self._results["out_file"] = os.path.join(runtime.cwd, "mat2itk.txt")
-        with open(self._results["out_file"], "w") as f:
-            f.write(tfms)
+        affarray.to_filename(self._results["out_file"])
 
         return runtime
 
@@ -206,30 +176,6 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def _mat2itk(args):
-    from nipype.interfaces.c3 import C3dAffineTool
-    from nipype.utils.filemanip import fname_presuffix
-
-    in_file, in_ref, in_src, index, newpath = args
-    # Generate a temporal file name
-    out_file = fname_presuffix(in_file, suffix="_itk-%05d.txt" % index, newpath=newpath)
-
-    # Run c3d_affine_tool
-    C3dAffineTool(
-        transform_file=in_file,
-        reference_file=in_ref,
-        source_file=in_src,
-        fsl2ras=True,
-        itk_transform=out_file,
-        resource_monitor=False,
-    ).run()
-    transform = "#Transform %d\n" % index
-    with open(out_file) as itkfh:
-        transform += "".join(itkfh.readlines()[2:])
-
-    return (index, transform)
-
-
 def _applytfms(args):
     """
     Applies ANTs' antsApplyTransforms to the input image.

niworkflows/interfaces/tests/data/MAT_0098

@@ -0,0 +1,4 @@
+0.999999 -0.000272 0.001561 -0.071358 
+0.000272 1.000000 -0.000210 0.053746 
+-0.001561 0.000210 0.999999 0.153994 
+0.000000 0.000000 0.000000 1.000000 

niworkflows/interfaces/tests/data/MAT_0099

@@ -0,0 +1,4 @@
+0.999999 -0.000320 0.001330 -0.045489 
+0.000320 1.000000 -0.000390 0.064510 
+-0.001329 0.000390 0.999999 0.111615 
+0.000000 0.000000 0.000000 1.000000 

niworkflows/interfaces/tests/data/boldref.nii

@@ -20,19 +20,26 @@
 #
 #     https://www.nipreps.org/community/licensing/
 #
+import re
+from pathlib import Path
+
+import numpy as np
 import pytest
-from ..itk import _applytfms
-from ... import data
 from nipype.interfaces.ants.base import Info
+from nipype.pipeline import engine as pe
+
+from ... import data
+from ..itk import MCFLIRT2ITK, _applytfms
+from .data import load_test_data
 
 
 @pytest.mark.skipif(Info.version() is None, reason="Missing ANTs")
 @pytest.mark.parametrize("ext", (".nii", ".nii.gz"))
 @pytest.mark.parametrize("copy_dtype", (True, False))
 @pytest.mark.parametrize("in_dtype", ("i2", "f4"))
 def test_applytfms(tmpdir, ext, copy_dtype, in_dtype):
-    import numpy as np
     import nibabel as nb
+    import numpy as np
 
     in_file = str(tmpdir / ("src" + ext))
     nii = nb.Nifti1Image(np.zeros((5, 5, 5), dtype=np.float32), np.eye(4))
@@ -52,3 +59,81 @@ def test_applytfms(tmpdir, ext, copy_dtype, in_dtype):
     assert np.allclose(nii.get_fdata(), out_nii.get_fdata())
     if copy_dtype:
         assert nii.get_data_dtype() == out_nii.get_data_dtype()
+
+
+def test_MCFLIRT2ITK(tmp_path):
+    # Test that MCFLIRT2ITK produces output that is consistent with convert3d
+    test_data = load_test_data()
+
+    fsl2itk = pe.Node(
+        MCFLIRT2ITK(
+            in_files=[str(test_data / "MAT_0098"), str(test_data / "MAT_0099")],
+            in_reference=str(test_data / "boldref.nii"),
+            in_source=str(test_data / "boldref.nii"),
+        ),
+        name="fsl2itk",
+        base_dir=str(tmp_path),
+    )
+
+    res = fsl2itk.run()
+    out_file = Path(res.outputs.out_file)
+
+    assert out_file.exists()
+    lines = out_file.read_text().splitlines()
+
+    assert lines[:2] == [
+        "#Insight Transform File V1.0",
+        "#Transform 0",
+    ]
+    assert re.match(
+        r"Transform: (MatrixOffsetTransformBase|AffineTransform)_(float|double)_3_3",
+        lines[2],
+    )
+    assert lines[3].startswith("Parameters: ")
+    assert lines[4] == "FixedParameters: 0 0 0"
+    offset = 1 if lines[5] == "" else 0
+    assert lines[5 + offset] == "#Transform 1"
+    assert lines[6 + offset] == lines[2]
+    assert lines[7 + offset].startswith("Parameters: ")
+
+    params0 = np.array([float(p) for p in lines[3].split(" ")[1:]])
+    params1 = np.array([float(p) for p in lines[7 + offset].split(" ")[1:]])
+    # Empirically determined
+    assert np.allclose(
+        params0,
+        np.array(
+            [
+                9.99998489e-01,
+                -4.36657508e-04,
+                -1.52316526e-03,
+                4.36017740e-04,
+                9.99999777e-01,
+                -2.10558666e-04,
+                1.52334852e-03,
+                2.09440681e-04,
+                9.99998624e-01,
+                -1.28961869e-03,
+                6.93155516e-02,
+                -1.12375673e-02,
+            ]
+        ),
+    )
+    assert np.allclose(
+        params1,
+        np.array(
+            [
+                9.99999130e-01,
+                -4.60021530e-04,
+                -1.28828576e-03,
+                4.58910652e-04,
+                9.99999648e-01,
+                -3.90485877e-04,
+                1.28764980e-03,
+                3.89513646e-04,
+                9.99999178e-01,
+                -9.19541650e-03,
+                7.45419094e-02,
+                -8.95843238e-03,
+            ]
+        ),
+    )