Add evaluation scripts for cooper data

Author: constantinpape

scripts/cooper/.gitignore

@@ -1 +1,4 @@
 pwd.txt
+debug/
+mito/
+synapse-examples/

scripts/cooper/analysis/az_postprocessing.py

@@ -0,0 +1,116 @@
+import os
+from glob import glob
+
+import h5py
+import napari
+import numpy as np
+
+from magicgui import magicgui
+from scipy.ndimage import binary_dilation, binary_opening
+from skimage.measure import label
+
+
+def postprocess_az(thin_az_seg):
+    # seg = binary_dilation(thin_az_seg, iterations=1)
+    # seg = binary_opening(seg)
+    seg = label(thin_az_seg)
+
+    ids, sizes = np.unique(seg, return_counts=True)
+    ids, sizes = ids[1:], sizes[1:]
+    seg = seg == ids[np.argmax(sizes)].astype("uint8")
+    return seg
+
+
+def process_az(raw_path, az_path):
+    with h5py.File(raw_path, "r") as f:
+        raw = f["raw"][:]
+
+    with h5py.File(az_path, "r") as f:
+        seg = f["thin_az"][:]
+
+    seg_pp = postprocess_az(seg)
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_labels(seg, opacity=1, visible=True)
+    segl = v.add_labels(seg_pp, opacity=1)
+    segl.new_colormap()
+    v.title = raw_path
+    napari.run()
+
+
+def check_all_postprocessed():
+    raw_paths = sorted(glob(os.path.join("imig_data/**/*.h5"), recursive=True))
+    seg_paths = sorted(glob(os.path.join("az_segmentation/**/*.h5"), recursive=True))
+    assert len(raw_paths) == len(seg_paths)
+    for raw_path, seg_path in zip(raw_paths, seg_paths):
+        process_az(raw_path, seg_path)
+
+
+def proofread_file(raw_path, az_path, out_root):
+    ds, fname = os.path.split(raw_path)
+    ds = os.path.basename(ds)
+
+    out_folder = os.path.join(out_root, ds)
+    os.makedirs(out_folder, exist_ok=True)
+    out_path = os.path.join(out_folder, fname)
+
+    if os.path.exists(out_path):
+        return
+
+    with h5py.File(raw_path, "r") as f:
+        raw = f["raw"][:]
+
+    with h5py.File(az_path, "r") as f:
+        seg = f["thin_az"][:]
+
+    seg_pp = postprocess_az(seg)
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_labels(seg, opacity=1, visible=True, name="original")
+    segl = v.add_labels(seg_pp, opacity=1, name="postprocessed")
+    segl.new_colormap()
+
+    v.title = raw_path
+
+    @magicgui(call_button="Postprocess")
+    def postprocess():
+        seg = v.layers["postprocessed"].data
+        seg = postprocess_az(seg)
+        v.layers["postprocessed"].data = seg
+
+    @magicgui(call_button="Save")
+    def save():
+        seg = v.layers["postprocessed"].data
+        with h5py.File(out_path, "a") as f:
+            f.create_dataset("az_thin_proofread", data=seg, compression="gzip")
+        print("Save done!")
+
+    v.window.add_dock_widget(postprocess)
+    v.window.add_dock_widget(save)
+
+    napari.run()
+
+
+def proofread_az(out_folder):
+    raw_paths = sorted(glob(os.path.join("imig_data/**/*.h5"), recursive=True))
+    seg_paths = sorted(glob(os.path.join("az_segmentation/**/*.h5"), recursive=True))
+    assert len(raw_paths) == len(seg_paths)
+    os.makedirs(out_folder, exist_ok=True)
+    for i, (raw_path, seg_path) in enumerate(zip(raw_paths, seg_paths)):
+        print(i, "/", len(seg_paths))
+        proofread_file(raw_path, seg_path, out_folder)
+
+
+def main():
+    # check_all_postprocessed()
+    # process_az(
+    #     "./imig_data/Munc13DKO/A_M13DKO_060212_DKO1.1_crop.h5",
+    #     "./az_segmentation/Munc13DKO/A_M13DKO_060212_DKO1.1_crop.h5"
+    # )
+    proofread_az("./proofread_az")
+
+
+if __name__ == "__main__":
+    main()

