enh: implement a unit-test for SDC-SyN

oesteban · oesteban · commit 72219d36653a · 2020-12-19T02:19:35.000-08:00
Resolves: #163.
diff --git a/.circleci/config.yml b/.circleci/config.yml
@@ -109,7 +109,7 @@ jobs:
           command: |
             datalad install -r https://github.com/nipreps-data/ds000054.git
             datalad update --merge -d ds000054/
-            datalad get -r -d ds000054/ ds000054/sub-100185/fmap/*
+            datalad get -r -J 2 -d ds000054/ ds000054/*
 
       - save_cache:
           key: data-v4-{{ .Branch }}-{{ .BuildNum }}
diff --git a/sdcflows/interfaces/bspline.py b/sdcflows/interfaces/bspline.py
@@ -240,10 +240,10 @@ def _run_interface(self, runtime):
 
         # Calculate the physical coordinates of target grid
         targetnii = nb.load(self.inputs.in_target)
-        targetaff = targetnii.affine
         allmask = np.ones_like(targetnii.dataobj, dtype="uint8")
-        voxels = np.argwhere(allmask == 1).astype("float32")
-        points = apply_affine(targetaff.astype("float32"), voxels)
+        # targetaff = targetnii.affine
+        # voxels = np.argwhere(allmask == 1).astype("float32")
+        # points = apply_affine(targetaff.astype("float32"), voxels)
 
         weights = []
         coeffs = []
diff --git a/sdcflows/workflows/fit/syn.py b/sdcflows/workflows/fit/syn.py
@@ -3,6 +3,15 @@
 """
 Estimating the susceptibility distortions without fieldmaps.
 
+.. testsetup::
+
+    >>> tmpdir = getfixture('tmpdir')
+    >>> tmp = tmpdir.chdir() # changing to a temporary directory
+    >>> data = np.zeros((10, 10, 10, 1, 3))
+    >>> data[..., 1] = 1
+    >>> nb.Nifti1Image(data, None, None).to_filename(
+    ...     tmpdir.join('field.nii.gz').strpath)
+
 .. _sdc_fieldmapless :
 
 Fieldmap-less approaches
@@ -55,7 +64,11 @@
 
 
 def init_syn_sdc_wf(
-    *, atlas_threshold=3, debug=False, name="syn_sdc_wf", omp_nthreads=1,
+    *,
+    atlas_threshold=3,
+    debug=False,
+    name="syn_sdc_wf",
+    omp_nthreads=1,
 ):
     """
     Build the *fieldmap-less* susceptibility-distortion estimation workflow.
@@ -99,7 +112,7 @@ def init_syn_sdc_wf(
         A preprocessed, skull-stripped anatomical (T1w or T2w) image.
     std2anat_xfm : :obj:`str`
         inverse registration transform of T1w image to MNI template
-    anat2bold_xfm : :obj:`str`
+    anat2epi_xfm : :obj:`str`
         transform mapping coordinates from the EPI space to the anatomical
         space (i.e., the transform to resample anatomical info into EPI space.)
 
@@ -120,7 +133,14 @@ def init_syn_sdc_wf(
         FixHeaderRegistration as Registration,
     )
     from niworkflows.interfaces.nibabel import Binarize
+    from ...interfaces.epi import EPIMask
     from ...utils.misc import front as _pop
+    from ...interfaces.bspline import (
+        BSplineApprox,
+        DEFAULT_LF_ZOOMS_MM,
+        DEFAULT_HF_ZOOMS_MM,
+        DEFAULT_ZOOMS_MM,
+    )
 
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""\
@@ -137,12 +157,13 @@ def init_syn_sdc_wf(
 """
     inputnode = pe.Node(
         niu.IdentityInterface(
-            ["epi_ref", "epi_mask", "anat_brain", "std2anat_xfm", "anat2bold_xfm"]
+            ["epi_ref", "epi_mask", "anat_brain", "std2anat_xfm", "anat2epi_xfm"]
         ),
         name="inputnode",
     )
     outputnode = pe.Node(
-        niu.IdentityInterface(["fmap", "fmap_ref", "fmap_mask"]), name="outputnode",
+        niu.IdentityInterface(["fmap", "fmap_ref", "fmap_coeff", "fmap_mask"]),
+        name="outputnode",
     )
 
     invert_t1w = pe.Node(Rescale(invert=True), name="invert_t1w", mem_gb=0.3)
@@ -174,29 +195,54 @@ def init_syn_sdc_wf(
         n_procs=omp_nthreads,
     )
 
-    unwarp_ref = pe.Node(ApplyTransforms(interpolation="BSpline"), name="unwarp_ref",)
+    unwarp_ref = pe.Node(
+        ApplyTransforms(interpolation="BSpline"),
+        name="unwarp_ref",
+    )
+
+    epi_mask = pe.Node(EPIMask(), name="epi_mask")
+
+    # Extract nonzero component
+    extract_field = pe.Node(niu.Function(function=_extract_field), name="extract_field")
+
+    # Regularize with B-Splines
+    bs_filter = pe.Node(BSplineApprox(), n_procs=omp_nthreads, name="bs_filter")
+    bs_filter.interface._always_run = debug
+    bs_filter.inputs.bs_spacing = (
+        [DEFAULT_LF_ZOOMS_MM, DEFAULT_HF_ZOOMS_MM] if not debug else [DEFAULT_ZOOMS_MM]
+    )
+    bs_filter.inputs.extrapolate = not debug
 
     # fmt: off
     workflow.connect([
-        (inputnode, transform_list, [("anat2bold_xfm", "in1"),
+        (inputnode, transform_list, [("anat2epi_xfm", "in1"),
                                      ("std2anat_xfm", "in2")]),
         (inputnode, invert_t1w, [("anat_brain", "in_file"),
                                  (("epi_ref", _pop), "ref_file")]),
-        (inputnode, anat2epi, [(("epi_ref", _pop), "reference_image")]),
+        (inputnode, anat2epi, [(("epi_ref", _pop), "reference_image"),
+                               ("anat2epi_xfm", "transforms")]),
         (inputnode, syn, [(("epi_ref", _pop), "moving_image"),
                           ("epi_mask", "moving_image_masks"),
                           (("epi_ref", _warp_dir), "restrict_deformation")]),
+        (inputnode, unwarp_ref, [(("epi_ref", _pop), "reference_image"),
+                                 (("epi_ref", _pop), "input_image")]),
         (inputnode, prior2epi, [(("epi_ref", _pop), "reference_image")]),
+        (inputnode, extract_field, [("epi_ref", "epi_meta")]),
         (invert_t1w, anat2epi, [("out_file", "input_image")]),
         (transform_list, prior2epi, [("out", "transforms")]),
         (prior2epi, atlas_msk, [("output_image", "in_file")]),
         (anat2epi, syn, [("output_image", "fixed_image")]),
         (atlas_msk, syn, [(("out_mask", _fixed_masks_arg), "fixed_image_masks")]),
-        (syn, outputnode, [("forward_transforms", "fmap")]),
+        (syn, extract_field, [("forward_transforms", "in_file")]),
         (syn, unwarp_ref, [("forward_transforms", "transforms")]),
-        (inputnode, unwarp_ref, [(("epi_ref", _pop), "reference_image"),
-                                 (("epi_ref", _pop), "input_image")]),
+        (unwarp_ref, epi_mask, [("output_image", "in_file")]),
+        (extract_field, bs_filter, [("out", "in_data")]),
+        (epi_mask, bs_filter, [("out_file", "in_mask")]),
         (unwarp_ref, outputnode, [("output_image", "fmap_ref")]),
+        (epi_mask, outputnode, [("out_file", "fmap_mask")]),
+        (bs_filter, outputnode, [
+            ("out_extrapolated" if not debug else "out_field", "fmap"),
+            ("out_coeff", "fmap_coeff")]),
     ])
     # fmt: on
 
@@ -231,3 +277,41 @@ def _fixed_masks_arg(mask):
 
     """
     return ["NULL", mask]
+
+
+def _extract_field(in_file, epi_meta):
+    """
+    Extract the nonzero component of the deformation field estimated by ANTs.
+
+    Examples
+    --------
+    >>> nii = nb.load(
+    ...     _extract_field(
+    ...         ["field.nii.gz"],
+    ...         ("epi.nii.gz", {"PhaseEncodingDirection": "j-", "TotalReadoutTime": 0.005}))
+    ... )
+    >>> nii.shape
+    (10, 10, 10)
+
+    >>> np.allclose(nii.get_fdata(), 200)
+    True
+
+    """
+    from nipype.utils.filemanip import fname_presuffix
+    import numpy as np
+    import nibabel as nb
+    from sdcflows.utils.epimanip import get_trt
+
+    fieldnii = nb.load(in_file[0])
+    trt = get_trt(epi_meta[1], in_file=epi_meta[0])
+    data = (
+        np.squeeze(fieldnii.get_fdata(dtype="float32"))[
+            ..., "ijk".index(epi_meta[1]["PhaseEncodingDirection"][0])
+        ]
+        / trt
+    )
+    out_file = fname_presuffix(in_file[0], suffix="_fieldmap")
+    nii = nb.Nifti1Image(data, fieldnii.affine, None)
+    nii.header.set_xyzt_units(fieldnii.header.get_xyzt_units()[0])
+    nii.to_filename(out_file)
+    return out_file
diff --git a/sdcflows/workflows/fit/tests/test_syn.py b/sdcflows/workflows/fit/tests/test_syn.py
@@ -0,0 +1,104 @@
+"""Test fieldmap-less SDC-SyN."""
+import os
+import pytest
+from nipype.pipeline import engine as pe
+from niworkflows.interfaces.nibabel import ApplyMask
+
+from ..syn import init_syn_sdc_wf
+
+
+@pytest.mark.skipif(os.getenv("TRAVIS") == "true", reason="this is TravisCI")
+@pytest.mark.skipif(os.getenv("GITHUB_ACTIONS") == "true", reason="this is GH Actions")
+def test_syn_wf(tmpdir, datadir, workdir, outdir):
+    """Build and run an SDC-SyN workflow."""
+    derivs_path = datadir / "ds000054" / "derivatives"
+    smriprep = derivs_path / "smriprep-0.6" / "sub-100185" / "anat"
+
+    wf = pe.Workflow(name="syn_test")
+
+    syn_wf = init_syn_sdc_wf(debug=True, omp_nthreads=4)
+    syn_wf.inputs.inputnode.epi_ref = (
+        str(
+            derivs_path
+            / "sdcflows-tests"
+            / "sub-100185_task-machinegame_run-1_boldref.nii.gz"
+        ),
+        {"PhaseEncodingDirection": "j", "TotalReadoutTime": 0.005},
+    )
+    syn_wf.inputs.inputnode.epi_mask = str(
+        derivs_path
+        / "sdcflows-tests"
+        / "sub-100185_task-machinegame_run-1_desc-brain_mask.nii.gz"
+    )
+    syn_wf.inputs.inputnode.anat2epi_xfm = str(
+        derivs_path / "sdcflows-tests" / "t1w2bold.txt"
+    )
+    syn_wf.inputs.inputnode.std2anat_xfm = str(
+        smriprep / "sub-100185_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5"
+    )
+
+    t1w_mask = pe.Node(
+        ApplyMask(
+            in_file=str(smriprep / "sub-100185_desc-preproc_T1w.nii.gz"),
+            in_mask=str(smriprep / "sub-100185_desc-brain_mask.nii.gz"),
+        ),
+        name="t1w_mask",
+    )
+
+    wf.connect(
+        [
+            (t1w_mask, syn_wf, [("out_file", "inputnode.anat_brain")]),
+        ]
+    )
+
+    if outdir:
+        from ...outputs import init_fmap_derivatives_wf, init_fmap_reports_wf
+
+        outdir = outdir / "unittests" / "syn_test"
+        fmap_derivatives_wf = init_fmap_derivatives_wf(
+            output_dir=str(outdir),
+            write_coeff=True,
+            custom_entities={"est": "sdcsyn"},
+            bids_fmap_id="sdcsyn",
+        )
+        fmap_derivatives_wf.inputs.inputnode.source_files = [
+            str(
+                derivs_path
+                / "sdcflows-tests"
+                / "sub-100185_task-machinegame_run-1_boldref.nii.gz"
+            )
+        ]
+        fmap_derivatives_wf.inputs.inputnode.fmap_meta = {
+            "PhaseEncodingDirection": "j-"
+        }
+
+        fmap_reports_wf = init_fmap_reports_wf(
+            output_dir=str(outdir),
+            fmap_type="pepolar",
+        )
+        fmap_reports_wf.inputs.inputnode.source_files = [
+            str(
+                derivs_path
+                / "sdcflows-tests"
+                / "sub-100185_task-machinegame_run-1_boldref.nii.gz"
+            )
+        ]
+
+        # fmt: off
+        wf.connect([
+            (syn_wf, fmap_reports_wf, [
+                ("outputnode.fmap", "inputnode.fieldmap"),
+                ("outputnode.fmap_ref", "inputnode.fmap_ref"),
+                ("outputnode.fmap_mask", "inputnode.fmap_mask")]),
+            (syn_wf, fmap_derivatives_wf, [
+                ("outputnode.fmap", "inputnode.fieldmap"),
+                ("outputnode.fmap_ref", "inputnode.fmap_ref"),
+                ("outputnode.fmap_coeff", "inputnode.fmap_coeff"),
+            ]),
+        ])
+        # fmt: on
+
+    if workdir:
+        wf.base_dir = str(workdir)
+
+    wf.run(plugin="Linear")
