Tonotopic mapping for SGN (filtered) and IHC

Author: schilling40

flamingo_tools/segmentation/cochlea_mapping.py

@@ -1,12 +1,15 @@
 import math
+from typing import List, Optional, Tuple
 
 import networkx as nx
+import numpy as np
+import pandas as pd
 from networkx.algorithms.approximation import steiner_tree
 
 from flamingo_tools.segmentation.postprocessing import graph_connected_components
 
 
-def find_most_distant_nodes(G, weight='weight'):
+def find_most_distant_nodes(G: nx.classes.graph.Graph, weight: str = 'weight') -> Tuple[float, float]:
     all_lengths = dict(nx.all_pairs_dijkstra_path_length(G, weight=weight))
     max_dist = 0
     farthest_pair = (None, None)
@@ -21,32 +24,29 @@ def find_most_distant_nodes(G, weight='weight'):
     return u, v
 
 
-def steiner_path_between_distant_nodes(G, weight='weight'):
-    # Step 1: Find the most distant pair of nodes
-    u, v = find_most_distant_nodes(G, weight=weight)
-    terminals = set(G.nodes())  # All nodes are required
-
-    # Step 2: Approximate Steiner Tree over all nodes
-    T = steiner_tree(G, terminals, weight=weight)
-
-    # Step 3: Find the shortest path between u and v in the Steiner Tree
-    path = nx.shortest_path(T, source=u, target=v, weight=weight)
-    total_weight = nx.path_weight(T, path, weight=weight)
-
-    return {
-        "start": u,
-        "end": v,
-        "path": path,
-        "total_weight": total_weight,
-        "steiner_tree": T
-    }
-
-
-def tonotopic_mapping(table, component_label=[1], min_edge_distance=30, min_component_length=50,
-                      cell_type="ihc", weight='weight'):
+def tonotopic_mapping(
+    table: pd.DataFrame,
+    component_label: List[int] = [1],
+    max_edge_distance: float = 30,
+    min_component_length: int = 50,
+    cell_type: str = "ihc",
+    filter_factor: Optional[float] = None
+) -> pd.DataFrame:
     """Tonotopic mapping of IHCs by supplying a table with component labels.
     The mapping assigns a tonotopic label to each IHC according to the position along the length of the cochlea.
+
+    Args:
+        table: Dataframe of segmentation table.
+        component_label: List of component labels to evaluate.
+        max_edge_distance: Maximal edge distance to connect nodes.
+        min_component_length: Minimal number of nodes in component.
+        cell_type: Cell type of segmentation.
+        Filter factor: Fraction of nodes to remove before mapping.
+
+    Returns:
+        Table with tonotopic label for cells.
     """
+    weight = "weight"
     # subset of centroids for given component label(s)
     new_subset = table[table["component_labels"].isin(component_label)]
     comp_label_ids = list(new_subset["label_id"])
