Repository: nipreps/nibabies

Issue:
Piecewise subcortical alignment for CIFTI mapping
To have a persistent link, here's the current script from DCAN: [SubcorticalAlign_ROIs.sh](https://github.com/DCAN-Labs/dcan-infant-pipeline/blob/247e19e5441cc814cea2f23720caeeb6c6aeadf8/fMRISurface/scripts/SubcorticalAlign_ROIs.sh)

Note that this supersedes an older script: [SubcorticalAlign_MNI_Affine.sh](https://github.com/DCAN-Labs/dcan-infant-pipeline/blob/247e19e5441cc814cea2f23720caeeb6c6aeadf8/fMRISurface/scripts/SubcorticalAlign_MNI_Affine.sh)

Discussion:
pep8speaks: Hello @mgxd! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`nibabies/workflows/bold/alignment.py`](https://github.com/nipreps/nibabies/blob/18b38200536df7aab9c2ffe307d6ac9372ce035b/nibabies/workflows/bold/alignment.py):

> [Line 67:5](https://github.com/nipreps/nibabies/blob/18b38200536df7aab9c2ffe307d6ac9372ce035b/nibabies/workflows/bold/alignment.py#L67): [E303](https://duckduckgo.com/?q=pep8%20E303) too many blank lines (2)
> [Line 80:5](https://github.com/nipreps/nibabies/blob/18b38200536df7aab9c2ffe307d6ac9372ce035b/nibabies/workflows/bold/alignment.py#L80): [E266](https://duckduckgo.com/?q=pep8%20E266) too many leading '#' for block comment
> [Line 81:5](https://github.com/nipreps/nibabies/blob/18b38200536df7aab9c2ffe307d6ac9372ce035b/nibabies/workflows/bold/alignment.py#L81): [E266](https://duckduckgo.com/?q=pep8%20E266) too many leading '#' for block comment


codecov-commenter: # [Codecov](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) Report
> Merging [#72](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) (8748be7) into [master](https://codecov.io/gh/nipreps/nibabies/commit/06b6a4b6c730a65e4774ec2bb841fec86355839f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) (06b6a4b) will **increase** coverage by `0.38%`.
> The diff coverage is `50.35%`.

[![Impacted file tree graph](https://codecov.io/gh/nipreps/nibabies/pull/72/graphs/tree.svg?width=650&height=150&src=pr&token=LFZGJHJSDF&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps)](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps)

```diff
@@            Coverage Diff             @@
##           master      #72      +/-   ##
==========================================
+ Coverage   31.75%   32.13%   +0.38%     
==========================================
  Files          42       43       +1     
  Lines        3228     3348     +120     
==========================================
+ Hits         1025     1076      +51     
- Misses       2203     2272      +69     
```


| [Impacted Files](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) | Coverage Δ | |
|---|---|---|
| [nibabies/workflows/bold/alignment.py](https://codecov.io/gh/nipreps/nibabies/pull/72/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps#diff-bmliYWJpZXMvd29ya2Zsb3dzL2JvbGQvYWxpZ25tZW50LnB5) | `15.38% <15.38%> (ø)` | |
| [nibabies/interfaces/workbench.py](https://codecov.io/gh/nipreps/nibabies/pull/72/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps#diff-bmliYWJpZXMvaW50ZXJmYWNlcy93b3JrYmVuY2gucHk=) | `90.38% <71.26%> (-9.62%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps). Last update [06b6a4b...8748be7](https://codecov.io/gh/nipreps/nibabies/pull/72?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps).


Patch:
diff --git a/nibabies/config.py b/nibabies/config.py
index 2d8f4a1..ef56498 100644
--- a/nibabies/config.py
+++ b/nibabies/config.py
@@ -684,6 +684,12 @@ def init_spaces(checkpoint=True):
     if checkpoint and not spaces.is_cached():
         spaces.checkpoint()
 
+    if workflow.age_months is not None:
+        from .utils.misc import cohort_by_months
+
+        if "MNIInfant" not in spaces.get_spaces(nonstandard=False, dim=(3,)):
+            cohort = cohort_by_months("MNIInfant", workflow.age_months)
+            spaces.add(Reference("MNIInfant", {"cohort": cohort}))
     # # Add the default standard space if not already present (required by several sub-workflows)
     # if "MNI152NLin2009cAsym" not in spaces.get_spaces(nonstandard=False, dim=(3,)):
     #     spaces.add(Reference("MNI152NLin2009cAsym", {}))
@@ -695,10 +701,9 @@ def init_spaces(checkpoint=True):
     #     # Make sure there's a normalization to FSL for AROMA to use.
     #     spaces.add(Reference("MNI152NLin6Asym", {"res": "2"}))
 
-    cifti_output = workflow.cifti_output
     if workflow.cifti_output:
         # CIFTI grayordinates to corresponding FSL-MNI resolutions.
-        vol_res = "2" if cifti_output == "91k" else "1"
+        vol_res = "2" if workflow.cifti_output == "91k" else "1"
         spaces.add(Reference("fsaverage", {"den": "164k"}))
         spaces.add(Reference("MNI152NLin6Asym", {"res": vol_res}))
 
diff --git a/nibabies/interfaces/nibabel.py b/nibabies/interfaces/nibabel.py
index 7bf1e14..fc56f8d 100644
--- a/nibabies/interfaces/nibabel.py
+++ b/nibabies/interfaces/nibabel.py
@@ -5,6 +5,7 @@ from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
     File,
     SimpleInterface,
+    InputMultiObject,
 )
 
 
@@ -34,6 +35,24 @@ class BinaryDilation(SimpleInterface):
         return runtime
 
 
+class MergeROIsInputSpec(BaseInterfaceInputSpec):
+    in_files = InputMultiObject(File(exists=True), desc="ROI files to be merged")
+
+
+class MergeROIsOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="NIfTI containing all ROIs")
+
+
+class MergeROIs(SimpleInterface):
+    """Combine multiple region of interest files (3D or 4D) into a single file"""
+    input_spec = MergeROIsInputSpec
+    output_spec = MergeROIsOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results["out_file"] = _merge_rois(self.inputs.in_files, newpath=runtime.cwd)
+        return runtime
+
+
 def _dilate(in_file, radius=3, iterations=1, newpath=None):
     """Dilate (binary) input mask."""
     from pathlib import Path
@@ -55,3 +74,41 @@ def _dilate(in_file, radius=3, iterations=1, newpath=None):
     out_file = fname_presuffix(in_file, suffix="_dil", newpath=newpath or Path.cwd())
     mask.__class__(newdata.astype("uint8"), mask.affine, hdr).to_filename(out_file)
     return out_file
+
+
+def _merge_rois(in_files, newpath=None):
+    """
+    Aggregate individual 4D ROI files together into a single subcortical NIfTI.
+    All ROI images are sanity checked with regards to:
+    1) Shape
+    2) Affine
+    3) Overlap
+
+    If any of these checks fail, an ``AssertionError`` will be raised.
+    """
+    from pathlib import Path
+    import nibabel as nb
+    import numpy as np
+
+    img = nb.load(in_files[0])
+    data = np.array(img.dataobj)
+    affine = img.affine
+    header = img.header
+
+    nonzero = np.any(data, axis=3)
+    for roi in in_files[1:]:
+        img = nb.load(roi)
+        assert img.shape == data.shape, "Mismatch in image shape"
+        assert np.allclose(img.affine, affine), "Mismatch in affine"
+        roi_data = np.asanyarray(img.dataobj)
+        roi_nonzero = np.any(roi_data, axis=3)
+        assert not np.any(roi_nonzero & nonzero), "Overlapping ROIs"
+        nonzero |= roi_nonzero
+        data += roi_data
+        del roi_data
+
+    if newpath is None:
+        newpath = Path()
+    out_file = str((Path(newpath) / "combined.nii.gz").absolute())
+    img.__class__(data, affine, header).to_filename(out_file)
+    return out_file
diff --git a/nibabies/interfaces/workbench.py b/nibabies/interfaces/workbench.py
index edb8d77..5ac70bc 100644
--- a/nibabies/interfaces/workbench.py
+++ b/nibabies/interfaces/workbench.py
@@ -1,7 +1,12 @@
+import os
 from nipype.interfaces.base import CommandLineInputSpec, File, traits, TraitedSpec, Str
-from nipype.interfaces.base.traits_extension import InputMultiObject
+from nipype.interfaces.base.traits_extension import InputMultiObject, OutputMultiObject, isdefined
 from nipype.interfaces.workbench.base import WBCommand
-
+# patch
+from nipype.interfaces.workbench.cifti import (
+    CiftiSmoothInputSpec as _CiftiSmoothInputSpec,
+    CiftiSmooth as _CiftiSmooth
+)
 
 VALID_STRUCTURES = (
     "CORTEX_LEFT",
@@ -50,6 +55,8 @@ class CiftiCreateDenseFromTemplateInputSpec(CommandLineInputSpec):
     )
     out_file = File(
         name_source=["in_file"],
+        name_template="template_%s.nii",
+        keep_extension=True,
         argstr="%s",
         position=1,
         desc="The output CIFTI file",
@@ -155,19 +162,19 @@ class CiftiCreateDenseFromTemplate(WBCommand):
     >>> from nibabies.interfaces import workbench as wb
     >>> frmtpl = wb.CiftiCreateDenseFromTemplate()
     >>> frmtpl.inputs.in_file = data_dir / "func.dtseries.nii"
-    >>> frmtpl.inputs.out_file = "out.dtseries.nii"
     >>> frmtpl.inputs.series = True
     >>> frmtpl.inputs.series_step = 0.8
     >>> frmtpl.inputs.series_start = 0
     >>> frmtpl.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
-    'wb_command -cifti-create-dense-from-template .../func.dtseries.nii out.dtseries.nii \
-    -series 0.8 0.0'
+    'wb_command -cifti-create-dense-from-template .../func.dtseries.nii \
+    template_func.dtseries.nii -series 0.8 0.0'
 
     >>> frmtpl.inputs.volume = [("OTHER", data_dir / 'functional.nii', True), \
         ("PUTAMEN_LEFT", data_dir / 'functional.nii')]
     >>> frmtpl.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
-    'wb_command -cifti-create-dense-from-template .../func.dtseries.nii out.dtseries.nii \
-    -series 0.8 0.0 -volume OTHER .../functional.nii -from-cropped \
+    'wb_command -cifti-create-dense-from-template .../func.dtseries.nii \
+    template_func.dtseries.nii -series 0.8 0.0 \
+    -volume OTHER .../functional.nii -from-cropped \
     -volume PUTAMEN_LEFT .../functional.nii'
     """
 
@@ -177,36 +184,35 @@ class CiftiCreateDenseFromTemplate(WBCommand):
 
     def _format_arg(self, name, trait_spec, value):
         if name in ("metric", "label", "volume"):
-            argstr = ""
+            cmds = []
             for val in value:
                 if val[-1] is True:  # volume specific
-                    val = val[:2] + ("-from-cropped ",)
-                argstr += " ".join((f"-{name}",) + val)
-            return trait_spec.argstr % argstr
+                    val = val[:2] + ("-from-cropped",)
+                cmds.append(" ".join((f"-{name}",) + val))
+            return trait_spec.argstr % " ".join(cmds)
         return super()._format_arg(name, trait_spec, value)
 
 
 class CiftiCreateDenseTimeseriesInputSpec(CommandLineInputSpec):
     out_file = File(
-        name_source=["in_file"],
-        name_template="%s.dtseries.nii",
-        keep_extension=False,
+        value='out.dtseries.nii',
+        usedefault=True,
         argstr="%s",
         position=0,
         desc="The output CIFTI file",
     )
-    in_file = File(
+    volume_data = File(
         exists=True,
-        mandatory=True,
         argstr="-volume %s",
         position=1,
+        requires=["volume_structure_labels"],
         desc="volume file containing all voxel data for all volume structures",
     )
-    structure_label_volume = File(
+    volume_structure_labels = File(
         exists=True,
-        mandatory=True,
         argstr="%s",
         position=2,
+        requires=["volume_data"],
         desc="label volume file containing labels for cifti structures",
     )
     left_metric = File(
@@ -248,25 +254,23 @@ class CiftiCreateDenseTimeseriesInputSpec(CommandLineInputSpec):
         requires=["cerebellum_metric"],
         desc="ROI (as metric file) of vertices to use from cerebellum",
     )
-    timestep = traits.Float(
-        1.0,
-        usedefault=True,
-        argstr="-timestep %g",
-        desc="the timestep, in seconds",
-    )
     timestart = traits.Float(
-        0.0,
-        usedefault=True,
         argstr="-timestart %g",
+        position=9,
         desc="the time at the first frame, in seconds",
     )
+    timestep = traits.Float(
+        argstr="-timestep %g",
+        position=10,
+        desc="the timestep, in seconds",
+    )
     unit = traits.Enum(
         "SECOND",
         "HERTZ",
         "METER",
         "RADIAN",
-        usedefault=True,
         argstr="-unit %s",
+        position=11,
         desc="use a unit other than time",
     )
 
@@ -323,121 +327,30 @@ class CiftiCreateDenseTimeseries(WBCommand):
 
     >>> from nibabies.interfaces.workbench import CiftiCreateDenseTimeseries
     >>> createdts = CiftiCreateDenseTimeseries()
-    >>> createdts.inputs.in_file = data_dir /'functional.nii'
-    >>> createdts.inputs.structure_label_volume = data_dir /'atlas.nii'
+    >>> createdts.inputs.volume_data = data_dir /'functional.nii'
+    >>> createdts.inputs.volume_structure_labels = data_dir / 'atlas.nii'
     >>> createdts.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
-    'wb_command -cifti-create-dense-timeseries functional.dtseries.nii \
-    -volume .../functional.nii .../atlas.nii -timestart 0 -timestep 1 -unit SECOND'
+    'wb_command -cifti-create-dense-timeseries out.dtseries.nii \
+    -volume .../functional.nii .../atlas.nii'
     """
 
     input_spec = CiftiCreateDenseTimeseriesInputSpec
     output_spec = CiftiCreateDenseTimeseriesOutputSpec
     _cmd = "wb_command -cifti-create-dense-timeseries"
 
-
-class CiftiDilateInputSpec(CommandLineInputSpec):
-    in_file = File(
-        exists=True,
-        mandatory=True,
-        argstr="%s",
-        position=0,
-        desc="The input CIFTI file",
-    )
-    direction = traits.Enum(
-        "ROW",
-        "COLUMN",
-        mandatory=True,
-        argstr="%s",
-        position=1,
-        desc="Which dimension to dilate along, ROW or COLUMN",
-    )
-    surface_distance = traits.Int(
-        mandatory=True,
-        argstr="%d",
-        position=2,
-        desc="The distance to dilate on surfaces, in mm",
-    )
-    volume_distance = traits.Int(
-        mandatory=True,
-        argstr="%d",
-        position=3,
-        desc="The distance to dilate in the volume, in mm",
-    )
-    out_file = File(
-        name_source=["in_file"],
-        name_template="dilated_%s.nii",
-        keep_extension=True,
-        argstr="%s",
-        position=4,
-        desc="The dilated CIFTI file",
-    )
-    left_surface = File(
-        exists=True,
-        position=5,
-        argstr="-left-surface %s",
-        desc="Specify the left surface to use",
-    )
-    left_corrected_areas = File(
-        exists=True,
-        position=6,
-        requires=["left_surface"],
-        argstr="-left-corrected-areas %s",
-        desc="vertex areas (as a metric) to use instead of computing them from the left surface.",
-    )
-    right_surface = File(
-        exists=True,
-        position=7,
-        argstr="-right-surface %s",
-        desc="Specify the right surface to use",
-    )
-    right_corrected_areas = File(
-        exists=True,
-        position=8,
-        requires=["right_surface"],
-        argstr="-right-corrected-areas %s",
-        desc="vertex areas (as a metric) to use instead of computing them from the right surface",
-    )
-    cerebellum_surface = File(
-        exists=True,
-        position=9,
-        argstr="-cerebellum-surface %s",
-        desc="specify the cerebellum surface to use",
-    )
-    cerebellum_corrected_areas = File(
-        exists=True,
-        position=10,
-        requires=["cerebellum_surface"],
-        argstr="-cerebellum-corrected-areas %s",
-        desc="vertex areas (as a metric) to use instead of computing them from the cerebellum "
-        "surface",
-    )
-    bad_brainordinate_roi = File(
-        exists=True,
-        position=11,
-        argstr="-bad-brainordinate-roi %s",
-        desc="CIFTI dscalar or dtseries file, positive values denote brainordinates to have their "
-        "values replaced",
-    )
-    nearest = traits.Bool(
-        position=12,
-        argstr="-nearest",
-        desc="Use nearest good value instead of a weighted average",
-    )
-    merged_volume = traits.Bool(
-        position=13,
-        argstr="-merged-volume",
-        desc="treat volume components as if they were a single component",
-    )
-    legacy_mode = traits.Bool(
-        position=14,
-        argstr="-legacy-mode",
-        desc="Use the math from v1.3.2 and earlier for weighted dilation",
-    )
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = os.path.abspath(self.inputs.out_file)
+        return outputs
 
 
 class CiftiCreateLabelInputSpec(CommandLineInputSpec):
     out_file = File(
-        mandatory=True, argstr="%s", position=0, desc="the output CIFTI file"
+        value="out.dlabel.nii",
+        usedefault=True,
+        argstr="%s",
+        position=0,
+        desc="the output CIFTI file",
     )
     volume_label = File(
         exists=True,
@@ -551,7 +464,6 @@ class CiftiCreateLabel(WBCommand):
 
     >>> from nibabies.interfaces import workbench as wb
     >>> lab = wb.CiftiCreateLabel()
-    >>> lab.inputs.out_file = "out.dlabel.nii"
     >>> lab.inputs.volume_label = data_dir / "functional.nii"
     >>> lab.inputs.structure_label_volume = data_dir / "functional.nii"
     >>> lab.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
@@ -562,6 +474,111 @@ class CiftiCreateLabel(WBCommand):
     output_spec = CiftiCreateLabelOutputSpec
     _cmd = "wb_command -cifti-create-label"
 
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = os.path.abspath(self.inputs.out_file)
+        return outputs
+
+
+class CiftiDilateInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=0,
+        desc="The input CIFTI file",
+    )
+    direction = traits.Enum(
+        "ROW",
+        "COLUMN",
+        mandatory=True,
+        argstr="%s",
+        position=1,
+        desc="Which dimension to dilate along, ROW or COLUMN",
+    )
+    surface_distance = traits.Int(
+        mandatory=True,
+        argstr="%d",
+        position=2,
+        desc="The distance to dilate on surfaces, in mm",
+    )
+    volume_distance = traits.Int(
+        mandatory=True,
+        argstr="%d",
+        position=3,
+        desc="The distance to dilate in the volume, in mm",
+    )
+    out_file = File(
+        name_source=["in_file"],
+        name_template="dilated_%s.nii",
+        keep_extension=True,
+        argstr="%s",
+        position=4,
+        desc="The dilated CIFTI file",
+    )
+    left_surface = File(
+        exists=True,
+        position=5,
+        argstr="-left-surface %s",
+        desc="Specify the left surface to use",
+    )
+    left_corrected_areas = File(
+        exists=True,
+        position=6,
+        requires=["left_surface"],
+        argstr="-left-corrected-areas %s",
+        desc="vertex areas (as a metric) to use instead of computing them from the left surface.",
+    )
+    right_surface = File(
+        exists=True,
+        position=7,
+        argstr="-right-surface %s",
+        desc="Specify the right surface to use",
+    )
+    right_corrected_areas = File(
+        exists=True,
+        position=8,
+        requires=["right_surface"],
+        argstr="-right-corrected-areas %s",
+        desc="vertex areas (as a metric) to use instead of computing them from the right surface",
+    )
+    cerebellum_surface = File(
+        exists=True,
+        position=9,
+        argstr="-cerebellum-surface %s",
+        desc="specify the cerebellum surface to use",
+    )
+    cerebellum_corrected_areas = File(
+        exists=True,
+        position=10,
+        requires=["cerebellum_surface"],
+        argstr="-cerebellum-corrected-areas %s",
+        desc="vertex areas (as a metric) to use instead of computing them from the cerebellum "
+        "surface",
+    )
+    bad_brainordinate_roi = File(
+        exists=True,
+        position=11,
+        argstr="-bad-brainordinate-roi %s",
+        desc="CIFTI dscalar or dtseries file, positive values denote brainordinates to have their "
+        "values replaced",
+    )
+    nearest = traits.Bool(
+        position=12,
+        argstr="-nearest",
+        desc="Use nearest good value instead of a weighted average",
+    )
+    merged_volume = traits.Bool(
+        position=13,
+        argstr="-merged-volume",
+        desc="treat volume components as if they were a single component",
+    )
+    legacy_mode = traits.Bool(
+        position=14,
+        argstr="-legacy-mode",
+        desc="Use the math from v1.3.2 and earlier for weighted dilation",
+    )
+
 
 class CiftiDilateOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc="Dilated CIFTI file")
@@ -938,6 +955,193 @@ class CiftiResample(WBCommand):
     _cmd = "wb_command -cifti-resample"
 
 
+class CiftiSeparateInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=0,
+        desc="the cifti to separate a component of",
+    )
+    direction = traits.Enum(
+        "ROW",
+        "COLUMN",
+        mandatory=True,
+        argstr="%s",
+        position=1,
+        desc="which dimension to smooth along, ROW or COLUMN",
+    )
+    volume_all_file = File(
+        argstr="-volume-all %s",
+        position=2,
+        desc="separate all volume structures into a volume file",
+    )
+    volume_all_roi_file = File(
+        argstr="-roi %s",
+        position=3,
+        requires=["volume_all_file"],
+        desc="output the roi of which voxels have data",
+    )
+    volume_all_label_file = File(
+        argstr="-label %s",
+        position=4,
+        requires=["volume_all_file"],
+        desc="output a volume label file indicating the location of structures",
+    )
+    volume_all_crop = traits.Bool(
+        argstr="-crop",
+        position=5,
+        requires=["volume_all_file"],
+        desc="crop volume to the size of the data rather than using the original volume size",
+    )
+    # the following can be repeated
+    label = InputMultiObject(
+        traits.Either(
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File()),
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File(), File()),
+        ),
+        argstr="%s",
+        position=6,
+        desc="separate one or more surface models into a surface label file",
+    )
+    metric = InputMultiObject(
+        traits.Either(
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File()),
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File(), File()),  # -roi
+        ),
+        argstr="%s",
+        position=7,
+        desc="separate one or more surface models into a metric file",
+    )
+    volume = InputMultiObject(
+        traits.Either(
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File()),
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File(), File()),  # -roi
+            traits.Tuple(traits.Enum(VALID_STRUCTURES, File(), traits.Bool)),  # -crop
+            traits.Tuple(traits.Enum(VALID_STRUCTURES), File(), File(), traits.Bool),  # -roi -crop
+        ),
+        argstr="%s",
+        position=8,
+        desc="separate one or more volume structures into a volume file",
+    )
+
+
+class CiftiSeparateOutputSpec(TraitedSpec):
+    volume_all_file = File(desc="File containing all volume structures")
+    volume_all_roi_file = File(desc="Output the roi of which voxels have data")
+    volume_all_label_file = File(
+        desc="output a volume label file indicating the location of structures"
+    )
+    label_files = OutputMultiObject(File(), desc="Output label files")
+    label_roi_files = OutputMultiObject(File(), desc="Output label rois files")
+    metric_files = OutputMultiObject(File(), desc="Output metric files")
+    metric_roi_files = OutputMultiObject(File(), desc="Output metric rois files")
+    volume_files = OutputMultiObject(File(), desc="Output volume files")
+    volume_roi_files = OutputMultiObject(File(), desc="Output volume roi files")
+
+
+class CiftiSeparate(WBCommand):
+    """
+    Extract left or right hemisphere surface from CIFTI file (.dtseries)
+    other structure can also be extracted
+    The input cifti file must have a brain models mapping on the chosen
+    dimension, columns for .dtseries.
+
+    >>> separate = CiftiSeparate()
+    >>> separate.inputs.in_file = data_dir / "func.dtseries.nii"
+    >>> separate.inputs.direction = "COLUMN"
+    >>> separate.inputs.volume_all_file = "volume_all.nii.gz"
+    >>> separate.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
+    'wb_command -cifti-separate .../func.dtseries.nii COLUMN \
+    -volume-all volume_all.nii.gz'
+
+    Metrics, labels, and volumes can also be freely extracted
+    >>> separate.inputs.metric = [("CORTEX_LEFT", "cortexleft.func.gii")]
+    >>> separate.inputs.volume = [("HIPPOCAMPUS_LEFT", "hippoL.nii.gz"), \
+        ("HIPPOCAMPUS_RIGHT", "hippoR.nii.gz", "hippoR.roi.nii.gz")]
+    >>> separate.cmdline  #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
+    'wb_command -cifti-separate .../func.dtseries.nii COLUMN \
+    -volume-all volume_all.nii.gz -metric CORTEX_LEFT cortexleft.func.gii \
+    -volume HIPPOCAMPUS_LEFT hippoL.nii.gz \
+    -volume HIPPOCAMPUS_RIGHT hippoR.nii.gz -roi hippoR.roi.nii.gz'
+
+    """
+    input_spec = CiftiSeparateInputSpec
+    output_spec = CiftiSeparateOutputSpec
+    _cmd = "wb_command -cifti-separate"
+    _label_roi_files = []
+    _metric_roi_files = []
+    _volume_roi_files = []
+
+    def _format_arg(self, name, trait_spec, value):
+        if name in ("label", "metric", "volume"):
+            cmds = []
+            for i, val in enumerate(value):
+                if len(val) == 3:
+                    if val[-1] is True:
+                        val = val[:-1] + ("-crop",)
+                    else:
+                        val = val[:-1] + ("-roi", val[-1])
+                        self._set_roi_file(name, val[-1])
+                elif len(val) == 4:
+                    val = val[:-2] + ("-roi", val[-2]) + ("crop") if val[-1] is True else ()
+                    self._set_roi_file(name, val[-2])
+                cmds.append(" ".join((f"-{name}",) + val))
+            return trait_spec.argstr % " ".join(cmds)
+        return super()._format_arg(name, trait_spec, value)
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        if self.inputs.volume_all_file:
+            outputs["volume_all_file"] = os.path.abspath(self.inputs.volume_all_file)
+        if self.inputs.volume_all_roi_file:
+            outputs["volume_all_roi_file"] = os.path.abspath(self.inputs.volume_all_roi_file)
+        if self.inputs.volume_all_label_file:
+            outputs["volume_all_label_file"] = os.path.abspath(self.inputs.volume_all_label_file)
+        if self.inputs.label:
+            for label in self.inputs.label:
+                outputs["label_files"] = (outputs["label_files"] or []) + \
+                    self._gen_filename(label[2])
+            if self._label_roi_files:
+                outputs["label_roi_files"] = self._label_roi_files
+        if self.inputs.metric:
+            for metric in self.inputs.metric:
+                outputs["metric_files"] = (outputs["metric_files"] or []) + \
+                    self._gen_filename(metric[2])
+            if self._metric_roi_files:
+                outputs["metric_roi_files"] = self._metric_roi_files
+        if self.inputs.volume:
+            for volume in self.inputs.volume:
+                outputs["volume_files"] = (outputs["volume_files"] or []) + \
+                    self._gen_filename(volume[2])
+            if self._volume_roi_files:
+                outputs["volume_roi_files"] = self._volume_roi_files
+        return outputs
+
+    def _set_roi_file(self, name, file):
+        rois = getattr(self, f"_{name}_roi_files")
+        rois.append(self._gen_filename(file))
+
+
+class CiftiSmoothInputSpec(_CiftiSmoothInputSpec):
+    left_surf = File(
+        exists=True,
+        position=5,
+        argstr="-left-surface %s",
+        desc="Specify the left surface to use",
+    )
+    right_surf = File(
+        exists=True,
+        position=7,
+        argstr="-right-surface %s",
+        desc="Specify the right surface to use",
+    )
+
+
+class CiftiSmooth(_CiftiSmooth):
+    input_spec = CiftiSmoothInputSpec
+
+
 class VolumeAffineResampleInputSpec(CommandLineInputSpec):
     in_file = File(
         exists=True,
diff --git a/nibabies/utils/misc.py b/nibabies/utils/misc.py
index ee727a6..a61fc82 100644
--- a/nibabies/utils/misc.py
+++ b/nibabies/utils/misc.py
@@ -41,9 +41,29 @@ def check_deps(workflow):
 
 
 def cohort_by_months(template, months):
+    """
+    Produce a recommended cohort based on partipants age
+    """
     cohort_key = {
-        'MNIInfant': (2, 5, 8, 11, 14, 17, 21, 27, 33, 44, 60),
-        'UNCInfant': (8, 12, 24),
+        'MNIInfant': (
+            # upper bound of template | cohort
+            2,  # 1
+            5,  # 2
+            8,  # 3
+            11,  # 4
+            14,  # 5
+            17,  # 6
+            21,  # 7
+            27,  # 8
+            33,  # 9
+            44,  # 10
+            60,  # 11
+        ),
+        'UNCInfant': (
+            8,  # 1
+            12,  # 2
+            24,  # 3
+        ),
     }
     ages = cohort_key.get(template)
     if ages is None:
@@ -52,7 +72,6 @@ def cohort_by_months(template, months):
     for cohort, age in enumerate(ages, 1):
         if months <= age:
             return cohort
-
     raise KeyError("Age exceeds all cohorts!")
 
 
diff --git a/nibabies/workflows/bold/alignment.py b/nibabies/workflows/bold/alignment.py
new file mode 100644
index 0000000..1f06a49
--- /dev/null
+++ b/nibabies/workflows/bold/alignment.py
@@ -0,0 +1,269 @@
+"""
+Subcortical alignment into MNI space
+"""
+
+from nibabies.interfaces.nibabel import MergeROIs
+
+
+def init_subcortical_mni_alignment_wf(*, vol_sigma=0.8, name='subcortical_mni_alignment_wf'):
+    """
+    Align individual subcortical structures into MNI space.
+
+    This is a nipype workflow port of the DCAN infant pipeline:
+    https://github.com/DCAN-Labs/dcan-infant-pipeline/blob\
+    /master/fMRISurface/scripts/SubcorticalAlign_ROIs.sh
+
+
+    Parameters
+    ----------
+    name : :obj:`str`
+        Name of the workflow
+    vol_sigma : :obj:`float`
+        The sigma for the gaussian volume smoothing kernel, in mm
+
+    Inputs
+    ------
+    bold_file : :obj:`str`
+        BOLD file
+    bold_roi : :obj:`str`
+        File containing ROIs in BOLD space
+    atlas_roi : :obj:`str`
+        File containing ROIs in atlas space
+    std_xfm : :obj:`str`
+        File containing transform to the standard (MNI) space
+
+    Outputs
+    -------
+    subcortical_file : :obj:`str`
+        Volume file containing all ROIs individually aligned to standard
+    """
+    from nipype.pipeline import engine as pe
+    from nipype.interfaces import utility as niu, fsl
+    from ...interfaces.workbench import (
+        CiftiCreateDenseTimeseries,
+        CiftiCreateLabel,
+        CiftiDilate,
+        CiftiResample,
+        CiftiSeparate,
+        CiftiSmooth,
+        VolumeAffineResample,
+        VolumeAllLabelsToROIs,
+        VolumeLabelExportTable,
+        VolumeLabelImport,
+    )
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=["bold_file", "bold_roi", "atlas_roi", "atlas_xfm"]),
+        name="inputnode",
+    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=["subcortical_file"]), name='outputnode')
+
+    applyxfm_atlas = pe.Node(fsl.FLIRT(), name="applyxfm_atlas")
+    vol_resample = pe.Node(
+        VolumeAffineResample(method="ENCLOSING_VOXEL", flirt=True),
+        name="vol_resample"
+    )
+    subj_rois = pe.Node(VolumeAllLabelsToROIs(label_map=1), name="subj_rois")
+    split_rois = pe.Node(fsl.Split(dimension="t"), name="split_rois")
+    atlas_rois = pe.Node(VolumeAllLabelsToROIs(label_map=1), name="atlas_rois")
+    split_atlas_rois = pe.Node(fsl.Split(dimension="t"), name="split_atlas_rois")
+    atlas_labels = pe.Node(VolumeLabelExportTable(label_map=1), name="atlas_labels")
+    parse_labels = pe.Node(
+        niu.Function(function=parse_roi_labels, output_names=["structures", "label_ids"]),
+        name="parse_labels",
+    )
+
+    # The following is wrapped in a for-loop, iterating across each roi
+    # Instead, we will use MapNodes and iter across the varying inputs
+    roi2atlas = pe.MapNode(
+        fsl.FLIRT(
+            searchr_x=[-20, 20],
+            searchr_y=[-20, 20],
+            searchr_z=[-20, 20],
+            interp="nearestneighbour",
+        ),
+        name="roi2atlas",
+        iterfield=["in_file", "reference"],
+    )
+    applyxfm_roi = pe.MapNode(
+        fsl.ApplyXFM(interp="spline"),
+        iterfield=["reference", "in_matrix_file"],
+        name='applyxfm_roi',
+        mem_gb=4,
+    )
+    bold_mask_roi = pe.MapNode(
+        fsl.ApplyMask(),
+        iterfield=["in_file", "mask_file"],
+        name='bold_mask_roi',
+    )
+    mul_roi = pe.MapNode(
+        fsl.BinaryMaths(operation="mul"),
+        iterfield=["in_file", "operand_value"],
+        name='mul_roi',
+    )
+    mul_atlas_roi = pe.MapNode(
+        fsl.BinaryMaths(operation="mul"),
+        iterfield=["in_file", "operand_value"],
+        name='mul_atlas_roi',
+    )
+    vol_label = pe.MapNode(
+        VolumeLabelImport(drop_unused_labels=True),
+        iterfield=["in_file"],
+        name='vol_label',
+    )
+    vol_atlas_label = pe.MapNode(
+        VolumeLabelImport(drop_unused_labels=True),
+        iterfield=["in_file"],
+        name='vol_atlas_label',
+    )
+    create_dtseries = pe.MapNode(
+        CiftiCreateDenseTimeseries(),
+        iterfield=["volume_data", "volume_structure_labels"],
+        name='create_dtseries',
+    )
+    create_label = pe.MapNode(
+        CiftiCreateLabel(),
+        iterfield=["volume_label", "structure_label_volume"],
+        name='create_label',
+    )
+    dilate = pe.MapNode(
+        CiftiDilate(direction="COLUMN", surface_distance=0, volume_distance=10),
+        iterfield=["in_file"],
+        name="dilate"
+    )
+    resample = pe.MapNode(
+        CiftiResample(
+            direction="COLUMN",
+            template_direction="COLUMN",
+            surface_method="ADAP_BARY_AREA",
+            volume_method="CUBIC",
+            volume_predilate=10,
+        ),
+        iterfield=["in_file", "template"],
+        name='resample',
+    )
+    smooth = pe.MapNode(
+        CiftiSmooth(direction="COLUMN", fix_zeros_vol=True, sigma_surf=0, sigma_vol=vol_sigma),
+        iterfield=["in_file"],
+        name="smooth"
+    )
+    separate = pe.MapNode(
+        CiftiSeparate(direction="COLUMN", volume_all_file='volume_all.nii.gz'),
+        iterfield=["in_file"],
+        name="separate"
+    )
+    fmt_agg_rois = pe.Node(
+        niu.Function(
+            function=format_agg_rois,
+            output_names=["first_image", "op_files", "op_string"],
+        ),
+        name='fmt_agg_rois',
+    )
+    agg_rois = pe.Node(fsl.MultiImageMaths(), name='agg_rois')
+    merge_rois = pe.Node(MergeROIs(), name="merge_rois")
+
+    workflow = Workflow(name=name)
+    # fmt: off
+    workflow.connect([
+        (inputnode, applyxfm_atlas, [
+            ("bold_file", "in_file"),
+            ("atlas_roi", "reference"),
+            ("atlas_xfm", "in_matrix_file")]),
+        (inputnode, vol_resample, [
+            ("bold_roi", "in_file"),
+            ("atlas_xfm", "affine"),
+            ("bold_roi", "flirt_source_volume")]),
+        (applyxfm_atlas, vol_resample, [
+            ("out_file", "volume_space"),
+            ("out_file", "flirt_target_volume")]),
+        (inputnode, subj_rois, [("bold_roi", "in_file")]),
+        (inputnode, atlas_rois, [("atlas_roi", "in_file")]),
+        (subj_rois, split_rois, [("out_file", "in_file")]),
+        (atlas_rois, split_atlas_rois, [("out_file", "in_file")]),
+        (inputnode, atlas_labels, [("atlas_roi", "in_file")]),
+        (atlas_labels, parse_labels, [("out_file", "label_file")]),
+        # for loop across ROIs
+        (split_rois, roi2atlas, [("out_files", "in_file")]),
+        (split_atlas_rois, roi2atlas, [("out_files", "reference")]),
+        (inputnode, applyxfm_roi, [("bold_file", "in_file")]),
+        (split_atlas_rois, applyxfm_roi, [("out_files", "reference")]),
+        (roi2atlas, applyxfm_roi, [("out_matrix_file", "in_matrix_file")]),
+        (applyxfm_roi, bold_mask_roi, [("out_file", "in_file")]),
+        (roi2atlas, bold_mask_roi, [("out_file", "mask_file")]),
+        (roi2atlas, mul_roi, [("out_file", "in_file")]),
+        (parse_labels, mul_roi, [("label_ids", "operand_value")]),
+        (split_atlas_rois, mul_atlas_roi, [("out_files", "in_file")]),
+        (parse_labels, mul_atlas_roi, [("label_ids", "operand_value")]),
+        (mul_roi, vol_label, [("out_file", "in_file")]),
+        (atlas_labels, vol_label, [("out_file", "label_list_file")]),
+        (mul_atlas_roi, vol_atlas_label, [("out_file", "in_file")]),
+        (atlas_labels, vol_atlas_label, [("out_file", "label_list_file")]),
+        (bold_mask_roi, create_dtseries, [("out_file", "volume_data")]),
+        (vol_label, create_dtseries, [("out_file", "volume_structure_labels")]),
+        (vol_atlas_label, create_label, [
+            ("out_file", "volume_label"),
+            ("out_file", "structure_label_volume")]),
+        (create_dtseries, dilate, [("out_file", "in_file")]),
+        (dilate, resample, [("out_file", "in_file")]),
+        (create_label, resample, [("out_file", "template")]),
+        (resample, smooth, [("out_file", "in_file")]),
+        (smooth, separate, [("out_file", "in_file")]),
+        # end loop
+        (mul_roi, fmt_agg_rois, [("out_file", "rois")]),
+        (fmt_agg_rois, agg_rois, [
+            ("first_image", "in_file"),
+            ("op_files", "operand_files"),
+            ("op_string", "op_string")]),
+        (separate, merge_rois, [("volume_all_file", "in_files")]),
+        (merge_rois, outputnode, [("out_file", "subcortical_file")]),
+    ])
+    # fmt: on
+    return workflow
+
+
+def parse_roi_labels(label_file):
+    """
+    Parse a label file composed of one or more sets of:
+    <labelname>
+    <key> <red> <green> <blue> <alpha>
+
+    Return a list of structure names and label keys.
+
+    >>> structs, ids = parse_roi_labels(test_data / "labelfile.txt")
+    >>> structs
+    ['CEREBELLUM_LEFT', 'THALAMUS_LEFT', 'CAUDATE_LEFT']
+    >>> ids
+    [8, 10, 11]
+    """
+
+    with open(label_file) as fp:
+        lines = fp.readlines()
+    if len(lines) % 2 == 1:
+        raise RuntimeError("Label file is incomplete or invalid")
+    structs, label_ids = [], []
+    for idx, line in enumerate(lines):
+        if idx % 2 == 0:
+            structs.append(line.strip())
+        else:
+            label_ids.append(int(line.split(' ', 1)[0]))
+    return structs, label_ids
+
+
+def format_agg_rois(rois):
+    """
+    Helper function to format MultiImageMaths command.
+
+    Parameters
+    ----------
+    rois : `list` of `str`s
+        List of files
+
+    Returns
+    -------
+    first_image
+    op_files
+    op_string
+
+    """
+    return rois[0], rois[1:], ("-add %s " * (len(rois) - 1)).strip()
diff --git a/scripts/bold_subcortical.py b/scripts/bold_subcortical.py
new file mode 100644
index 0000000..2d50c77
--- /dev/null
+++ b/scripts/bold_subcortical.py
@@ -0,0 +1,66 @@
+"""Script for testing the subcortical MNI alignment"""
+from pathlib import Path
+
+
+def init_workflow(bold_file, bold_roi, bold_atlas_roi, atlas_xfm, vol_sigma):
+    from nibabies.workflows.bold.alignment import init_subcortical_mni_alignment_wf
+
+    wf = init_subcortical_mni_alignment_wf(vol_sigma=vol_sigma)
+    wf.inputs.inputnode.bold_file = bold_file
+    wf.inputs.inputnode.bold_roi = bold_roi
+    wf.inputs.inputnode.atlas_roi = bold_atlas_roi
+    wf.inputs.inputnode.atlas_xfm = atlas_xfm
+
+    wf.base_dir = Path('workdir').absolute()
+    return wf
+
+
+if __name__ == "__main__":
+    from argparse import ArgumentParser, RawTextHelpFormatter
+
+    parser = ArgumentParser(
+        description="DCAN subcortical MNI alignment",
+        formatter_class=RawTextHelpFormatter,
+    )
+    parser.add_argument(
+        "bold_file",
+        type=Path,
+        help="the input BOLD file",
+    )
+    parser.add_argument(
+        "bold_roi",
+        type=Path,
+        help="segmentations in BOLD space",
+    )
+    parser.add_argument(
+        "bold_atlas_roi",
+        type=Path,
+        help="segmentations in ROI space, unrefined",
+    )
+    parser.add_argument(
+        "atlas_xfm",
+        type=Path,
+        help="transformation of input BOLD file to MNI space",
+    )
+    parser.add_argument(
+        "--vol-sigma",
+        type=float,
+        default=0.8,
+        help="The sigma for the gaussian volume smoothing kernel, in mm",
+    )
+    parser.add_argument(
+        "--nipype-plugin",
+        default="MultiProc",
+        help="Nipype plugin to run workflow with",
+    )
+    opts = parser.parse_args()
+    wf = init_workflow(
+        opts.bold_file.absolute(),
+        opts.bold_roi.absolute(),
+        opts.bold_atlas_roi.absolute(),
+        opts.atlas_xfm.absolute(),
+        vol_sigma=opts.vol_sigma,
+    )
+
+    wf.config['execution']['crashfile_format'] = 'txt'
+    wf.run(plugin=opts.nipype_plugin)