scripts/cooper/analysis/check_size_export.py

@@ -0,0 +1,22 @@
+from elf.io import open_file
+
+
+def test_export():
+    from synaptic_reconstruction.imod.to_imod import write_segmentation_to_imod_as_points
+    from subprocess import run
+
+    mrc_path = "20241108_3D_Imig_DATA_2014/!_M13DKO_TOMO_DATA_Imig2014_mrc-mod-FM/A_M13DKO_080212_CTRL4.8_crop/A_M13DKO_080212_CTRL4.8_crop.mrc"  # noqa
+    seg_path = "imig_data/Munc13DKO/A_M13DKO_080212_CTRL4.8_crop.h5"
+    out_path = "exported_vesicles.mod"
+
+    with open_file(seg_path, "r") as f:
+        seg = f["vesicles/segment_from_combined_vesicles"][:]
+
+    # !!!! 0.7
+    write_segmentation_to_imod_as_points(
+        mrc_path, seg, out_path, min_radius=10, radius_factor=0.7
+    )
+    run(["imod", mrc_path, out_path])
+
+
+test_export()

scripts/cooper/analysis/correct_manual_azs.py

@@ -0,0 +1,48 @@
+import os
+
+import h5py
+import napari
+
+from magicgui import magicgui
+
+
+def correct_manual_az(raw_path, seg_path):
+    with h5py.File(raw_path, "r") as f:
+        raw = f["raw"][:]
+
+    seg_key = "az_thin_proofread"
+    with h5py.File(seg_path, "r") as f:
+        seg = f[seg_key][:]
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_labels(seg)
+
+    @magicgui(call_button="save")
+    def save():
+        seg = v.layers["seg"].data
+        with h5py.File(seg_path, "a") as f:
+            f[seg_key][:] = seg
+
+    v.window.add_dock_widget(save)
+
+    napari.run()
+
+
+def main():
+    to_correct = [
+        # ("Munc13DKO", "B_M13DKO_080212_CTRL4.8_crop"),
+        # ("SNAP25", "A_SNAP25_12082_KO1.2_6_crop"),
+        # ("SNAP25", "B_SNAP25_120812_CTRL1.3_13_crop"),
+        ("SNAP25", "B_SNAP25_12082_KO1.2_6_crop")
+    ]
+    for ds, fname in to_correct:
+        raw_path = os.path.join("imig_data", ds, f"{fname}.h5")
+        seg_path = os.path.join("proofread_az", ds, f"{fname}.h5")
+        assert os.path.exists(raw_path)
+        assert os.path.exists(seg_path)
+        correct_manual_az(raw_path, seg_path)
+
+
+if __name__ == "__main__":
+    main()

scripts/cooper/analysis/export_az_to_imod.py

@@ -16,10 +16,12 @@ def check_imod(tomo_path, mod_path):
 
 
 def export_all_to_imod(check_input=True, check_export=True):
-    files = sorted(glob("./az_segmentation/**/*.h5", recursive=True))
+    files = sorted(glob("./proofread_az/**/*.h5", recursive=True))
     mrc_root = "./mrc_files"
     output_folder = "./az_export/initial_model"
 
+    ratings = pd.read_excel("quality_ratings/az_quality_clean_FM.xlsx")
+
     for ff in files:
         ds, fname = os.path.split(ff)
         ds = os.path.basename(ds)
@@ -28,12 +30,26 @@ def export_all_to_imod(check_input=True, check_export=True):
         if os.path.exists(out_path):
             continue
 
