Repository: PennLINC/xcp_d

Issue:
Processing native-resolution runs with different voxel sizes will fail
## Summary
<!--What is the nature of the bug?-->
When processing multiple runs in the same standard space (MNI152NLin2009cAsym) with the native BOLD resolution, if the different runs have different native resolutions (in this case, 1.625 x 1.625 x 3.992 mm voxels in one run and 1.625 x 1.625 x 3.99 mm in another), then XCP-D will fail with the following error:

```
Traceback:
	Traceback (most recent call last):
	  File "/usr/local/miniconda/lib/python3.10/site-packages/nipype/interfaces/base/core.py", line 397, in run
	    runtime = self._run_interface(runtime)
	  File "/usr/local/miniconda/lib/python3.10/site-packages/xcp_d/interfaces/connectivity.py", line 105, in _run_interface
	    n_voxels_in_masked_parcels = sum_masker_masked.fit_transform(atlas_img_bin)
	  File "/usr/local/miniconda/lib/python3.10/site-packages/sklearn/utils/_set_output.py", line 273, in wrapped
	    data_to_wrap = f(self, X, *args, **kwargs)
	  File "/usr/local/miniconda/lib/python3.10/site-packages/nilearn/maskers/nifti_labels_masker.py", line 455, in fit_transform
	    return self.fit().transform(imgs, confounds=confounds,
	  File "/usr/local/miniconda/lib/python3.10/site-packages/nilearn/maskers/nifti_labels_masker.py", line 376, in fit
	    raise ValueError(
	ValueError: Regions and mask do not have the same affine.
	labels_img: /out/xcp_d/atlases/atlas-4S156Parcels/space-MNI152NLin2009cAsym_atlas-4S156Parcels_dseg.nii.gz
	mask_img: /fmriprep/sub-002/func/sub-002_task-rest_run-2_space-MNI152NLin2009cAsym_desc-brain_mask.nii.gz
```

