From dc4a9f2111d2d0b98f6c02342ed5664c1cc33b2a Mon Sep 17 00:00:00 2001 From: Donisha Smith <112973674+donishadsmith@users.noreply.github.com> Date: Fri, 5 Jul 2024 00:28:14 -0400 Subject: [PATCH] add knn interpolation parameter to help with projecting from volume to surface space --- CHANGELOG.md | 7 ++ neurocaps/__init__.py | 2 +- .../_utils/_cap_internals/_cap2statmap.py | 71 +++++++++++++++-- neurocaps/analysis/cap.py | 77 ++++++++++++------- neurocaps/extraction/timeseriesextractor.py | 2 +- 5 files changed, 124 insertions(+), 35 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb17b18..a1cfa15 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,6 +41,13 @@ noted in the changelog (i.e new functions or parameters, changes in parameter de - *.patch* : Contains no new features, simply fixes any identified bugs. - *.postN* : Consists of only metadata-related changes, such as updates to type hints or doc strings/documentation. +## [0.13.2] - 2024-07-05 +### 🐛 Fixes +- Certain custom atlases may not project well from volume to surface space. A new parameter, `knn_dict` has been added to +`CAP.caps2surf()` and `CAP.caps2niftis()` to apply k-nearest neighbors (knn) interpolation while leveraging the +Schaefer atlas, which projects well from volumetric to surface space. +- No longer need to add `parcel_approach` when using `CAP.caps2surf()` with `fslr_giftis_dict`. + ## [0.13.1] - 2024-06-30 ### ♻ Changed - For `CAP.caps2radar()`, the `scattersize` kwarg can be used to control the size of the scatter/markers regardless diff --git a/neurocaps/__init__.py b/neurocaps/__init__.py index 14b5290..6a3333b 100644 --- a/neurocaps/__init__.py +++ b/neurocaps/__init__.py @@ -3,4 +3,4 @@ __all__=["analysis", "extraction"] # Version in single place -__version__ = "0.13.1" \ No newline at end of file +__version__ = "0.13.2" \ No newline at end of file diff --git a/neurocaps/_utils/_cap_internals/_cap2statmap.py b/neurocaps/_utils/_cap_internals/_cap2statmap.py index e5cb2d7..4f512d6 100644 --- a/neurocaps/_utils/_cap_internals/_cap2statmap.py +++ b/neurocaps/_utils/_cap_internals/_cap2statmap.py @@ -1,18 +1,75 @@ """Internal function for turning CAPs into NifTI Statistical Maps""" +import textwrap, warnings import nibabel as nib, numpy as np -from nilearn import image +from nilearn import datasets, image +from scipy.spatial import cKDTree -def _cap2statmap(atlas_file, cap_vector, fwhm): +def _cap2statmap(atlas_file, cap_vector, fwhm, knn_dict): atlas = nib.load(atlas_file) atlas_fdata = atlas.get_fdata() - # Get array containing all labels in atlas to avoid issue if atlas labels dont start at 1, like Nilearn's AAL map + # Get array containing all labels in atlas to avoid issue if the first non-zero atlas label is not 1 target_array = sorted(np.unique(atlas_fdata)) - for indx, value in enumerate(cap_vector): - actual_indx = indx + 1 - atlas_fdata[np.where(atlas_fdata == target_array[actual_indx])] = value + for indx, value in enumerate(cap_vector, start=1): + atlas_fdata[np.where(atlas_fdata == target_array[indx])] = value stat_map = nib.Nifti1Image(atlas_fdata, atlas.affine, atlas.header) + # Add smoothing to stat map to help mitigate potential coverage issues if fwhm is not None: stat_map = image.smooth_img(stat_map, fwhm=fwhm) - return stat_map \ No newline at end of file + # Knn implementation to aid in coverage issues + if knn_dict: + # Get target indices + target_indices = _get_target_indices(atlas=atlas, knn_dict=knn_dict) + # Get non-zero indices of the stat map + non_zero_indices = np.array(np.where(stat_map.get_fdata() != 0)).T + if "k" not in knn_dict: + warnings.warn("Defaulting to k=1 since 'k' was not specified in `knn_dict`.") + k = 1 + else: k = knn_dict["k"] + + # Build kdtree for nearest neighbors + kdtree = cKDTree(non_zero_indices) + + for target_indx in target_indices: + # Get the nearest non-zero index + _ , neighbor_indx = kdtree.query(target_indx, k = k) + nearest_neighbors = non_zero_indices[neighbor_indx] + + if k > 1: + # Values of nearest neighbors + neighbor_values = [stat_map.get_fdata()[tuple(nearest_neighbor)] for nearest_neighbor in nearest_neighbors] + # Majority vote + new_value = np.bincount(neighbor_values).argmax() + else: + nearest_neighbor = non_zero_indices[neighbor_indx] + new_value = stat_map.get_fdata()[tuple(nearest_neighbor)] + + # Assign the new value to the current index + stat_map.get_fdata()[tuple(target_indx)] = new_value + + return stat_map + +def _get_target_indices(atlas, knn_dict): + # Get schaefer atlas, which projects well onto cortical surface plots + if "resolution_mm" not in knn_dict: + warnings.warn(textwrap.dedent(""" + Defaulting to 1mm resolution for the Schaefer atlas since 'resolution_knn' was + not specified in `knn_dict`. + """)) + resolution_mm = "1mm" + else: + resolution_mm = knn_dict["resolution_mm"] + + schaefer_atlas = datasets.fetch_atlas_schaefer_2018(resolution_mm=resolution_mm)["maps"] + + # Resample schaefer to atlas file + resampled_schaefer = image.resample_to_img(schaefer_atlas, atlas) + # Get indices that equal zero in schaefer atlas to avoid interpolating background values, will also get the indices for subcortical + background_indices_schaefer = set(zip(*np.where(resampled_schaefer.get_fdata() == 0))) + # Get indices 0 indices for atlas + background_indices_atlas = set(zip(*np.where(atlas.get_fdata() == 0))) + # Get the non-background indices through subtraction + target_indices = list(background_indices_atlas - background_indices_schaefer) + + return target_indices diff --git a/neurocaps/analysis/cap.py b/neurocaps/analysis/cap.py index c2a0dff..23d19ae 100644 --- a/neurocaps/analysis/cap.py +++ b/neurocaps/analysis/cap.py @@ -1775,7 +1775,7 @@ def caps2corr(self, output_dir: Optional[os.PathLike]=None, suffix_title: Option if return_df: return corr_dict def caps2niftis(self, output_dir: os.PathLike, suffix_file_name: Optional[str]=None, - fwhm: Optional[float]=None) -> nib.Nifti1Image: + fwhm: Optional[float]=None, knn_dict: dict[str, Union[int, bool]]=None) -> nib.Nifti1Image: """ **Standalone Method to Convert CAPs to NifTI Statistical Maps** @@ -1787,17 +1787,11 @@ def caps2niftis(self, output_dir: os.PathLike, suffix_file_name: Optional[str]=N def _cap2statmap(atlas_file, cap_vector, fwhm): atlas = nib.load(atlas_file) atlas_fdata = atlas.get_fdata() - # Get array containing all labels in atlas to avoid issue if atlas labels dont start at 1, - # like Nilearn's AAL map + # Get array containing all labels in atlas to avoid issue if the first non-zero atlas label is not 1 target_array = sorted(np.unique(atlas_fdata)) - for indx, value in enumerate(cap_vector): - actual_indx = indx + 1 - atlas_fdata[np.where(atlas_fdata == target_array[actual_indx])] = value + for indx, value in enumerate(cap_vector, start=1): + atlas_fdata[np.where(atlas_fdata == target_array[indx])] = value stat_map = nib.Nifti1Image(atlas_fdata, atlas.affine, atlas.header) - # Add smoothing to stat map to help mitigate potential coverage issues - if fwhm is not None: - stat_map = image.smooth_img(stat_map, fwhm=fwhm) - retutn stat_map Parameters ---------- @@ -1812,6 +1806,24 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): from MNI152 space to fslr surface space. Note, this can assist with coverage issues in the plot. Uses ``nilearn.image.smooth_img``. + knn_dict : :obj:`dict[str, int | bool]`, default=None + Use KNN (k-nearest neighbors) interpolation to fill in non-background values that are assigned zero. This is + primarily used as a fix for when a custom parcellation does not project well from volumetric to surface + space. This method involves resampling the Schaefer parcellation, a volumetric parcellation that projects + well onto surface space, to the target parcellation specified in the "maps" sub-key in ``self.parcel_approach``. + The background indices are extracted from the Schaefer parcellation, and these indices are used to obtain + the non-background indices (parcels) that are set to zero in the target parcellation. + + These indices are then iterated over, and the zero values are replaced with the value of the nearest neighbor, + determined by the sub-key "k". The dictionary contains the following sub-keys: + + - "k" : An integer that determines the number of nearest neighbors to consider, with the majority vote determining the new value. If not specified, the default is 1. + - "resolution_mm" : An integer (1 or 2) that determines the resolution of the Schaefer parcellation. If not specified, the default is 1. + + This method is applied after the `fwhm` method. + + .. versionadded:: 0.13.2 + Returns ------- `NifTI1Image` @@ -1837,7 +1849,7 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): for group in self._caps: for cap in self._caps[group]: stat_map = _cap2statmap(atlas_file=self._parcel_approach[parcellation_name]["maps"], - cap_vector=self._caps[group][cap], fwhm=fwhm) + cap_vector=self._caps[group][cap], fwhm=fwhm, knn_dict=knn_dict) if suffix_file_name: save_name = f"{group.replace(' ', '_')}_{cap.replace('-', '_')}_{suffix_file_name}.nii.gz" @@ -1848,7 +1860,8 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): def caps2surf(self, output_dir: Optional[os.PathLike]=None, suffix_title: Optional[str]=None, show_figs: bool=True, fwhm: Optional[float]=None, fslr_density: Literal["4k", "8k", "32k", "164k"]="32k", method: Literal["linear", "nearest"]="linear", - save_stat_map: bool=False, fslr_giftis_dict: Optional[dict]=None, **kwargs) -> surfplot.Plot: + save_stat_map: bool=False, fslr_giftis_dict: Optional[dict]=None, + knn_dict: dict[str, Union[int, bool]]=None, **kwargs) -> surfplot.Plot: """ **Project CAPs onto Surface Plots** @@ -1866,17 +1879,11 @@ def caps2surf(self, output_dir: Optional[os.PathLike]=None, suffix_title: Option def _cap2statmap(atlas_file, cap_vector, fwhm): atlas = nib.load(atlas_file) atlas_fdata = atlas.get_fdata() - # Get array containing all labels in atlas to avoid issue if atlas labels dont start at 1, - # like Nilearn's AAL map + # Get array containing all labels in atlas to avoid issue if the first non-zero atlas label is not 1 target_array = sorted(np.unique(atlas_fdata)) - for indx, value in enumerate(cap_vector): - actual_indx = indx + 1 - atlas_fdata[np.where(atlas_fdata == target_array[actual_indx])] = value + for indx, value in enumerate(cap_vector, start=1): + atlas_fdata[np.where(atlas_fdata == target_array[indx])] = value stat_map = nib.Nifti1Image(atlas_fdata, atlas.affine, atlas.header) - # Add smoothing to stat map to help mitigate potential coverage issues - if fwhm is not None: - stat_map = image.smooth_img(stat_map, fwhm=fwhm) - return stat_map Parameters ---------- @@ -1926,6 +1933,24 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): parameter allows plotting without re-running the analysis. Initialize the CAP class and use this method if using this parameter. + knn_dict : :obj:`dict[str, int | bool]`, default=None + Use KNN (k-nearest neighbors) interpolation to fill in non-background values that are assigned zero. This is + primarily used as a fix for when a custom parcellation does not project well from volumetric to surface + space. This method involves resampling the Schaefer parcellation, a volumetric parcellation that projects + well onto surface space, to the target parcellation specified in the "maps" sub-key in ``self.parcel_approach``. + The background indices are extracted from the Schaefer parcellation, and these indices are used to obtain + the non-background indices (parcels) that are set to zero in the target parcellation. + + These indices are then iterated over, and the zero values are replaced with the value of the nearest neighbor, + determined by the sub-key "k". The dictionary contains the following sub-keys: + + - "k": An integer that determines the number of nearest neighbors to consider, with the majority vote determining the new value. If not specified, the default is 1. + - "resolution_mm": An integer (1 or 2) that determines the resolution of the Schaefer parcellation. If not specified, the default is 1. + + This method is applied after the `fwhm` method. + + .. versionadded:: 0.13.2 + kwargs : :obj:`dict` Additional parameters to pass to modify certain plot parameters. Options include: @@ -1985,7 +2010,7 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): centroid vector is the first nonzero label, which is assumed to be at the first index of the array in ``sorted(np.unique(atlas_fdata))``. """ - if not self._parcel_approach: + if not self._parcel_approach and fslr_giftis_dict is None: raise AttributeError(textwrap.dedent(""" `self.parcel_approach` is None. Add parcel_approach using `self.parcel_approach=parcel_approach` to use this @@ -2010,14 +2035,14 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): groups = self._caps if hasattr(self,"_caps") and fslr_giftis_dict is None else fslr_giftis_dict - parcellation_name = list(self._parcel_approach)[0] + if fslr_giftis_dict is None: parcellation_name = list(self._parcel_approach)[0] for group in groups: caps = self._caps[group] if hasattr(self,"_caps") and fslr_giftis_dict is None else fslr_giftis_dict[group] for cap in caps: if fslr_giftis_dict is None: stat_map = _cap2statmap(atlas_file=self._parcel_approach[parcellation_name]["maps"], - cap_vector=self._caps[group][cap], fwhm=fwhm) + cap_vector=self._caps[group][cap], fwhm=fwhm, knn_dict=knn_dict) # Fix for python 3.12, saving stat map so that it is path instead of a NifTi try: @@ -2040,8 +2065,8 @@ def _cap2statmap(atlas_file, cap_vector, fwhm): # Delete os.unlink(temp_nifti.name) else: - gii_lh, gii_rh = fslr_to_fslr([fslr_giftis_dict[group][cap]["lh"], - fslr_giftis_dict[group][cap]["rh"]], + gii_lh, gii_rh = fslr_to_fslr((fslr_giftis_dict[group][cap]["lh"], + fslr_giftis_dict[group][cap]["rh"]), target_density=fslr_density, method=method) # Code slightly adapted from surfplot example 2 surfaces = fetch_fslr() diff --git a/neurocaps/extraction/timeseriesextractor.py b/neurocaps/extraction/timeseriesextractor.py index 774ec46..5f40cef 100644 --- a/neurocaps/extraction/timeseriesextractor.py +++ b/neurocaps/extraction/timeseriesextractor.py @@ -335,7 +335,7 @@ def get_bold(self, bids_dir: os.PathLike, task: str, session: Optional[Union[int **This pipeline is most optimized towards BOLD data preprocessed by fMRIPrep.** When extracting specific conditions, this function uses ``math.ceil`` when calculating the duration of a - condition to round up. + condition to round up and ``int`` to round down for the onset. :: onset_scan = int(condition_df.loc[i,"onset"]/tr)