add knn interpolation parameter to help with projecting from volume t…

…o surface space

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

neurocaps/__init__.py

@@ -3,4 +3,4 @@
 __all__=["analysis", "extraction"]
 
 # Version in single place
-__version__ = "0.13.1"
+__version__ = "0.13.2"

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
+    # 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

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()

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)