diff --git a/docs/changes/latest.inc b/docs/changes/latest.inc index 8a91eaeb3..412b4b97b 100644 --- a/docs/changes/latest.inc +++ b/docs/changes/latest.inc @@ -67,6 +67,8 @@ Enhancements - Add more documentation on registering parcellations, coordinates, and masks (:gh:`166` by `Leonard Sasse`_) +- Expose a :func:`junifer.data.parcellations.merge_parcellations` function to merge a list of parcellations (:gh:`146` by Leonard Sasse`_). + Bugs ~~~~ diff --git a/junifer/data/__init__.py b/junifer/data/__init__.py index 4aabf6a9e..56894d7b4 100644 --- a/junifer/data/__init__.py +++ b/junifer/data/__init__.py @@ -13,6 +13,7 @@ from .parcellations import ( list_parcellations, load_parcellation, register_parcellation, + merge_parcellations, ) from .masks import ( diff --git a/junifer/data/parcellations.py b/junifer/data/parcellations.py index b4bf1993e..490cf6408 100644 --- a/junifer/data/parcellations.py +++ b/junifer/data/parcellations.py @@ -16,9 +16,9 @@ import nibabel as nib import numpy as np import pandas as pd import requests -from nilearn import datasets +from nilearn import datasets, image -from ..utils.logging import logger, raise_error +from ..utils.logging import logger, raise_error, warn_with_log from .utils import closest_resolution @@ -700,3 +700,84 @@ def _retrieve_suit( ].to_list() return parcellation_fname, labels + + +def merge_parcellations( + parcellations_list: List["Nifti1Image"], + parcellations_names: List[str], + labels_lists: List[List[str]], +) -> Tuple["Nifti1Image", List[str]]: + """Merge all parcellations from a list into one parcellation. + + Parameters + ---------- + parcellations_list : list of niimg-like object + List of parcellations to merge. + parcellations_names: list of str + List of names for parcellations at the corresponding indices. + labels_lists : list of list of str + A list of lists. Each list in the list contains the labels for the + parcellation at the corresponding index. + + Returns + ------- + parcellation : niimg-like object + The parcellation that results from merging the list of input + parcellations. + labels : list of str + List of labels for the resultant parcellation. + + """ + # Check for duplicated labels + labels_lists_flat = [item for sublist in labels_lists for item in sublist] + if len(labels_lists_flat) != len(set(labels_lists_flat)): + warn_with_log( + "The parcellations have duplicated labels. " + "Each label will be prefixed with the parcellation name." + ) + for i_parcellation, t_labels in enumerate(labels_lists): + labels_lists[i_parcellation] = [ + f"{parcellations_names[i_parcellation]}_{t_label}" + for t_label in t_labels + ] + overlapping_voxels = False + ref_parc = parcellations_list[0] + parc_data = ref_parc.get_fdata() + + labels = labels_lists[0] + + for t_parc, t_labels in zip(parcellations_list[1:], labels_lists[1:]): + if t_parc.shape != ref_parc.shape: + warn_with_log( + "The parcellations have different resolutions!" + "Resampling all parcellations to the first one in the list." + ) + t_parc = image.resample_to_img( + t_parc, ref_parc, interpolation="nearest", copy=True + ) + + # Get the data from this parcellation + t_parc_data = t_parc.get_fdata().copy() # must be copied + # Increase the values of each ROI to match the labels + t_parc_data[t_parc_data != 0] += len(labels) + + # Only set new values for the voxels that are 0 + # This makes sure that the voxels that are in multiple + # parcellations are assigned to the parcellation that was + # first in the list. + if np.any(parc_data[t_parc_data != 0] != 0): + overlapping_voxels = True + + parc_data[parc_data == 0] += t_parc_data[parc_data == 0] + labels.extend(t_labels) + + if overlapping_voxels: + warn_with_log( + "The parcellations have overlapping voxels. " + "The overlapping voxels will be assigned to the " + "parcellation that was first in the list." + ) + + parcellation_img_res = image.new_img_like(parcellations_list[0], parc_data) + + return parcellation_img_res, labels diff --git a/junifer/data/tests/test_parcellations.py b/junifer/data/tests/test_parcellations.py index fcc01c5e2..b68e95fe5 100644 --- a/junifer/data/tests/test_parcellations.py +++ b/junifer/data/tests/test_parcellations.py @@ -9,6 +9,7 @@ from pathlib import Path from typing import List import nibabel as nib +import numpy as np import pytest from nilearn.image import new_img_like from numpy.testing import assert_array_almost_equal, assert_array_equal @@ -20,6 +21,7 @@ from junifer.data.parcellations import ( _retrieve_tian, list_parcellations, load_parcellation, + merge_parcellations, register_parcellation, ) @@ -236,9 +238,7 @@ def test_schaefer_parcellation(tmp_path: Path) -> None: ) # Load parcellation img2, lbl, fname = load_parcellation( - name="Schaefer100x7", - parcellations_dir=tmp_path, - resolution=3, + name="Schaefer100x7", parcellations_dir=tmp_path, resolution=3 ) # Check parcellation values assert fname.name == fname2 @@ -247,9 +247,7 @@ def test_schaefer_parcellation(tmp_path: Path) -> None: assert_array_equal(img2.header["pixdim"][1:4], [2, 2, 2]) # type: ignore # Load parcellation img2, lbl, fname = load_parcellation( - "Schaefer100x7", - parcellations_dir=tmp_path, - resolution=2.1, + "Schaefer100x7", parcellations_dir=tmp_path, resolution=2.1 ) # Check parcellation values assert fname.name == fname2 @@ -258,9 +256,7 @@ def test_schaefer_parcellation(tmp_path: Path) -> None: assert_array_equal(img2.header["pixdim"][1:4], [2, 2, 2]) # type: ignore # Load parcellation img2, lbl, fname = load_parcellation( - "Schaefer100x7", - parcellations_dir=tmp_path, - resolution=1.99, + "Schaefer100x7", parcellations_dir=tmp_path, resolution=1.99 ) # Check parcellation values assert fname.name == fname1 @@ -269,9 +265,7 @@ def test_schaefer_parcellation(tmp_path: Path) -> None: assert_array_equal(img2.header["pixdim"][1:4], [1, 1, 1]) # type: ignore # Load parcellation img2, lbl, fname = load_parcellation( - "Schaefer100x7", - parcellations_dir=tmp_path, - resolution=0.5, + "Schaefer100x7", parcellations_dir=tmp_path, resolution=0.5 ) # Check parcellation values assert fname.name == fname1 @@ -383,18 +377,10 @@ def test_retrieve_suit_incorrect_space(tmp_path: Path) -> None: @pytest.mark.parametrize( - "scale, n_label", - [ - (1, 16), - (2, 32), - (3, 50), - (4, 54), - ], + "scale, n_label", [(1, 16), (2, 32), (3, 50), (4, 54)] ) def test_tian_3T_6thgeneration( - tmp_path: Path, - scale: int, - n_label: int, + tmp_path: Path, scale: int, n_label: int ) -> None: """Test Tian parcellation. @@ -415,8 +401,7 @@ def test_tian_3T_6thgeneration( assert "TianxS4x3TxMNI6thgeneration" in parcellations # Load parcellation img, lbl, fname = load_parcellation( - name=f"TianxS{scale}x3TxMNI6thgeneration", - parcellations_dir=tmp_path, + name=f"TianxS{scale}x3TxMNI6thgeneration", parcellations_dir=tmp_path ) fname1 = f"Tian_Subcortex_S{scale}_3T_1mm.nii.gz" assert img is not None @@ -437,18 +422,10 @@ def test_tian_3T_6thgeneration( @pytest.mark.parametrize( - "scale, n_label", - [ - (1, 16), - (2, 32), - (3, 50), - (4, 54), - ], + "scale, n_label", [(1, 16), (2, 32), (3, 50), (4, 54)] ) def test_tian_3T_nonlinear2009cAsym( - tmp_path: Path, - scale: int, - n_label: int, + tmp_path: Path, scale: int, n_label: int ) -> None: """Test Tian parcellation. @@ -480,18 +457,10 @@ def test_tian_3T_nonlinear2009cAsym( @pytest.mark.parametrize( - "scale, n_label", - [ - (1, 16), - (2, 34), - (3, 54), - (4, 62), - ], + "scale, n_label", [(1, 16), (2, 34), (3, 54), (4, 62)] ) def test_tian_7T_6thgeneration( - tmp_path: Path, - scale: int, - n_label: int, + tmp_path: Path, scale: int, n_label: int ) -> None: """Test Tian parcellation. @@ -534,10 +503,7 @@ def test_retrieve_tian_incorrect_space(tmp_path: Path) -> None: """ with pytest.raises(ValueError, match=r"The parameter `space`"): _retrieve_tian( - parcellations_dir=tmp_path, - resolution=1, - scale=1, - space="wrong", + parcellations_dir=tmp_path, resolution=1, scale=1, space="wrong" ) with pytest.raises(ValueError, match=r"MNI6thgeneration"): @@ -584,3 +550,162 @@ def test_retrieve_tian_incorrect_scale(tmp_path: Path) -> None: scale=5, space="MNI6thgeneration", ) + + +def test_merge_parcellations() -> None: + """Test merging parcellations.""" + # load some parcellations for testing + schaefer_parcellation, schaefer_labels, _ = load_parcellation( + "Schaefer100x17" + ) + tian_parcellation, tian_labels, _ = load_parcellation( + "TianxS2x3TxMNInonlinear2009cAsym" + ) + # prepare the list of the actual parcellations + parcellation_list = [schaefer_parcellation, tian_parcellation] + # prepare a list of names + names = ["Schaefer100x17", "TianxS2x3TxMNInonlinear2009cAsym"] + # prepare a list of label lists + labels_lists = [schaefer_labels, tian_labels] + # merge the parcellations + merged_parc, labels = merge_parcellations( + parcellation_list, names, labels_lists + ) + + # we should have 132 integer labels plus 1 for background + parc_data = merged_parc.get_fdata() + assert len(np.unique(parc_data)) == 133 + # no background label, so labels is one less + assert len(labels) == 132 + + +def test_merge_parcellations_3D_multiple_non_overlapping( + tmp_path: Path, +) -> None: + """Test merge_parcellations with multiple non-overlapping parcellations. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + # Get the testing parcellation + parcellation, labels, _ = load_parcellation("Schaefer100x7") + + assert parcellation is not None + + # Create two parcellations from it + parcellation_data = parcellation.get_fdata() + parcellation1_data = parcellation_data.copy() + parcellation1_data[parcellation1_data > 50] = 0 + parcellation2_data = parcellation_data.copy() + parcellation2_data[parcellation2_data <= 50] = 0 + parcellation2_data[parcellation2_data > 0] -= 50 + labels1 = labels[:50] + labels2 = labels[50:] + + parcellation1_img = new_img_like(parcellation, parcellation1_data) + parcellation2_img = new_img_like(parcellation, parcellation2_data) + + parcellation_list = [parcellation1_img, parcellation2_img] + names = ["high", "low"] + labels_lists = [labels1, labels2] + + merged_parc, merged_labels = merge_parcellations( + parcellation_list, names, labels_lists + ) + + parc_data = parcellation.get_fdata() + assert_array_equal(parc_data, merged_parc.get_fdata()) + assert len(labels) == 100 + assert len(np.unique(parc_data)) == 101 # 100 + 1 because background 0 + + +def test_merge_parcellations_3D_multiple_overlapping(tmp_path: Path) -> None: + """Test merge_parcellations with multiple overlapping parcellations. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + + # Get the testing parcellation + parcellation, labels, _ = load_parcellation("Schaefer100x7") + + assert parcellation is not None + + # Create two parcellations from it + parcellation_data = parcellation.get_fdata() + parcellation1_data = parcellation_data.copy() + parcellation1_data[parcellation1_data > 50] = 0 + parcellation2_data = parcellation_data.copy() + + # Make the second parcellation overlap with the first + parcellation2_data[parcellation2_data <= 45] = 0 + parcellation2_data[parcellation2_data > 0] -= 45 + labels1 = [f"low_{x}" for x in labels[:50]] # Change the labels + labels2 = [f"high_{x}" for x in labels[45:]] # Change the labels + + parcellation1_img = new_img_like(parcellation, parcellation1_data) + parcellation2_img = new_img_like(parcellation, parcellation2_data) + + parcellation_list = [parcellation1_img, parcellation2_img] + names = ["high", "low"] + labels_lists = [labels1, labels2] + + with pytest.warns(RuntimeWarning, match="overlapping voxels"): + merged_parc, merged_labels = merge_parcellations( + parcellation_list, names, labels_lists + ) + + parc_data = parcellation.get_fdata() + assert len(labels) == 100 + assert len(np.unique(parc_data)) == 101 # 100 + 1 because background 0 + + +def test_merge_parcellations_3D_multiple_duplicated_labels( + tmp_path: Path, +) -> None: + """Test merge_parcellations with two parcellations with duplicated labels. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + + # Get the testing parcellation + parcellation, labels, _ = load_parcellation("Schaefer100x7") + + assert parcellation is not None + + # Create two parcellations from it + parcellation_data = parcellation.get_fdata() + parcellation1_data = parcellation_data.copy() + parcellation1_data[parcellation1_data > 50] = 0 + parcellation2_data = parcellation_data.copy() + parcellation2_data[parcellation2_data <= 50] = 0 + parcellation2_data[parcellation2_data > 0] -= 50 + labels1 = labels[:50] + labels2 = labels[49:-1] # One label is duplicated + + parcellation1_img = new_img_like(parcellation, parcellation1_data) + parcellation2_img = new_img_like(parcellation, parcellation2_data) + + parcellation_list = [parcellation1_img, parcellation2_img] + names = ["high", "low"] + labels_lists = [labels1, labels2] + + with pytest.warns(RuntimeWarning, match="duplicated labels."): + merged_parc, merged_labels = merge_parcellations( + parcellation_list, names, labels_lists + ) + + parc_data = parcellation.get_fdata() + assert_array_equal(parc_data, merged_parc.get_fdata()) + assert len(labels) == 100 + assert len(np.unique(parc_data)) == 101 # 100 + 1 because background 0 diff --git a/junifer/datagrabber/datalad_base.py b/junifer/datagrabber/datalad_base.py index d6ad0e16f..cd732de75 100644 --- a/junifer/datagrabber/datalad_base.py +++ b/junifer/datagrabber/datalad_base.py @@ -178,9 +178,7 @@ class DataladDataGrabber(BaseDataGrabber): try: dl_out = self._dataset.get(to_get, result_renderer="disabled") except IncompleteResultsError as e: - raise_error( - f"Failed to get from dataset: {e.failed}" - ) + raise_error(f"Failed to get from dataset: {e.failed}") if not self._was_cloned: # If the dataset was already installed, check that the # file was actually downloaded to avoid removing a diff --git a/junifer/markers/parcel_aggregation.py b/junifer/markers/parcel_aggregation.py index 075464125..9cb5b69d6 100644 --- a/junifer/markers/parcel_aggregation.py +++ b/junifer/markers/parcel_aggregation.py @@ -7,13 +7,13 @@ from typing import Any, Dict, List, Optional, Union import numpy as np -from nilearn.image import math_img, new_img_like, resample_to_img +from nilearn.image import math_img, resample_to_img from nilearn.maskers import NiftiMasker from ..api.decorators import register_marker -from ..data import get_mask, load_parcellation +from ..data import get_mask, load_parcellation, merge_parcellations from ..stats import get_aggfunc_by_name -from ..utils import logger, warn_with_log +from ..utils import logger from .base import BaseMarker @@ -98,9 +98,7 @@ class ParcelAggregation(BaseMarker): raise ValueError(f"Unknown input kind for {input_type}") def compute( - self, - input: Dict[str, Any], - extra_input: Optional[Dict] = None, + self, input: Dict[str, Any], extra_input: Optional[Dict] = None ) -> Dict: """Compute. @@ -129,8 +127,7 @@ class ParcelAggregation(BaseMarker): t_input_img = input["data"] logger.debug(f"Parcel aggregation using {self.method}") agg_func = get_aggfunc_by_name( - name=self.method, - func_params=self.method_params, + name=self.method, func_params=self.method_params ) # Get the min of the voxels sizes and use it as the resolution resolution = np.min(t_input_img.header.get_zooms()[:3]) @@ -140,15 +137,11 @@ class ParcelAggregation(BaseMarker): all_labels = [] for t_parc_name in self.parcellation: t_parcellation, t_labels, _ = load_parcellation( - name=t_parc_name, - resolution=resolution, + name=t_parc_name, resolution=resolution ) # Resample all of them to the image t_parcellation_img_res = resample_to_img( - t_parcellation, - t_input_img, - interpolation="nearest", - copy=True, + t_parcellation, t_input_img, interpolation="nearest", copy=True ) all_parcelations.append(t_parcellation_img_res) all_labels.append(t_labels) @@ -159,56 +152,11 @@ class ParcelAggregation(BaseMarker): labels = all_labels[0] else: # Merge the parcellations - - # Check for duplicated labels - all_labels_flat = [ - item for sublist in all_labels for item in sublist - ] - if len(all_labels_flat) != len(set(all_labels_flat)): - warn_with_log( - "The parcellations have duplicated labels. " - "Each label will be prefixed with the parcellation name." - ) - for i_parcellation, t_labels in enumerate(all_labels): - all_labels[i_parcellation] = [ - f"{self.parcellation[i_parcellation]}_{t_label}" - for t_label in t_labels - ] - overlapping_voxels = False - parc_data = all_parcelations[0].get_fdata() - labels = all_labels[0] - for t_parc, t_labels in zip(all_parcelations[1:], all_labels[1:]): - # Get the data from this parcellation - t_parc_data = t_parc.get_fdata().copy() # must be copied - # Increase the values of each ROI to match the labels - t_parc_data[t_parc_data != 0] += len(labels) - - # Only set new values for the voxels that are 0 - # This makes sure that the voxels that are in multiple - # parcellations are assigned to the parcellation that was - # first in the list. - if np.any(parc_data[t_parc_data != 0] != 0): - overlapping_voxels = True - - parc_data[parc_data == 0] += t_parc_data[parc_data == 0] - labels.extend(t_labels) - - if overlapping_voxels: - warn_with_log( - "The parcellations have overlapping voxels. " - "The overlapping voxels will be assigned to the " - "parcellation that was first in the list." - ) - - parcellation_img_res = new_img_like( - all_parcelations[0], - parc_data, + parcellation_img_res, labels = merge_parcellations( + all_parcelations, self.parcellation, all_labels ) - parcellation_bin = math_img( - "img != 0", - img=parcellation_img_res, - ) + parcellation_bin = math_img("img != 0", img=parcellation_img_res) if self.masks is not None: logger.debug(f"Masking with {self.masks}")