TimeSeriesVolume and assignment

Author: AhmetNSimsek

e2e/volumes/test_time_series_volume.py

@@ -0,0 +1,35 @@
+import numpy as np
+import nibabel
+import siibra
+
+
+def create_synthetic_data():
+    julich_pmaps = siibra.get_map(
+        parcellation="julich 2.9",
+        space="mni152",
+        maptype="statistical"
+    )
+    length = 20
+    template_img = julich_pmaps.space.get_template().fetch()
+    arr = np.zeros(list(template_img.shape) + [length])
+    for i, img in enumerate(julich_pmaps.fetch_iter()):
+        if i == length:
+            break
+        arr[:, :, :, i] += img.dataobj
+    return siibra.volumes.from_nifti(
+        nibabel.nifti1.Nifti1Image(arr, affine=template_img.affine),
+        time_index=np.asanyarray(range(length)),
+        space='mni152',
+        name="synthetic timeseries volume"
+    )
+
+
+def test_timeseries_volume_assignment():
+    difumo128 = siibra.get_map(
+        parcellation="difumo 128",
+        space="mni152",
+        maptype="statistical"
+    )
+    synthetic_vol = create_synthetic_data()
+    assignments = difumo128.assign(synthetic_vol, split_components=False)
+    print(assignments)

siibra/volumes/parcellationmap.py

@@ -56,6 +56,7 @@ class MapAssignment:
     volume: int
     fragment: str
     map_value: np.ndarray
+    time_index: Union[int, None]
 
 
 @dataclass
