Evaluate LaVision-rbOtof thresholds

Author: schilling40

flamingo_tools/segmentation/chreef_utils.py

@@ -73,6 +73,9 @@ def get_roi(coord: tuple, roi_halo: tuple, resolution: float = 0.38) -> Tuple[in
     # reverse dimensions for correct extraction
     coords.reverse()
     coords = np.array(coords)
+    if not isinstance(resolution, float):
+        assert len(resolution) == 3
+        resolution = np.array(resolution)[::-1]
     coords = coords / resolution
     coords = np.round(coords).astype(np.int32)
 
@@ -144,12 +147,13 @@ def find_inbetween_ids(
     return inbetween_ids, allweak_positives, negexc_negatives
 
 
-def get_median_intensity(file_negexc, file_allweak, center, data_seg, table, column="median"):
+def get_median_intensity(file_negexc, file_allweak, center, data_seg, table, column="median",
+                         resolution=0.38):
     arr_negexc = tifffile.imread(file_negexc)
     arr_allweak = tifffile.imread(file_allweak)
 
     roi_halo = tuple([r // 2 for r in arr_negexc.shape])
-    roi = get_roi(center, roi_halo)
+    roi = get_roi(center, roi_halo, resolution=resolution)
 
     roi_seg = data_seg[roi]
     inbetween_ids, allweak_positives, negexc_negatives = find_inbetween_ids(arr_negexc, arr_allweak, roi_seg)
@@ -172,7 +176,8 @@ def get_median_intensity(file_negexc, file_allweak, center, data_seg, table, col
     return np.median(list(intensities))
 
 
-def localize_median_intensities(annotation_dir, cochlea, data_seg, table_measure, column="median", pattern=None):
+def localize_median_intensities(annotation_dir, cochlea, data_seg, table_measure, column="median", pattern=None,
+                                resolution=0.38):
     """Find median intensities in blocks and assign them to center positions of cropped block.
     """
     annotation_dic = find_annotations(annotation_dir, cochlea, pattern=pattern)
@@ -184,7 +189,7 @@ def localize_median_intensities(annotation_dir, cochlea, data_seg, table_measure
         file_pos = annotation_dic[center_str]["file_pos"]
         file_neg = annotation_dic[center_str]["file_neg"]
         median_intensity = get_median_intensity(file_neg, file_pos, center_coord, data_seg,
-                                                table_measure, column=column)
+                                                table_measure, column=column, resolution=resolution)
 
         if median_intensity is None:
             print(f"No threshold identified for {center_str}.")

flamingo_tools/segmentation/cochlea_mapping.py

@@ -314,8 +314,8 @@ def measure_run_length_ihcs_multi_component(
             print(f"Graph consists of {len(components)} connected components.")
             if len(component_label) != len(components):
                 raise ValueError(f"Length of graph components {len(components)} "
-                                f"does not match number of component labels {len(component_label)}. "
-                                "Check max_edge_distance and post-processing.")
+                                 f"does not match number of component labels {len(component_label)}. "
+                                 "Check max_edge_distance and post-processing.")
 
             # Order connected components in order of component labels
             # e.g. component_labels = [7, 4, 1, 11] and len_c = [600, 400, 300, 55]
@@ -410,6 +410,7 @@ def measure_run_length_ihcs_multi_component(
 
     return total_distance, path, path_dict
 
+
 def measure_run_length_ihcs(
     centroids: np.ndarray,
     max_edge_distance: float = 30,
@@ -583,7 +584,7 @@ def get_centers_from_path(
     try:
         f = interp1d(cum_len, path, axis=0)  # fill_value="extrapolate"
         centers = f(target_s)
-    except ValueError as ve:
+    except ValueError:
         print("Using extrapolation to fill values.")
         f = interp1d(cum_len, path, axis=0, fill_value="extrapolate")
         centers = f(target_s)
@@ -724,6 +725,7 @@ def tonotopic_mapping(
     component_mapping: Optional[List[int]] = None,
     cell_type: str = "ihc",
     animal: str = "mouse",
+    max_edge_distance: float = 30,
     apex_higher: bool = True,
 ) -> pd.DataFrame:
     """Tonotopic mapping of IHCs by supplying a table with component labels.
@@ -749,6 +751,7 @@ def tonotopic_mapping(
     if cell_type == "ihc":
         total_distance, _, path_dict = measure_run_length_ihcs(
             centroids, component_label=component_label, apex_higher=apex_higher,
+            max_edge_distance=max_edge_distance
         )
 
     else:

reproducibility/object_measures/repro_object_measures.py

@@ -82,8 +82,8 @@ def repro_object_measures(
                     )
 
                 compute_object_measures(
-                    image_path=img_path,
-                    segmentation_path=seg_path,
+                    image_path=img_s3,
+                    segmentation_path=seg_s3,
                     segmentation_table_path=seg_table_s3,
                     output_table_path=output_table_path,
                     image_key=input_key,

reproducibility/tonotopic_mapping/LaVision_OTOF23R.json

@@ -0,0 +1,14 @@
+[
+	{
+		"cochlea": "LaVision-OTOF23R",
+		"segmentation_channel": "IHC_LOWRES-v3",
+		"component_list":
+		[
+			4,
+			18, 
+			7
+		],
+		"max_edge_distance": 50,
+		"type": "ihc"
+	}
+]

reproducibility/tonotopic_mapping/LaVision_OTOF25R.json

@@ -0,0 +1,12 @@
+[	
+	{
+		"cochlea": "LaVision-OTOF25R",
+		"segmentation_channel": "IHC_LOWRES-v3",
+		"component_list":
+		[
+			1
+		],
+		"max_edge_distance": 300,
+		"type": "ihc"
+	}
+]

reproducibility/tonotopic_mapping/MLR184L.json

@@ -4,7 +4,7 @@
 		"component_list": [
 			1
 		],
-		"segmentation_channel": "SGN_v2",
+		"segmentation_channel": "SGN_v2b",
 		"type": "sgn"
 	}
 ]									

reproducibility/tonotopic_mapping/MLR184R.json

@@ -4,7 +4,7 @@
 		"component_list": [
 			1
 		],
-		"segmentation_channel": "SGN_v2",
+		"segmentation_channel": "SGN_v2b",
 		"type": "sgn"
 	}
 ]									

reproducibility/tonotopic_mapping/repro_tonotopic_mapping.py

@@ -18,6 +18,7 @@ def repro_tonotopic_mapping(
 ):
     default_cell_type = "ihc"
     default_apex_position = "apex_higher"
+    default_max_edge_distance = 30
     default_component_list = [1]
 
     remove_columns = ["tonotopic_label",
@@ -39,7 +40,8 @@ def repro_tonotopic_mapping(
         elif cochlea[0] in ["G", "g"]:
             animal = "gerbil"
         else:
-            raise ValueError("Cochlea does not have expected name format 'M_[...]' or 'G_[...]'.")
+            animal = "mouse"
+            # raise ValueError("Cochlea does not have expected name format 'M_[...]' or 'G_[...]'.")
 
         cochlea_str = "-".join(cochlea.split("_"))
         seg_str = "-".join(seg_channel.split("_"))
@@ -58,6 +60,7 @@ def repro_tonotopic_mapping(
         component_list = dic["component_list"] if "component_list" in dic else default_component_list
         component_mapping = dic["component_mapping"] if "component_mapping" in dic else component_list
         apex_position = dic["apex_position"] if "apex_position" in dic else default_apex_position
+        max_edge_distance = dic["max_edge_distance"] if "max_edge_distance" in dic else default_max_edge_distance
 
         apex_higher = (apex_position == "apex_higher")
 
@@ -68,7 +71,7 @@ def repro_tonotopic_mapping(
         if not os.path.isfile(output_table_path) or force_overwrite:
             table = tonotopic_mapping(table, component_label=component_list, animal=animal,
                                       cell_type=cell_type, component_mapping=component_mapping,
-                                      apex_higher=apex_higher)
+                                      apex_higher=apex_higher, max_edge_distance=max_edge_distance)
 
             table.to_csv(output_table_path, sep="\t", index=False)
 

scripts/measurements/evaluate_marker_annotations.py

@@ -29,34 +29,33 @@
 ]
 
 
-def get_length_fraction_from_center(table, center_str):
+def get_length_fraction_from_center(table, center_str, halo_size=20):
     """Get 'length_fraction' parameter for center coordinate by averaging nearby segmentation instances.
     """
     center_coord = tuple([int(c) for c in center_str.split("-")])
     (cx, cy, cz) = center_coord
-    offset = 20
     subset = table[
-        (cx - offset < table["anchor_x"]) &
-        (table["anchor_x"] < cx + offset) &
-        (cy - offset < table["anchor_y"]) &
-        (table["anchor_y"] < cy + offset) &
-        (cz - offset < table["anchor_z"]) &
-        (table["anchor_z"] < cz + offset)
+        (cx - halo_size < table["anchor_x"]) &
+        (table["anchor_x"] < cx + halo_size) &
+        (cy - halo_size < table["anchor_y"]) &
+        (table["anchor_y"] < cy + halo_size) &
+        (cz - halo_size < table["anchor_z"]) &
+        (table["anchor_z"] < cz + halo_size)
     ]
     length_fraction = list(subset["length_fraction"])
     length_fraction = float(sum(length_fraction) / len(length_fraction))
     return length_fraction
 
 
-def apply_nearest_threshold(intensity_dic, table_seg, table_measurement):
+def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, halo_size=20):
     """Apply threshold to nearest segmentation instances.
     Crop centers are transformed into the "length fraction" parameter of the segmentation table.
     This avoids issues with the spiral shape of the cochlea and maps the assignment onto the Rosenthal"s canal.
     """
     # assign crop centers to length fraction of Rosenthal"s canal
     lf_intensity = {}
     for key in intensity_dic.keys():
-        length_fraction = get_length_fraction_from_center(table_seg, key)
+        length_fraction = get_length_fraction_from_center(table_seg, key, halo_size=halo_size)
         intensity_dic[key]["length_fraction"] = length_fraction
         lf_intensity[length_fraction] = {"threshold": intensity_dic[key]["median_intensity"]}
 
@@ -94,13 +93,14 @@ def apply_nearest_threshold(intensity_dic, table_seg, table_measurement):
     return table_seg
 
 
-def find_thresholds(cochlea_annotations, cochlea, data_seg, table_measurement):
+def find_thresholds(cochlea_annotations, cochlea, data_seg, table_measurement, resolution=0.38):
     # Find the median intensities by averaging the individual annotations for specific crops
     annotation_dics = {}
     annotated_centers = []
     for annotation_dir in cochlea_annotations:
         print(f"Localizing threshold with median intensities for {os.path.basename(annotation_dir)}.")
-        annotation_dic = localize_median_intensities(annotation_dir, cochlea, data_seg, table_measurement)
+        annotation_dic = localize_median_intensities(annotation_dir, cochlea, data_seg, table_measurement,
+                                                     resolution=resolution)
         annotated_centers.extend(annotation_dic["center_strings"])
         annotation_dics[annotation_dir] = annotation_dic
 
@@ -168,6 +168,13 @@ def evaluate_marker_annotation(
     """
     input_key = "s0"
 
+    if marker_name == "rbOtof":
+        halo_size = 150
+        resolution = [1.887779, 1.887779, 3.0]
+    else:
+        halo_size = 20
+        resolution = 0.38
+
     if annotation_dirs is None:
         if "MARKER_DIR" in globals():
             marker_dir = MARKER_DIR
@@ -199,16 +206,17 @@ def evaluate_marker_annotation(
         with fs.open(table_path_s3, "r") as f:
             table_measurement = pd.read_csv(f, sep="\t")
 
-        # Find the threholds from the annotated blocks and save it if specified.
-        intensity_dic = find_thresholds(cochlea_annotations, cochlea, data_seg, table_measurement)
+        # Find the thresholds from the annotated blocks and save it if specified.
+        intensity_dic = find_thresholds(cochlea_annotations, cochlea, data_seg, table_measurement,
+                                        resolution=resolution)
         if threshold_save_dir is not None:
             os.makedirs(threshold_save_dir, exist_ok=True)
             threshold_out_path = os.path.join(threshold_save_dir, f"{cochlea_str}_{marker_name}_{seg_string}.json")
             with open(threshold_out_path, "w") as f:
                 json.dump(intensity_dic, f, sort_keys=True, indent=4)
 
         # Apply the threshold to all SGNs.
-        table_seg = apply_nearest_threshold(intensity_dic, table_seg, table_measurement)
+        table_seg = apply_nearest_threshold(intensity_dic, table_seg, table_measurement, halo_size=halo_size)
 
         # Save the table with positives / negatives for all SGNs.
         os.makedirs(output_dir, exist_ok=True)
@@ -224,12 +232,15 @@ def main():
     parser.add_argument("-o", "--output", type=str, required=True, help="Output directory.")
     parser.add_argument("-a", "--annotation_dirs", type=str, nargs="+", default=None,
                         help="Directories containing marker annotations.")
-    parser.add_argument("--threshold_save_dir", "-t")
+    parser.add_argument("-t", "--threshold_save_dir")
+    parser.add_argument("-s", "--seg_name", type=str, default="SGN_v2")
+    parser.add_argument("-m", "--marker_name", type=str, default="GFP")
     parser.add_argument("-f", "--force", action="store_true")
 
     args = parser.parse_args()
     evaluate_marker_annotation(
-        args.cochlea, args.output, args.annotation_dirs, threshold_save_dir=args.threshold_save_dir, force=args.force
+        args.cochlea, args.output, args.annotation_dirs, threshold_save_dir=args.threshold_save_dir,
+        seg_name=args.seg_name, marker_name=args.marker_name, force=args.force
     )