This was originally brought up in [this NeuroStars post](https://neurostars.org/t/regions-and-mask-do-not-have-the-same-affine-while-running-xcp-d/28341).


Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) Report
Attention: Patch coverage is `96.00000%` with `2 lines` in your changes are missing coverage. Please review.
> Project coverage is 82.36%. Comparing base [(`b2d16c1`)](https://app.codecov.io/gh/PennLINC/xcp_d/commit/b2d16c17b006cdf756bf0cb0f6b08dcab06e0f1d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) to head [(`bb0186b`)](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).

| [Files](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) | Patch % | Lines |
|---|---|---|
| [xcp\_d/interfaces/bids.py](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC#diff-eGNwX2QvaW50ZXJmYWNlcy9iaWRzLnB5) | 95.83% | [0 Missing and 2 partials :warning: ](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #1075      +/-   ##
==========================================
+ Coverage   82.33%   82.36%   +0.03%     
==========================================
  Files          48       48              
  Lines        5897     5916      +19     
  Branches      796      798       +2     
==========================================
+ Hits         4855     4873      +18     
  Misses        839      839              
- Partials      203      204       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1075?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).


Patch:
diff --git a/xcp_d/interfaces/bids.py b/xcp_d/interfaces/bids.py
index 06e1c18..bfdc4bc 100644
--- a/xcp_d/interfaces/bids.py
+++ b/xcp_d/interfaces/bids.py
@@ -1,8 +1,12 @@
 """Adapted interfaces from Niworkflows."""
 
+import os
+import shutil
 from json import loads
 from pathlib import Path
 
+import nibabel as nb
+import numpy as np
 from bids.layout import Config
 from nipype import logging
 from nipype.interfaces.base import (
@@ -16,6 +20,8 @@ from nipype.interfaces.base import (
 from niworkflows.interfaces.bids import DerivativesDataSink as BaseDerivativesDataSink
 from pkg_resources import resource_filename as pkgrf
 
+from xcp_d.utils.bids import get_entity
+
 # NOTE: Modified for xcpd's purposes
 xcp_d_spec = loads(Path(pkgrf("xcp_d", "data/xcp_d_bids_config.json")).read_text())
 bids_config = Config.load("bids")
@@ -190,3 +196,117 @@ class CollectRegistrationFiles(SimpleInterface):
             )
 
         return runtime
+
+
+class _CopyAtlasInputSpec(BaseInterfaceInputSpec):
+    name_source = traits.Str(
+        desc="The source file's name.",
+        mandatory=False,
+    )
+    in_file = File(
+        exists=True,
+        desc="The atlas file to copy.",
+        mandatory=True,
+    )
+    output_dir = Directory(
+        exists=True,
+        desc="The output directory.",
+        mandatory=True,
+    )
+    atlas = traits.Str(
+        desc="The atlas name.",
+        mandatory=True,
+    )
+
+
+class _CopyAtlasOutputSpec(TraitedSpec):
+    out_file = File(
+        exists=True,
+        desc="The copied atlas file.",
+    )
+
+
+class CopyAtlas(SimpleInterface):
+    """Copy atlas file to output directory.
+
+    Parameters
+    ----------
+    name_source : :obj:`str`
+        The source name of the atlas file.
+    in_file : :obj:`str`
+        The atlas file to copy.
+    output_dir : :obj:`str`
+        The output directory.
+    atlas : :obj:`str`
+        The name of the atlas.
+
+    Returns
+    -------
+    out_file : :obj:`str`
+        The path to the copied atlas file.
+
+    Notes
+    -----
+    I can't use DerivativesDataSink because it has a problem with dlabel CIFTI files.
+    It gives the following error:
+    "AttributeError: 'Cifti2Header' object has no attribute 'set_data_dtype'"
+
+    I can't override the CIFTI atlas's data dtype ahead of time because setting it to int8 or int16
+    somehow converts all of the values in the data array to weird floats.
+    This could be a version-specific nibabel issue.
+
+    I've also updated this function to handle JSON and TSV files as well.
+    """
+
+    input_spec = _CopyAtlasInputSpec
+    output_spec = _CopyAtlasOutputSpec
+
+    def _run_interface(self, runtime):
+        output_dir = self.inputs.output_dir
+        in_file = self.inputs.in_file
+        name_source = self.inputs.name_source
+        atlas = self.inputs.atlas
+
+        atlas_out_dir = os.path.join(output_dir, f"xcp_d/atlases/atlas-{atlas}")
+
+        if in_file.endswith(".json"):
+            out_basename = f"atlas-{atlas}_dseg.json"
+        elif in_file.endswith(".tsv"):
+            out_basename = f"atlas-{atlas}_dseg.tsv"
+        else:
+            extension = ".nii.gz" if name_source.endswith(".nii.gz") else ".dlabel.nii"
+            space = get_entity(name_source, "space")
+            res = get_entity(name_source, "res")
+            den = get_entity(name_source, "den")
+            cohort = get_entity(name_source, "cohort")
+
+            cohort_str = f"_cohort-{cohort}" if cohort else ""
+            res_str = f"_res-{res}" if res else ""
+            den_str = f"_den-{den}" if den else ""
+            if extension == ".dlabel.nii":
+                out_basename = f"space-{space}_atlas-{atlas}{den_str}{cohort_str}_dseg{extension}"
+            elif extension == ".nii.gz":
+                out_basename = f"space-{space}_atlas-{atlas}{res_str}{cohort_str}_dseg{extension}"
+
+        os.makedirs(atlas_out_dir, exist_ok=True)
+        out_file = os.path.join(atlas_out_dir, out_basename)
+
+        if out_file.endswith(".nii.gz") and os.path.isfile(out_file):
+            # Check that native-resolution atlas doesn't have a different resolution from the last
+            # run's atlas.
+            old_img = nb.load(out_file)
+            new_img = nb.load(in_file)
+            if not np.allclose(old_img.affine, new_img.affine):
+                raise ValueError(
+                    f"Existing '{atlas}' atlas affine ({out_file}) is different from the input "
+                    f"file affine ({in_file})."
+                )
+
+        # Don't copy the file if it exists, to prevent any race conditions between parallel
+        # processes.
+        if not os.path.isfile(out_file):
+            shutil.copyfile(in_file, out_file)
+
+        self._results["out_file"] = out_file
+
+        return runtime
diff --git a/xcp_d/utils/atlas.py b/xcp_d/utils/atlas.py
index 561973b..ac1b70d 100644
--- a/xcp_d/utils/atlas.py
+++ b/xcp_d/utils/atlas.py
@@ -155,68 +155,3 @@ def get_atlas_cifti(atlas):
         )
 
     return atlas_file, atlas_labels_file, atlas_metadata_file
-
-
-def copy_atlas(name_source, in_file, output_dir, atlas):
-    """Copy atlas file to output directory.
-
-    Parameters
-    ----------
-    name_source : :obj:`str`
-        The source name of the atlas file.
-    in_file : :obj:`str`
-        The atlas file to copy.
-    output_dir : :obj:`str`
-        The output directory.
-    atlas : :obj:`str`
-        The name of the atlas.
-
-    Returns
-    -------
-    out_file : :obj:`str`
-        The path to the copied atlas file.
-
-    Notes
-    -----
-    I can't use DerivativesDataSink because it has a problem with dlabel CIFTI files.
-    It gives the following error:
-    "AttributeError: 'Cifti2Header' object has no attribute 'set_data_dtype'"
-
-    I can't override the CIFTI atlas's data dtype ahead of time because setting it to int8 or int16
-    somehow converts all of the values in the data array to weird floats.
-    This could be a version-specific nibabel issue.
-
-    I've also updated this function to handle JSON and TSV files as well.
-    """
-    import os
-    import shutil
-
-    from xcp_d.utils.bids import get_entity
-
-    if in_file.endswith(".json"):
-        out_basename = f"atlas-{atlas}_dseg.json"
-    elif in_file.endswith(".tsv"):
-        out_basename = f"atlas-{atlas}_dseg.tsv"
-    else:
-        extension = ".nii.gz" if name_source.endswith(".nii.gz") else ".dlabel.nii"
-        space = get_entity(name_source, "space")
-        res = get_entity(name_source, "res")
-        den = get_entity(name_source, "den")
-        cohort = get_entity(name_source, "cohort")
-
-        cohort_str = f"_cohort-{cohort}" if cohort else ""
-        res_str = f"_res-{res}" if res else ""
-        den_str = f"_den-{den}" if den else ""
-        if extension == ".dlabel.nii":
-            out_basename = f"space-{space}_atlas-{atlas}{den_str}{cohort_str}_dseg{extension}"
-        elif extension == ".nii.gz":
-            out_basename = f"space-{space}_atlas-{atlas}{res_str}{cohort_str}_dseg{extension}"
-
-    atlas_out_dir = os.path.join(output_dir, f"xcp_d/atlases/atlas-{atlas}")
-    os.makedirs(atlas_out_dir, exist_ok=True)
-    out_file = os.path.join(atlas_out_dir, out_basename)
-    # Don't copy the file if it exists, to prevent any race conditions between parallel processes.
-    if not os.path.isfile(out_file):
-        shutil.copyfile(in_file, out_file)
-
-    return out_file
diff --git a/xcp_d/utils/bids.py b/xcp_d/utils/bids.py
index 8230b3c..0938f70 100644
--- a/xcp_d/utils/bids.py
+++ b/xcp_d/utils/bids.py
@@ -317,11 +317,12 @@ def collect_data(
     # This probably works well for resolution (1 typically means 1x1x1,
     # 2 typically means 2x2x2, etc.), but probably doesn't work well for density.
     resolutions = layout.get_res(**queries["bold"])
-    densities = layout.get_den(**queries["bold"])
-    if len(resolutions) > 1:
-        queries["bold"]["resolution"] = resolutions[0]
+    if len(resolutions) >= 1:
+        # This will also select res-* when there are both res-* and native-resolution files.
+        queries["bold"]["res"] = resolutions[0]
 
-    if len(densities) > 1:
+    densities = layout.get_den(**queries["bold"])
+    if len(densities) >= 1:
         queries["bold"]["den"] = densities[0]
 
     # Check for anatomical images, and determine if T2w xfms must be used.
diff --git a/xcp_d/workflows/connectivity.py b/xcp_d/workflows/connectivity.py
index 0535d49..643651a 100644
--- a/xcp_d/workflows/connectivity.py
+++ b/xcp_d/workflows/connectivity.py
@@ -7,7 +7,7 @@ from nipype.pipeline import engine as pe
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
 from xcp_d.interfaces.ants import ApplyTransforms
-from xcp_d.interfaces.bids import DerivativesDataSink
+from xcp_d.interfaces.bids import CopyAtlas, DerivativesDataSink
 from xcp_d.interfaces.connectivity import (
     CiftiConnect,
     CiftiParcellate,
@@ -21,12 +21,7 @@ from xcp_d.interfaces.workbench import (
     CiftiCreateDenseFromTemplate,
     CiftiParcellateWorkbench,
 )
-from xcp_d.utils.atlas import (
-    copy_atlas,
-    get_atlas_cifti,
-    get_atlas_nifti,
-    select_atlases,
-)
+from xcp_d.utils.atlas import get_atlas_cifti, get_atlas_nifti, select_atlases
 from xcp_d.utils.doc import fill_doc
 from xcp_d.utils.utils import get_std2bold_xfms
 
@@ -214,21 +209,11 @@ def init_load_atlases_wf(
         ])  # fmt:skip
 
     ds_atlas = pe.MapNode(
-        Function(
-            function=copy_atlas,
-            input_names=[
-                "name_source",
-                "in_file",
-                "output_dir",
-                "atlas",
-            ],
-            output_names=["out_file"],
-        ),
+        CopyAtlas(output_dir=output_dir),
         name="ds_atlas",
         iterfield=["in_file", "atlas"],
         run_without_submitting=True,
     )
-    ds_atlas.inputs.output_dir = output_dir
     ds_atlas.inputs.atlas = atlases
 
     workflow.connect([
@@ -238,21 +223,11 @@ def init_load_atlases_wf(
     ])  # fmt:skip
 
     ds_atlas_labels_file = pe.MapNode(
-        Function(
-            function=copy_atlas,
-            input_names=[
-                "name_source",
-                "in_file",
-                "output_dir",
-                "atlas",
-            ],
-            output_names=["out_file"],
-        ),
+        CopyAtlas(output_dir=output_dir),
         name="ds_atlas_labels_file",
         iterfield=["in_file", "atlas"],
         run_without_submitting=True,
     )
-    ds_atlas_labels_file.inputs.output_dir = output_dir
     ds_atlas_labels_file.inputs.atlas = atlases
 
     workflow.connect([
@@ -262,21 +237,11 @@ def init_load_atlases_wf(
     ])  # fmt:skip
 
     ds_atlas_metadata = pe.MapNode(
-        Function(
-            function=copy_atlas,
-            input_names=[
-                "name_source",
-                "in_file",
-                "output_dir",
-                "atlas",
-            ],
-            output_names=["out_file"],
-        ),
+        CopyAtlas(output_dir=output_dir),
         name="ds_atlas_metadata",
         iterfield=["in_file", "atlas"],
         run_without_submitting=True,
     )
-    ds_atlas_metadata.inputs.output_dir = output_dir
     ds_atlas_metadata.inputs.atlas = atlases
 
     workflow.connect([


Test patch:
diff --git a/xcp_d/tests/test_interfaces_bids.py b/xcp_d/tests/test_interfaces_bids.py
new file mode 100644
index 0000000..89bf26e
--- /dev/null
+++ b/xcp_d/tests/test_interfaces_bids.py
@@ -0,0 +1,84 @@
+"""Tests for xcp_d.interfaces.bids."""
+
+import os
+
+import pytest
+
+from xcp_d.interfaces import bids
+from xcp_d.utils import atlas
+
+
+def test_copy_atlas(tmp_path_factory):
+    """Test xcp_d.interfaces.bids.CopyAtlas."""
+    tmpdir = tmp_path_factory.mktemp("test_copy_atlas")
+    os.makedirs(os.path.join(tmpdir, "xcp_d"), exist_ok=True)
+
+    # NIfTI
+    atlas_file, _, _ = atlas.get_atlas_nifti("Gordon")
+    name_source = "sub-01_task-A_run-01_space-MNI152NLin2009cAsym_res-2_desc-z_bold.nii.gz"
+    copyatlas = bids.CopyAtlas(
+        name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y"
+    )
+    result = copyatlas.run(cwd=tmpdir)
+    assert os.path.isfile(result.outputs.out_file)
+    assert (
+        os.path.basename(result.outputs.out_file)
+        == "space-MNI152NLin2009cAsym_atlas-Y_res-2_dseg.nii.gz"
+    )
+
+    # Check that the NIfTI file raises an error if the resolution varies
+    # Gordon atlas is 1mm, HCP is 2mm
+    atlas_file_diff_affine, _, _ = atlas.get_atlas_nifti("HCP")
+    with pytest.raises(ValueError, match="is different from the input file affine"):
+        copyatlas = bids.CopyAtlas(
+            name_source=name_source,
+            in_file=atlas_file_diff_affine,
+            output_dir=tmpdir,
+            atlas="Y",
+        )
+        copyatlas.run(cwd=tmpdir)
+
+    # CIFTI
+    atlas_file, atlas_labels_file, atlas_metadata_file = atlas.get_atlas_cifti("Gordon")
+    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
+    copyatlas = bids.CopyAtlas(
+        name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y"
+    )
+    result = copyatlas.run(cwd=tmpdir)
+    assert os.path.isfile(result.outputs.out_file)
+    assert (
+        os.path.basename(result.outputs.out_file) == "space-fsLR_atlas-Y_den-91k_dseg.dlabel.nii"
+    )
+
+    # TSV
+    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
+    copyatlas = bids.CopyAtlas(
+        name_source=name_source, in_file=atlas_labels_file, output_dir=tmpdir, atlas="Y"
+    )
+    result = copyatlas.run(cwd=tmpdir)
+    assert os.path.isfile(result.outputs.out_file)
+    assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.tsv"
+
+    # JSON
+    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
+    copyatlas = bids.CopyAtlas(
+        name_source=name_source, in_file=atlas_metadata_file, output_dir=tmpdir, atlas="Y"
+    )
+    result = copyatlas.run(cwd=tmpdir)
+    assert os.path.isfile(result.outputs.out_file)
+    assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.json"
+
+    # Ensure that out_file isn't overwritten if it already exists
+    fake_in_file = os.path.join(tmpdir, "fake.json")
+    with open(fake_in_file, "w") as fo:
+        fo.write("fake")
+
+    copyatlas = bids.CopyAtlas(
+        name_source=name_source, in_file=fake_in_file, output_dir=tmpdir, atlas="Y"
+    )
+    result = copyatlas.run(cwd=tmpdir)
+    assert os.path.isfile(result.outputs.out_file)
+    assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.json"
+    # The file should not be overwritten, so the contents shouldn't be "fake"
+    with open(result.outputs.out_file, "r") as fo:
+        assert fo.read() != "fake"
diff --git a/xcp_d/tests/test_utils_atlas.py b/xcp_d/tests/test_utils_atlas.py
index 43ffdf7..f877fd7 100644
--- a/xcp_d/tests/test_utils_atlas.py
+++ b/xcp_d/tests/test_utils_atlas.py
@@ -45,57 +45,3 @@ def test_get_atlas_cifti():
 
     with pytest.raises(FileNotFoundError, match="DNE"):
         atlas.get_atlas_cifti("tofail")
-
-
-def test_copy_atlas(tmp_path_factory):
-    """Test xcp_d.utils.atlas.copy_atlas."""
-    tmpdir = tmp_path_factory.mktemp("test_copy_atlas")
-    os.makedirs(os.path.join(tmpdir, "xcp_d"), exist_ok=True)
-
-    # NIfTI
-    atlas_file, _, _ = atlas.get_atlas_nifti("Gordon")
-    name_source = "sub-01_task-A_run-01_space-MNI152NLin2009cAsym_res-2_desc-z_bold.nii.gz"
-    out_file = atlas.copy_atlas(
-        name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y"
-    )
-    assert os.path.isfile(out_file)
-    assert os.path.basename(out_file) == "space-MNI152NLin2009cAsym_atlas-Y_res-2_dseg.nii.gz"
-
-    # CIFTI
-    atlas_file, atlas_labels_file, atlas_metadata_file = atlas.get_atlas_cifti("Gordon")
-    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
-    out_file = atlas.copy_atlas(
-        name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y"
-    )
-    assert os.path.isfile(out_file)
-    assert os.path.basename(out_file) == "space-fsLR_atlas-Y_den-91k_dseg.dlabel.nii"
-
-    # TSV
-    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
-    out_file = atlas.copy_atlas(
-        name_source=name_source, in_file=atlas_labels_file, output_dir=tmpdir, atlas="Y"
-    )
-    assert os.path.isfile(out_file)
-    assert os.path.basename(out_file) == "atlas-Y_dseg.tsv"
-
-    # JSON
-    name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii"
-    out_file = atlas.copy_atlas(
-        name_source=name_source, in_file=atlas_metadata_file, output_dir=tmpdir, atlas="Y"
-    )
-    assert os.path.isfile(out_file)
-    assert os.path.basename(out_file) == "atlas-Y_dseg.json"
-
-    # Ensure that out_file isn't overwritten if it already exists
-    fake_in_file = os.path.join(tmpdir, "fake.json")
-    with open(fake_in_file, "w") as fo:
-        fo.write("fake")
-
-    out_file = atlas.copy_atlas(
-        name_source=name_source, in_file=fake_in_file, output_dir=tmpdir, atlas="Y"
-    )
-    assert os.path.isfile(out_file)
-    assert os.path.basename(out_file) == "atlas-Y_dseg.json"
-    # The file should not be overwritten, so the contents shouldn't be "fake"
-    with open(out_file, "r") as fo:
-        assert fo.read() != "fake"