@@ -893,12 +894,20 @@ def _assign(
         if isinstance(item, pointcloud.PointCloud):
             return self._assign_points(item, lower_threshold)
         if isinstance(item, _volume.Volume):
-            return self._assign_volume(
-                queryvolume=item,
-                lower_threshold=lower_threshold,
-                minsize_voxel=minsize_voxel,
-                **kwargs
-            )
+            if isinstance(item, _volume.TimeSeriesVolume):
+                return self._assign_timeseries_volume(
+                    queryvolume=item,
+                    lower_threshold=lower_threshold,
+                    minsize_voxel=minsize_voxel,
+                    **kwargs
+                )
+            else:
+                return self._assign_volume(
+                    queryvolume=item,
+                    lower_threshold=lower_threshold,
+                    minsize_voxel=minsize_voxel,
+                    **kwargs
+                )
 
         raise RuntimeError(
             f"Items of type {item.__class__.__name__} cannot be used for region assignment."
@@ -957,6 +966,7 @@ def assign(
         # format assignments as pandas dataframe
         columns = [
             "input structure",
+            "time_index",
             "centroid",
             "volume",
             "fragment",
@@ -970,7 +980,7 @@ def assign(
             "input containedness"
         ]
         if len(assignments) == 0:
-            return pd.DataFrame(columns=columns)
+            return pd.DataFrame(columns=columns).dropna(axis='columns', how='all')
         # determine the unique set of observed indices in order to do region lookups
         # only once for each map index occurring in the point list
         labelled = self.is_labelled  # avoid calling this in a loop
@@ -1147,7 +1157,8 @@ def _assign_volume(
         queryvolume: "_volume.Volume",
         lower_threshold: float,
         split_components: bool = True,
-        **kwargs
+        time_index: int = None,
+        **kwargs,
     ) -> List[AssignImageResult]:
         """
         Assign an image volume to this parcellation map.
@@ -1214,12 +1225,32 @@ def _assign_volume(
                                 volume=index.volume,
                                 fragment=index.fragment,
                                 map_value=index.label,
+                                time_index=time_index,
                                 **asdict(scores)
                             )
                         )
 
         return assignments
 
+    def _assign_timeseries_volume(
+        self,
+        queryvolume: "_volume.TimeSeriesVolume",
+        lower_threshold: float,
+        split_components: bool = True,
+        **kwargs
+    ) -> List[AssignImageResult]:
+        assignments = []
+        for v_t in siibra_tqdm(queryvolume, unit='time_index'):
+            assignments_t = self._assign_volume(
+                v_t,
+                lower_threshold=lower_threshold,
+                split_components=split_components,
+                time_index=v_t.time_index,
+                **kwargs
+            )
+            assignments.extend(assignments_t)
+        return assignments
+
 
 def from_volume(
     name: str,

siibra/volumes/sparsemap.py

@@ -335,7 +335,8 @@ def _assign_volume(
         queryvolume: "_volume.Volume",
         minsize_voxel: int,
         lower_threshold: float,
-        split_components: bool = True
+        split_components: bool = True,
+        time_index: int = None,
     ) -> List[parcellationmap.AssignImageResult]:
         """
         Assign an image volume to this sparse map.
@@ -444,6 +445,7 @@ def _assign_volume(
 
                 assignments.append(
                     parcellationmap.AssignImageResult(
+                        time_index=time_index,
                         input_structure=mode,
                         centroid=tuple(position.round(2)),
                         volume=volume,

siibra/volumes/volume.py

@@ -14,7 +14,7 @@
 # limitations under the License.
 """A specific mesh or 3D array."""
 
-from typing import List, Dict, Union, Set, TYPE_CHECKING
+from typing import Iterable, List, Dict, Union, Set, TYPE_CHECKING
 from dataclasses import dataclass
 from time import sleep
 import json
@@ -633,6 +633,7 @@ def __init__(
         label: int = None,
         fragment: str = None,
         threshold: float = None,
+        time_index: int = None,
     ):
         """
         A prescribed Volume to fetch specified label and fragment.
@@ -647,6 +648,8 @@ def __init__(
             If a volume is fragmented, get a specified one.
         threshold : float, default None
             Provide a float value to threshold the image.
+        time_index: int = None,
+            If parent volume is a timeseries Nifti, filter a time index without fetching the full image.
         """
         name = parent_volume.name
         if label:
@@ -655,6 +658,8 @@ def __init__(
             name += f" - fragment: {fragment}"
         if threshold:
             name += f" - threshold: {threshold}"
+        if time_index:
+            name += f" - time index: {time_index}"
         Volume.__init__(
             self,
             space_spec=parent_volume._space_spec,
@@ -664,6 +669,7 @@ def __init__(
         self.fragment = fragment
         self.label = label
         self.threshold = threshold
+        self.time_index = time_index
 
     def fetch(
         self,
@@ -680,7 +686,8 @@ def fetch(
             kwargs["label"] = self.label
 
         result = super().fetch(format=format, **kwargs)
-
+        if self.time_index is not None:
+            result = result.slicer[:, :, :, self.time_index]
         if self.threshold is not None:
             assert self.label is None
             if not isinstance(result, Nifti1Image):
@@ -710,6 +717,29 @@ def get_boundingbox(
         )
 
 
+class TimeSeriesVolume(Volume):
+    def __init__(
+        self,
+        time_index: np.ndarray,
+        **kwargs,
+    ):
+        Volume.__init__(self, **kwargs)
+        self.time_index = time_index
+
+    def __iter__(self) -> Iterable[FilteredVolume]:
+        yield from (
+            FilteredVolume(parent_volume=self, time_index=t)
+            for t in self.time_index
+        )
+
+    def get_index(self, time_index: int):
+        return FilteredVolume(parent_volume=self, time_index=time_index)
+
+    def fetch(self, format: str = None, time_index: int = None, **kwargs):
+        img = super().fetch(format, **kwargs)
+        return img.slicer[:, :, :, time_index] if time_index else img
+
+
 class ReducedVolume(Volume):
     def __init__(
         self,
@@ -832,14 +862,17 @@ def from_file(filename: str, space: str, name: str) -> Volume:
     )
 
 
-def from_nifti(nifti: Nifti1Image, space: str, name: str) -> Volume:
+def from_nifti(nifti: Nifti1Image, space: str, name: str, time_index: np.ndarray = None) -> Union[Volume, TimeSeriesVolume]:
     """Builds a nifti volume from a Nifti image."""
     spaceobj = get_registry("Space").get(space)
-    return Volume(
+    kwargs = dict(
         space_spec={"@id": spaceobj.id},
         providers=[_providers.NiftiProvider((np.asanyarray(nifti.dataobj), nifti.affine))],
         name=name
     )
+    if time_index is None:
+        return Volume(**kwargs)
+    return TimeSeriesVolume(time_index=time_index, **kwargs)
 
 
 def from_array(