+        restrict_to_good_azs = False
+        if restrict_to_good_azs:
+            tomo_name = os.path.splitext(fname)[0]
+            rating = ratings[
+                (ratings["Dataset"] == ds) & (ratings["Tomogram"] == tomo_name)
+            ]["Rating"].values[0]
+            if rating != "Good":
+                print("Skipping", ds, tomo_name, "due to", rating)
+                continue
+
         os.makedirs(out_folder, exist_ok=True)
         mrc_path = os.path.join(mrc_root, ds, fname.replace(".h5", ".rec"))
         assert os.path.exists(mrc_path), mrc_path
 
         with h5py.File(ff, "r") as f:
-            seg = f["thin_az"][:]
+            if "thin_az_corrected" in f:
+                print("Loading corrected az!")
+                seg = f["thin_az_corrected"][:]
+            else:
+                seg = f["az_thin_proofread"][:]
 
         seg = binary_dilation(seg, iterations=2)
         seg = binary_closing(seg, iterations=2)
@@ -114,21 +130,21 @@ def measure_surfaces():
             result["AZ Surface"].append(area)
 
     result = pd.DataFrame(result)
-    result.to_excel("./az_measurements_all.xlsx", index=False)
+    result.to_excel("./results/az_areas_all.xlsx", index=False)
 
 
 def filter_surfaces():
-    all_results = pd.read_excel("./az_measurements_all.xlsx")
+    all_results = pd.read_excel("./results/az_areas_all.xlsx")
     man_tomos = pd.read_csv("./man_tomos.tsv")
 
     man_results = all_results.merge(man_tomos[["Dataset", "Tomogram"]], on=["Dataset", "Tomogram"], how="inner")
-    man_results.to_excel("./az_measuerements_manual.xlsx", index=False)
+    man_results.to_excel("./results/az_areas_manual.xlsx", index=False)
 
 
 def main():
     export_all_to_imod(False, False)
     smooth_all_surfaces(False)
-    # measure_surfaces()
+    measure_surfaces()
     filter_surfaces()
 
 

scripts/cooper/analysis/filter_vesicle_sizes.py

@@ -0,0 +1,50 @@
+import os
+from glob import glob
+
+import numpy as np
+import pandas as pd
+
+
+def filter_sizes_by_distance(size_table, distance_table, out_dir, max_distance=100):
+    fname = os.path.basename(size_table)
+    print("Filtering vesicles for", fname)
+
+    size_table = pd.read_csv(size_table)
+    distance_table = pd.read_csv(distance_table)
+    assert (size_table.columns == distance_table.columns).all()
+    out_columns = {}
+    n_tot, n_filtered = 0, 0
+    all_values = []
+    for col_name in size_table.columns:
+        size_values = size_table[col_name].values
+        distance_values = distance_table[col_name].values
+        size_values, distance_values = (
+            size_values[np.isfinite(size_values)],
+            distance_values[np.isfinite(distance_values)]
+        )
+        assert len(size_values) == len(distance_values)
+        n_tot += len(size_values)
+        mask = distance_values < max_distance
+        out_columns[col_name] = size_values[mask]
+        n_filtered += mask.sum()
+        all_values.extend(size_values[mask].tolist())
+
+    print("Total number of vesicles:", n_tot)
+    print("Number of vesicles after filtering:", n_filtered)
+    print("Average diameter:", np.mean(all_values))
+    os.makedirs(out_dir, exist_ok=True)
+    out_path = os.path.join(out_dir, fname)
+
+    filtered_sizes = pd.DataFrame.from_dict(out_columns, orient='index').transpose()
+    filtered_sizes.to_csv(out_path, index=False)
+
+
+def main():
+    size_tables = sorted(glob("./results/diameters/*.csv"))
+    distance_tables = sorted(glob("./results/distances/*.csv"))
+    for size_tab, distance_tab in zip(size_tables, distance_tables):
+        filter_sizes_by_distance(size_tab, distance_tab, "./results/filtered_diameters")
+
+
+if __name__ == "__main__":
+    main()