@@ -58,12 +58,25 @@ def tonotopic_mapping(table, component_label=[1], min_edge_distance=30, min_comp
     for index, element in zip(labels_subset, centroids_subset):
         coords[index] = element
 
-    components, graph = graph_connected_components(coords, min_edge_distance, min_component_length)
+    _, graph = graph_connected_components(coords, max_edge_distance, min_component_length)
 
-    # approximate Steiner tree and find shortest path between the two most distant nodes
+    unfiltered_graph = graph.copy()
+
+    if filter_factor is not None:
+        if 0 < filter_factor < 1:
+            rng = np.random.default_rng(seed=1234)
+            original_array = np.array(comp_label_ids)
+            target_length = int(len(original_array) * filter_factor)
+            filtered_list = list(rng.choice(original_array, size=target_length, replace=False))
+            for filter_id in filtered_list:
+                graph.remove_node(filter_id)
+        else:
+            raise ValueError(f"Invalid filter factor {filter_factor}. Choose a filter factor between 0 and 1.")
 
     u, v = find_most_distant_nodes(graph)
-    if cell_type == "ihc":
+
+    if not nx.has_path(graph, source=u, target=v) or cell_type == "ihc":
+        # approximate Steiner tree and find shortest path between the two most distant nodes
         terminals = set(graph.nodes())  # All nodes are required
         # Approximate Steiner Tree over all nodes
         T = steiner_tree(graph, terminals, weight=weight)
@@ -86,7 +99,7 @@ def tonotopic_mapping(table, component_label=[1], min_edge_distance=30, min_comp
     path_list[path[-1]] = {"label_id": path[-1], "tonotopic": 1}
 
     # add missing nodes from component
-    pos = nx.get_node_attributes(graph, 'pos')
+    pos = nx.get_node_attributes(unfiltered_graph, 'pos')
     for c in comp_label_ids:
         if c not in path:
             min_dist = float('inf')
@@ -102,8 +115,8 @@ def tonotopic_mapping(table, component_label=[1], min_edge_distance=30, min_comp
 
     tonotopic = [0 for _ in range(len(table))]
     # be aware of 'label_id' of dataframe starting at 1
-    for d in path_list:
-        tonotopic[d["label_id"] - 1] = d["value"] * total_distance
+    for key in list(path_list.keys()):
+        tonotopic[int(path_list[key]["label_id"] - 1)] = path_list[key]["tonotopic"] * total_distance
 
     table.loc[:, "tonotopic_label"] = tonotopic
 

flamingo_tools/segmentation/postprocessing.py

@@ -319,12 +319,12 @@ def downscaled_centroids(
     return new_array
 
 
-def graph_connected_components(coords: dict, min_edge_distance: float, min_component_length: int):
+def graph_connected_components(coords: dict, max_edge_distance: float, min_component_length: int):
     """Create a list of IDs for each connected component of a graph.
 
     Args:
         coords: Dictionary containing label IDs as keys and their position as value.
-        min_edge_distance: Maximal edge distance between graph nodes to create an edge between nodes.
+        max_edge_distance: Maximal edge distance between graph nodes to create an edge between nodes.
         min_component_length: Minimal length of nodes of connected component. Filtered out if lower.
 
     Returns:
@@ -335,12 +335,12 @@ def graph_connected_components(coords: dict, min_edge_distance: float, min_compo
     for num, pos in coords.items():
         graph.add_node(num, pos=pos)
 
-    # create edges between points whose distance is less than threshold min_edge_distance
+    # create edges between points whose distance is less than threshold max_edge_distance
     for num_i, pos_i in coords.items():
         for num_j, pos_j in coords.items():
             if num_i < num_j:
                 dist = math.dist(pos_i, pos_j)
-                if dist <= min_edge_distance:
+                if dist <= max_edge_distance:
                     graph.add_edge(num_i, num_j, weight=dist)
 
     components = list(nx.connected_components(graph))
@@ -360,7 +360,7 @@ def components_sgn(
     keyword: str = "distance_nn100",
     threshold_erode: Optional[float] = None,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
     iterations_erode: Optional[int] = None,
     postprocess_threshold: Optional[float] = None,
     postprocess_components: Optional[List[int]] = None,
@@ -372,7 +372,7 @@ def components_sgn(
         keyword: Keyword of the dataframe column for erosion.
         threshold_erode: Threshold of column value after erosion step with spatial statistics.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
         iterations_erode: Number of iterations for erosion, normally determined automatically.
         postprocess_threshold: Post-process graph connected components by searching for points closer than threshold.
         postprocess_components: Post-process specific graph connected components ([0] for largest component only).
@@ -412,7 +412,7 @@ def components_sgn(
     for index, element in zip(labels_subset, centroids_subset):
         coords[index] = element
 
-    components, _ = graph_connected_components(coords, min_edge_distance, min_component_length)
+    components, _ = graph_connected_components(coords, max_edge_distance, min_component_length)
 
     length_components = [len(c) for c in components]
     length_components, components = zip(*sorted(zip(length_components, components), reverse=True))
@@ -448,7 +448,7 @@ def label_components_sgn(
     min_size: int = 1000,
     threshold_erode: Optional[float] = None,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
     iterations_erode: Optional[int] = None,
     postprocess_threshold: Optional[float] = None,
     postprocess_components: Optional[List[int]] = None,
@@ -460,7 +460,7 @@ def label_components_sgn(
         min_size: Minimal number of pixels for filtering small instances.
         threshold_erode: Threshold of column value after erosion step with spatial statistics.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
         iterations_erode: Number of iterations for erosion, normally determined automatically.
         postprocess_threshold: Post-process graph connected components by searching for points closer than threshold.
         postprocess_components: Post-process specific graph connected components ([0] for largest component only).
@@ -474,7 +474,7 @@ def label_components_sgn(
     table = table[table.n_pixels >= min_size]
 
     components = components_sgn(table, threshold_erode=threshold_erode, min_component_length=min_component_length,
-                                min_edge_distance=min_edge_distance, iterations_erode=iterations_erode,
+                                max_edge_distance=max_edge_distance, iterations_erode=iterations_erode,
                                 postprocess_threshold=postprocess_threshold,
                                 postprocess_components=postprocess_components)
 
@@ -496,7 +496,7 @@ def postprocess_sgn_seg(
     min_size: int = 1000,
     threshold_erode: Optional[float] = None,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
     iterations_erode: Optional[int] = None,
 ) -> pd.DataFrame:
     """Postprocessing SGN segmentation of cochlea.
@@ -506,7 +506,7 @@ def postprocess_sgn_seg(
         min_size: Minimal number of pixels for filtering small instances.
         threshold_erode: Threshold of column value after erosion step with spatial statistics.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
         iterations_erode: Number of iterations for erosion, normally determined automatically.
 
     Returns:
@@ -515,7 +515,7 @@ def postprocess_sgn_seg(
 
     comp_labels = label_components_sgn(table, min_size=min_size, threshold_erode=threshold_erode,
                                        min_component_length=min_component_length,
-                                       min_edge_distance=min_edge_distance, iterations_erode=iterations_erode)
+                                       max_edge_distance=max_edge_distance, iterations_erode=iterations_erode)
 
     table.loc[:, "component_labels"] = comp_labels
 
@@ -525,14 +525,14 @@ def postprocess_sgn_seg(
 def components_ihc(
     table: pd.DataFrame,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
 ):
     """Create connected components for IHC segmentation.
 
     Args:
         table: Dataframe of segmentation table.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
 
     Returns:
         Subgraph components as lists of label_ids of dataframe.
@@ -543,23 +543,23 @@ def components_ihc(
     for index, element in zip(labels, centroids):
         coords[index] = element
 
-    components, _ = graph_connected_components(coords, min_edge_distance, min_component_length)
+    components, _ = graph_connected_components(coords, max_edge_distance, min_component_length)
     return components
 
 
 def label_components_ihc(
     table: pd.DataFrame,
     min_size: int = 1000,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
 ) -> List[int]:
     """Label components using graph connected components.
 
     Args:
         table: Dataframe of segmentation table.
         min_size: Minimal number of pixels for filtering small instances.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
 
     Returns:
         List of component label for each point in dataframe. 0 - background, then in descending order of size
@@ -570,7 +570,7 @@ def label_components_ihc(
     table = table[table.n_pixels >= min_size]
 
     components = components_ihc(table, min_component_length=min_component_length,
-                                min_edge_distance=min_edge_distance)
+                                max_edge_distance=max_edge_distance)
 
     # add size-filtered objects to have same initial length
     table = pd.concat([table, entries_filtered], ignore_index=True)
@@ -592,23 +592,23 @@ def postprocess_ihc_seg(
     table: pd.DataFrame,
     min_size: int = 1000,
     min_component_length: int = 50,
-    min_edge_distance: float = 30,
+    max_edge_distance: float = 30,
 ) -> pd.DataFrame:
     """Postprocessing IHC segmentation of cochlea.
 
     Args:
         table: Dataframe of segmentation table.
         min_size: Minimal number of pixels for filtering small instances.
         min_component_length: Minimal length for filtering out connected components.
-        min_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
+        max_edge_distance: Maximal distance in micrometer between points to create edges for connected components.
 
     Returns:
         Dataframe with component labels.
     """
 
     comp_labels = label_components_ihc(table, min_size=min_size,
                                        min_component_length=min_component_length,
-                                       min_edge_distance=min_edge_distance)
+                                       max_edge_distance=max_edge_distance)
 
     table.loc[:, "component_labels"] = comp_labels
 

reproducibility/tonotopic_mapping/2025-07-IHC_fig2.json

@@ -0,0 +1,22 @@
+[
+	{
+		"cochlea": "M_LR_000226_L",
+		"segmentation_channel": "IHC_v3",
+		"type": "ihc"
+	},
+	{
+		"cochlea": "M_LR_000226_R",
+		"segmentation_channel": "IHC_v3",
+		"type": "ihc"
+	},
+	{
+		"cochlea": "M_LR_000227_L",
+		"segmentation_channel": "IHC_v3",
+		"type": "ihc"
+	},
+	{
+		"cochlea": "M_LR_000227_R",
+		"segmentation_channel": "IHC_v3",
+		"type": "ihc"
+	}
+]

reproducibility/tonotopic_mapping/2025-07-SGN.json

@@ -0,0 +1,57 @@
+[
+	{
+		"cochlea": "M_AMD_000058_L",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000144_L",
+		"segmentation_channel": "SGN_resized_v2",
+		"max_edge_distance": 70,
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000144_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000145_L",
+		"segmentation_channel": "SGN_resized_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000151_R",
+		"segmentation_channel": "SGN_resized_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000155_L",
+		"segmentation_channel": "SGN_resized_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000155_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000167_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	},
+	{
+		"cochlea": "M_LR_000184_L",
+		"segmentation_channel": "SGN_resized_v2",
+		"type": "sgn",
+		"filter_factor": 0.75
+	}
+]