Implementation of Geometric Subsegmentation into Feature Detection

environment-examples.yml

@@ -37,3 +37,4 @@ dependencies:
   - trame-vuetify
   - ipywidgets
   - vtk>=9.4.2
+  - SimpleITK

example_requirements.txt

@@ -29,3 +29,5 @@ ipywidgets
 trame
 trame-vtk
 trame-vuetify
+SimpleITK
+mahotas

pyproject.toml

@@ -63,6 +63,12 @@ dev = [
     "nbconvert",
 ]
 
+# optional dependencies for advanced segmentation features
+segmentation = [
+    "SimpleITK",
+    "mahotas",
+]
+
 # required to run example notebooks
 examples = [
     "jupyter",
@@ -79,6 +85,7 @@ examples = [
     "intake-xarray==0.7.0",
     "healpix",
     "easygems",
+    "SimpleITK",
 ]
 
 # All extras combined
@@ -87,6 +94,8 @@ all = [
     "black",
     "pytest",
     "nbconvert",
+    "SimpleITK",
+    "mahotas",
     "jupyter",
     "notebook",
     "pytables",

tobac/feature_detection.py

@@ -36,6 +36,7 @@
 from tobac.utils import periodic_boundaries as pbc_utils
 from tobac.utils.general import spectral_filtering
 from tobac.utils.generators import field_and_features_over_time
+from tobac.segmentation import subsegmentation
 
 
 def feature_position(
@@ -478,6 +479,10 @@ def feature_detection_threshold(
     from skimage.morphology import binary_erosion
     from copy import deepcopy
 
+    feature_detection_method = kwargs.get(
+        "feature_detection_method", "connected_components"
+    )
+
     if min_num != 0:
         warnings.warn(
             "min_num parameter has no effect and will be deprecated in a future version of tobac. "
@@ -501,18 +506,43 @@ def feature_detection_threshold(
         # if looking for minima, set values above threshold to 0 and scale by data minimum:
     elif target == "minimum":
         mask = data_i <= threshold
-    # only include values greater than threshold
-    # erode selected regions by n pixels
-    if n_erosion_threshold > 0:
-        if is_3D:
-            selem = np.ones(
-                (n_erosion_threshold, n_erosion_threshold, n_erosion_threshold)
-            )
+
+    # selecting feature detection implementations
+    # ===========================================
+    if feature_detection_method == "connected_components":
+        print("... using connected components for feature detection")
+
+        # only include values greater than threshold
+        # erode selected regions by n pixels
+        if n_erosion_threshold > 0:
+            if is_3D:
+                selem = np.ones(
+                    (n_erosion_threshold, n_erosion_threshold, n_erosion_threshold)
+                )
+            else:
+                selem = np.ones((n_erosion_threshold, n_erosion_threshold))
+            mask = binary_erosion(mask, selem)
+            # detect individual regions, label  and count the number of pixels included:
+
+        labels, num_labels = label(mask, background=0, return_num=True)
+
+    elif feature_detection_method == "subsegmentation":
+        print("... using subsegmentation for feature detection")
+
+        # we keep data and just turn sign for minimum criterium
+        if target == "minimum":
+            data_for_detection = -data_i
+            threshold_for_detection = -threshold
         else:
-            selem = np.ones((n_erosion_threshold, n_erosion_threshold))
-        mask = binary_erosion(mask, selem)
-        # detect individual regions, label  and count the number of pixels included:
-    labels, num_labels = label(mask, background=0, return_num=True)
+            data_for_detection = data_i
+            threshold_for_detection = threshold
+
+        kwargs["min_size"] = min_num  # same things, but different option names
+        labels = subsegmentation(data_for_detection, threshold_for_detection, **kwargs)
+
+        # quickfix (dirty and potentially buggy)
+        num_labels = labels.max()
+
     if not is_3D:
         # let's transpose labels to a 1,y,x array to make calculations etc easier.
         labels = labels[np.newaxis, :, :]
@@ -1127,6 +1157,7 @@ def feature_detection_multithreshold_timestep(
             idx_start=idx_start,
             PBC_flag=PBC_flag,
             vertical_axis=vertical_axis,
+            **kwargs,
         )
         if any([x is not None for x in features_threshold_i]):
             features_thresholds = pd.concat(
@@ -1469,6 +1500,7 @@ def feature_detection_multithreshold(
             statistic=statistic,
             statistics_unsmoothed=statistics_unsmoothed,
             return_labels=return_labels,
+            **kwargs,
         )
         # Process the returned data depending on the flags
         if return_labels:

tobac/segmentation/__init__.py

@@ -8,3 +8,5 @@
     watershedding_2D,
     watershedding_3D,
 )
+
+from .geometric_subsegmentation import clustering as subsegmentation