Merge pull request #50 from computational-cell-analytics/fix_ihc_tonotopic_mapping

Fix tonotopic mapping of IHC segmentation.
Documentation and updated dictionaries have been added for tonotopic mapping of inner hair cells (IHCs) and spiral ganglion neurons (SGNs).
For tonotopic mapping of IHCs, a graph with a reasonably high maximal edge distance is created based on all IHC centroid locations. If this approach fails because the components of the IHC spiral are too far apart, the parameter will need to be adjusted manually. While not ideal, this solution seems satisfactory for the near future.

Author: schilling40

flamingo_tools/segmentation/cochlea_mapping.py

@@ -4,7 +4,6 @@
 import networkx as nx
 import numpy as np
 import pandas as pd
-from networkx.algorithms.approximation import steiner_tree
 from scipy.ndimage import distance_transform_edt, binary_dilation, binary_closing
 from scipy.interpolate import interp1d
 
@@ -146,28 +145,48 @@ def measure_run_length_sgns(centroids: np.ndarray, scale_factor=10):
     return total_distance, path, path_dict
 
 
-def measure_run_length_ihcs(centroids):
+def measure_run_length_ihcs(centroids, max_edge_distance=50):
     """Measure the run lengths of the IHC segmentation
-    by finding the shortest path between the most distant nodes in a Steiner Tree.
+    by determining the shortest path between the most distant nodes of a graph.
+    The graph is created based on a maximal edge distance between nodes.
+    However, this value may not correspond to the max_edge_distance parameter used to process the IHC segmentation,
+    which may consist of more than one component.
+    That's why a large max_edge_distance is used to connect different components and identify
+    the two most distant nodes and the shortest path between them.
+    The path is then edited using a moving average filter.
 
     Args:
         centroids: Centroids of SGN segmentation.
+        max_edge_distance: Maximal edge distance between graph nodes to create an edge between nodes.
 
     Returns:
         Total distance of the path.
         Path as an nd.array of positions.
         A dictionary containing the position and the length fraction of each point in the path.
     """
     graph = nx.Graph()
-    for num, pos in enumerate(centroids):
+    coords = {}
+    labels = [int(i) for i in range(len(centroids))]
+    for index, element in zip(labels, centroids):
+        coords[index] = element
+
+    for num, pos in coords.items():
         graph.add_node(num, pos=pos)
-    # 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)
-    u, v = find_most_distant_nodes(T)
-    path = nx.shortest_path(T, source=u, target=v)
-    total_distance = nx.path_weight(T, path, weight="weight")
+
+    # TODO: Find cleaner option than the creation of a new graph with edges
+    # to identify start / end point and shortest path.
+
+    # 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 <= max_edge_distance:
+                    graph.add_edge(num_i, num_j, weight=dist)
+
+    u, v = find_most_distant_nodes(graph)
+    path = nx.shortest_path(graph, source=u, target=v)
+    total_distance = nx.path_weight(graph, path, weight="weight")
 
     # assign relative distance to points on path
     path_dict = {}
@@ -180,6 +199,9 @@ def measure_run_length_ihcs(centroids):
         path_dict[num + 1] = {"pos": graph.nodes[p]["pos"], "length_fraction": rel_dist}
     path_dict[len(path)] = {"pos": graph.nodes[path[-1]]["pos"], "length_fraction": 1}
 
+    path_pos = np.array([graph.nodes[p]["pos"] for p in path])
+    path = moving_average_3d(path_pos, window=5)
+
     return total_distance, path, path_dict
 
 

reproducibility/block_extraction/Alphatag_MAMDOTOF1L.json

@@ -0,0 +1,54 @@
+[
+	{
+		"cochlea": "M_AMD_OTOF1_L",
+		"image_channel": [
+			"Apha",
+			"Vglut3",
+			"IHC_v4b"
+		],
+		"segmentation_channel": "IHC_v4b",
+		"type": "ihc",
+		"n_blocks": 6,
+		"component_list": [
+			3,
+			11
+		],
+		"halo_size": [
+			256,
+			256,
+			128
+		],
+		"crop_centers": [
+			[
+				981,
+				1326,
+				947
+			],
+			[
+				1281,
+				731,
+				699
+			],
+			[
+				645,
+				872,
+				377
+			],
+			[
+				614,
+				1126,
+				1048
+			],
+			[
+				1017,
+				453,
+				1106
+			],
+			[
+				715,
+				37,
+				527
+			]
+		]
+	}
+]

reproducibility/block_extraction/repro_equidistant_centers.py

@@ -32,7 +32,6 @@ def repro_equidistant_centers(
 
     for dic in param_dicts:
         cochlea = dic["cochlea"]
-        img_channel = dic["image_channel"]
         seg_channel = dic["segmentation_channel"]
 
         s3_path = os.path.join(f"{cochlea}", "tables", f"{seg_channel}", "default.tsv")
@@ -50,8 +49,7 @@ def repro_equidistant_centers(
 
         centers = equidistant_centers(table, component_label=component_list, cell_type=cell_type, n_blocks=n_blocks)
         centers = [[int(c) for c in center] for center in centers]
-        ddict = {"cochlea": cochlea}
-        ddict["image_channel"] = img_channel
+        ddict = dic.copy()
         ddict["crop_centers"] = centers
         ddict["halo_size"] = halo_size
         out_dict.append(ddict)

reproducibility/tonotopic_mapping/2025-07-IHC_fig2.json

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

reproducibility/tonotopic_mapping/2025-07-SGN.json

@@ -0,0 +1,47 @@
+[
+	{
+		"cochlea": "M_LR_000143_L",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000144_L",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000145_L",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000189_L",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000144_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000145_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000153_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000155_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	},
+	{
+		"cochlea": "M_LR_000189_R",
+		"segmentation_channel": "SGN_v2",
+		"type": "sgn"
+	}
+]

reproducibility/tonotopic_mapping/2025-07-SGN_ChReef.json

@@ -2,26 +2,21 @@
 	{
 		"cochlea": "M_LR_000226_L",
 		"segmentation_channel": "SGN_v2",
-		"type": "sgn",
-		"filter_factor": 0.75
+		"type": "sgn"
 	},
 	{
 		"cochlea": "M_LR_000226_R",
 		"segmentation_channel": "SGN_v2",
-		"type": "sgn",
-		"filter_factor": 0.75
+		"type": "sgn"
 	},
 	{
 		"cochlea": "M_LR_000227_L",
 		"segmentation_channel": "SGN_v2",
-		"type": "sgn",
-		"max_edge_distance": 70,
-		"filter_factor": 0.75
+		"type": "sgn"
 	},
 	{
 		"cochlea": "M_LR_000227_R",
 		"segmentation_channel": "SGN_v2",
-		"type": "sgn",
-		"filter_factor": 0.75
+		"type": "sgn"
 	}
 ]