Add volume analysis

Author: constantinpape

flamingo_tools/measurements.py

@@ -43,7 +43,7 @@ def _get_bounding_box_and_center(table, seg_id, resolution, shape, dilation):
     if dilation is not None and dilation > 0:
         bb_extension = dilation + 1
     else:
-        bb_extension = 1
+        bb_extension = 2
 
     bb_min = np.array([
         row.bb_min_z.item(), row.bb_min_y.item(), row.bb_min_x.item()
@@ -166,6 +166,21 @@ def _default_object_features(
     return measures
 
 
+def _morphology_features(seg_id, table, image, segmentation, resolution, **kwargs):
+    measures = {"label_id": seg_id}
+
+    bb, center = _get_bounding_box_and_center(table, seg_id, resolution, image.shape, dilation=0)
+    mask = segmentation[bb] == seg_id
+
+    # Hard-coded value for LaVision cochleae. This is a hack for the wrong voxel size in MoBIE.
+    # resolution = (3.0, 0.76, 0.76)
+
+    volume, surface = _measure_volume_and_surface(mask, resolution)
+    measures["volume"] = volume
+    measures["surface"] = surface
+    return measures
+
+
 def _regionprops_features(seg_id, table, image, segmentation, resolution, background_mask=None, dilation=None):
     bb, _ = _get_bounding_box_and_center(table, seg_id, resolution, image.shape, dilation)
 
@@ -219,6 +234,7 @@ def get_object_measures_from_table(arr_seg, table, keyword="median"):
     "skimage": _regionprops_features,
     "default_background_norm": partial(_default_object_features, background_radius=75, norm=np.divide),
     "default_background_subtract": partial(_default_object_features, background_radius=75, norm=np.subtract),
+    "morphology": _morphology_features,
 }
 """The different feature functions that are supported in `compute_object_measures` and
 that can be selected via the feature_set argument. Currently this supports:

scripts/measurements/density_analysis.py

@@ -73,8 +73,7 @@ def check_implementation():
 
 def compare_gerbil_cochleae():
     # We need the tonotopic mapping for G_EK_000233_L
-    # cochleae = ["G_EK_000233_L", "G_EK_000049_L", "G_EK_000049_R"]
-    cochleae = ["G_EK_000049_L", "G_EK_000049_R"]
+    cochleae = ["G_EK_000233_L", "G_EK_000049_L", "G_EK_000049_R"]
 
     plt.figure(figsize=(6, 3))
     for cochlea in cochleae:
@@ -88,6 +87,10 @@ def compare_gerbil_cochleae():
     plt.show()
 
 
+# TODO: implement the same for mouse cochleae (healthy vs. opto treatment)
+# also show this in tonotopic mapping
+
+
 def main():
     # check_implementation()
     compare_gerbil_cochleae()

scripts/measurements/volume_analysis.py

@@ -0,0 +1,104 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from flamingo_tools.s3_utils import get_s3_path
+from flamingo_tools.measurements import compute_object_measures
+
+
+def compute_sgn_volume(cochlea, output_path, voxel_spacing, seg_name="SGN_v2"):
+    img_path = f"{cochlea}/images/ome-zarr/PV.ome.zarr"
+    seg_path = f"{cochlea}/images/ome-zarr/{seg_name}.ome.zarr"
+
+    img_path, _ = get_s3_path(img_path)
+    seg_path, _ = get_s3_path(seg_path)
+
+    segmentation_table_path = f"{cochlea}/tables/{seg_name}/default.tsv"
+    feature_set = "morphology"
+    compute_object_measures(
+        image_path=img_path,
+        segmentation_path=seg_path,
+        segmentation_table_path=segmentation_table_path,
+        output_table_path=output_path,
+        resolution=voxel_spacing,
+        feature_set=feature_set,
+        component_list=[1],
+        s3_flag=True,
+        image_key="s0",
+        segmentation_key="s0",
+    )
+
+
+def compute_volumes_flamingo():
+    output_folder = "./data/volumes_flamingo"
+    os.makedirs(output_folder, exist_ok=True)
+    cochleae = ["M_LR_000226_L", "M_LR_000226_R", "M_LR_000227_L", "M_LR_000227_R"]
+    voxel_spacing = (0.38, 0.38, 0.38)
+    for cochlea in cochleae:
+        output_path = os.path.join(output_folder, f"{cochlea}.tsv")
+        if os.path.exists(output_path):
+            continue
+        compute_sgn_volume(cochlea, output_path, voxel_spacing)
+
+
+def compute_volumes_lavision():
+    output_folder = "./data/volumes_lavision"
+    os.makedirs(output_folder, exist_ok=True)
+    cochleae = ["LaVision-M02", "LaVision-M03"]
+    # Note: we used a wrong voxel spacing in MoBIE (1.9 micron in-plane instead of 0.76)
+    # We solved this here in a hacky fashion by hard-coding the resolution for the volume
+    # calculation temporariliy in the measurement function.
+    voxel_spacing = (3.0, 1.887779, 1.887779)
+    # voxel_spacing = (3.0, 0.76, 0.76)
+    for cochlea in cochleae:
+        output_path = os.path.join(output_folder, f"{cochlea}.tsv")
+        if os.path.exists(output_path):
+            continue
+        compute_sgn_volume(cochlea, output_path, voxel_spacing, seg_name="SGN_LOWRES-v5c")
+
+
+def compare_volumes():
+    cochleae_flamingo = ["M_LR_000226_L", "M_LR_000226_R", "M_LR_000227_L", "M_LR_000227_R"]
+    cochleae_lavision = ["LaVision-M02", "LaVision-M03"]
+
+    folder_flamingo = "./data/volumes_flamingo"
+    folder_lavision = "./data/volumes_lavision"
+
+    data_flamingo = []
+    for cochlea in cochleae_flamingo:
+        x = pd.read_csv(os.path.join(folder_flamingo, f"{cochlea}.tsv"), sep="\t")
+        volumes = x["volume"].values
+        volumes = volumes[~np.isnan(volumes)]
+        data_flamingo.append(volumes)
+
+    data_lavision = []
+    for cochlea in cochleae_lavision:
+        x = pd.read_csv(os.path.join(folder_lavision, f"{cochlea}.tsv"), sep="\t")
+        volumes = x["volume"].values
+        volumes = volumes[~np.isnan(volumes)]
+        data_lavision.append(volumes)
+
+    fig, axes = plt.subplots(2, sharey=True)
+
+    ax = axes[0]
+    ax.boxplot(data_flamingo, tick_labels=cochleae_flamingo)
+    ax.set_ylabel("SGN Volume [µm^3]")
+
+    ax = axes[1]
+    ax.boxplot(data_lavision, tick_labels=cochleae_lavision)
+    ax.set_ylabel("SGN Volume [µm^3]")
+
+    plt.show()
+
+
+def main():
+    compute_volumes_flamingo()
+    compute_volumes_lavision()
+
+    compare_volumes()
+
+
+if __name__ == "__main__":
+    main()