Test patch:
diff --git a/nibabies/interfaces/conftest.py b/nibabies/conftest.py
similarity index 85%
rename from nibabies/interfaces/conftest.py
rename to nibabies/conftest.py
index e3c3caf..096df92 100644
--- a/nibabies/interfaces/conftest.py
+++ b/nibabies/conftest.py
@@ -1,8 +1,9 @@
 """py.test configuration"""
 from pathlib import Path
-import pytest
 from tempfile import TemporaryDirectory
+from pkg_resources import resource_filename
 
+import pytest
 
 FILES = (
     'functional.nii',
@@ -32,3 +33,4 @@ def data_dir():
 @pytest.fixture(autouse=True)
 def set_namespace(doctest_namespace, data_dir):
     doctest_namespace["data_dir"] = data_dir
+    doctest_namespace["test_data"] = Path(resource_filename("nibabies", "tests/data"))
diff --git a/nibabies/interfaces/tests/__init__.py b/nibabies/interfaces/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/nibabies/interfaces/tests/test_nibabel.py b/nibabies/interfaces/tests/test_nibabel.py
new file mode 100644
index 0000000..4b40d38
--- /dev/null
+++ b/nibabies/interfaces/tests/test_nibabel.py
@@ -0,0 +1,73 @@
+import uuid
+
+import nibabel as nb
+import numpy as np
+import pytest
+
+from ..nibabel import MergeROIs
+
+
+@pytest.fixture
+def create_roi(tmp_path):
+    files = []
+
+    def _create_roi(affine, img_data, roi_index):
+        img_data[tuple(roi_index)] = 1
+        nii = nb.Nifti1Image(img_data, affine)
+        filename = tmp_path / f"{str(uuid.uuid4())}.nii.gz"
+        files.append(filename)
+        nii.to_filename(filename)
+        return filename
+
+    yield _create_roi
+
+    for f in files:
+        f.unlink()
+
+
+# create a slightly off affine
+bad_affine = np.eye(4)
+bad_affine[0, -1] = -1
+
+
+@pytest.mark.parametrize(
+    "affine, data, roi_index, error, err_message",
+    [
+        (np.eye(4), np.zeros((2, 2, 2, 2), dtype=int), [1, 0], None, None),
+        (
+            np.eye(4),
+            np.zeros((2, 2, 3, 2), dtype=int),
+            [1, 0],
+            True,
+            "Mismatch in image shape",
+        ),
+        (
+            bad_affine,
+            np.zeros((2, 2, 2, 2), dtype=int),
+            [1, 0],
+            True,
+            "Mismatch in affine",
+        ),
+        (
+            np.eye(4),
+            np.zeros((2, 2, 2, 2), dtype=int),
+            [0, 0, 0],
+            True,
+            "Overlapping ROIs",
+        ),
+    ],
+)
+def test_merge_rois(tmpdir, create_roi, affine, data, roi_index, error, err_message):
+    tmpdir.chdir()
+    roi0 = create_roi(np.eye(4), np.zeros((2, 2, 2, 2), dtype=int), [0, 0])
+    roi1 = create_roi(np.eye(4), np.zeros((2, 2, 2, 2), dtype=int), [0, 1])
+    test_roi = create_roi(affine, data, roi_index)
+
+    merge = MergeROIs(in_files=[roi0, roi1, test_roi])
+    if error is None:
+        merge.run()
+        return
+    # otherwise check expected exceptions
+    with pytest.raises(AssertionError) as err:
+        merge.run()
+    assert err_message in str(err.value)
diff --git a/nibabies/tests/data/labelfile.txt b/nibabies/tests/data/labelfile.txt
new file mode 100644
index 0000000..6a990b5
--- /dev/null
+++ b/nibabies/tests/data/labelfile.txt
@@ -0,0 +1,6 @@
+CEREBELLUM_LEFT
+8 230 148 34 255
+THALAMUS_LEFT
+10 0 118 14 255
+CAUDATE_LEFT
+11 122 186 220 255