Add missing scripts and gitignores

.gitignore

@@ -10,4 +10,5 @@ slurm/
 scripts/cooper/evaluation_results/
 scripts/cooper/training/copy_testset.py
 scripts/rizzoli/upsample_data.py
-scripts/cooper/training/find_rec_testset.py
+scripts/cooper/training/find_rec_testset.py
+synapse-net-models/

pyproject.toml

@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools>=64.0", "wheel"]
+build-backend = "setuptools.build_meta"

scripts/cooper/analysis/export_vesicles_to_imod.py

@@ -0,0 +1,37 @@
+import os
+from glob import glob
+
+import h5py
+
+from synaptic_reconstruction.imod.to_imod import write_segmentation_to_imod_as_points
+
+
+def export_all_to_imod(check_input=True, check_export=True):
+    files = sorted(glob("./proofread_az/**/*.h5", recursive=True))
+    mrc_root = "./mrc_files"
+    output_folder = "./vesicle_export"
+
+    for ff in files:
+        ds, fname = os.path.split(ff)
+        ds = os.path.basename(ds)
+        out_folder = os.path.join(output_folder, ds)
+        out_path = os.path.join(out_folder, fname.replace(".h5", ".mod"))
+        if os.path.exists(out_path):
+            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["vesicles"][:]
+
+        write_segmentation_to_imod_as_points(mrc_path, seg, out_path, min_radius=7, radius_factor=0.7)
+
+
+def main():
+    export_all_to_imod()
+
+
+if __name__ == "__main__":
+    main()

scripts/cooper/ground_truth/compartments/.gitignore

@@ -0,0 +1 @@
+annotations/

scripts/cryo/vesicles/.gitignore

@@ -1 +1,2 @@
 *.tif
+prediction/

scripts/cryo/vesicles/check_vesicle_gt.py

@@ -0,0 +1,28 @@
+import h5py
+import napari
+
+
+def check_gt(path):
+
+    with h5py.File(path, "r") as f:
+        tomo = f["raw"][:]
+        vesicles = f["labels/vesicles"][:]
+        mask = f["labels/mask"][:]
+
+    v = napari.Viewer()
+    v.add_image(tomo)
+    v.add_labels(vesicles)
+    v.add_labels(mask)
+    napari.run()
+
+
+def main():
+    gt_path1 = "/home/pape/Work/data/fernandez-busnadiego/vesicle_gt/v2/vesicles-33K-L1.h5"
+    check_gt(gt_path1)
+
+    gt_path2 = "/home/pape/Work/data/fernandez-busnadiego/vesicle_gt/v2/vesicles-64K-LAM12.h5"
+    check_gt(gt_path2)
+
+
+if __name__ == "__main__":
+    main()

scripts/cryo/vesicles/debug_vesicle_seg.py

@@ -0,0 +1,35 @@
+import h5py
+import napari
+
+from synaptic_reconstruction.inference.vesicles import distance_based_vesicle_segmentation
+
+
+def debug_vesicle_seg(path, pred_path):
+    with h5py.File(path, "r") as f:
+        raw = f["raw"][:]
+
+    with h5py.File(pred_path, "r") as f:
+        seg = f["/vesicles/segment_from_DA_cryo_v2_masked"][:]
+        fg = f["/prediction_DA_cryo_v2_masked/foreground"][:]
+        bd = f["/prediction_DA_cryo_v2_masked/boundaries"][:]
+
+    vesicles = distance_based_vesicle_segmentation(
+        fg, bd, verbose=True, min_size=500, distance_threshold=4,
+    )
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_image(fg, visible=False)
+    v.add_image(bd, visible=False)
+    v.add_labels(seg)
+    v.add_labels(vesicles)
+    napari.run()
+
+
+def main():
+    path = "/home/pape/Work/data/fernandez-busnadiego/vesicle_gt/v3/vesicles-33K-L1.h5"
+    pred_path = "./prediction/vesicles-33K-L1.h5"
+    debug_vesicle_seg(path, pred_path)
+
+
+main()

scripts/cryo/vesicles/update_vesicle_gt.py

@@ -0,0 +1,71 @@
+import os
+from glob import glob
+
+import h5py
+import napari
+import numpy as np
+from magicgui import magicgui
+
+from skimage.segmentation import watershed
+
+
+INPUT_ROOT = "/home/pape/Work/data/fernandez-busnadiego/vesicle_gt/v2"
+OUTPUT_ROOT = "/home/pape/Work/data/fernandez-busnadiego/vesicle_gt/v3"
+PRED_ROOT = "./prediction"
+
+
+def update_vesicle_gt(path, pred_path):
+    os.makedirs(OUTPUT_ROOT, exist_ok=True)
+    fname = os.path.basename(path)
+    out_path = os.path.join(OUTPUT_ROOT, fname)
+
+    if os.path.exists(out_path):
+        return
+
+    with h5py.File(path, "r") as f:
+        raw = f["raw"][:]
+        gt = f["labels/vesicles"][:]
+        mask = f["labels/mask"][:]
+
+    with h5py.File(pred_path, "r") as f:
+        # pred = f["/vesicles/segment_from_DA_cryo_v2_masked"][:]
+        fg = f["/prediction_DA_cryo_v2_masked/foreground"][:]
+        bd = f["/prediction_DA_cryo_v2_masked/boundaries"][:]
+
+    print("Run watershed")
+    ws_mask = np.logical_or(fg > 0.5, gt != 0)
+    updated_vesicles = watershed(bd, markers=gt, mask=ws_mask)
+    updated_vesicles[mask == 0] = 0
+    print("done")
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_image(fg, visible=False)
+    v.add_image(bd, visible=False)
+    v.add_labels(gt, name="ground-truth")
+    v.add_labels(updated_vesicles, name="updated-gt")
+    v.add_labels(mask, visible=False, name="mask")
+
+    @magicgui(call_button="Save Vesicles")
+    def save_vesicles():
+        with h5py.File(out_path, "a") as f:
+            f.create_dataset("raw", data=raw, compression="gzip")
+            f.create_dataset("labels/vesicles", data=updated_vesicles, compression="gzip")
+            f.create_dataset("labels/mask", data=mask, compression="gzip")
+
+    v.window.add_dock_widget(save_vesicles)
+
+    napari.run()
+
+
+def main():
+    paths = glob(os.path.join(INPUT_ROOT, "*.h5"))
+    for path in paths:
+        fname = os.path.basename(path)
+        pred_path = os.path.join(PRED_ROOT, fname)
+        assert os.path.exists(pred_path)
+        update_vesicle_gt(path, pred_path)
+
+
+if __name__ == "__main__":
+    main()

scripts/data_summary/.gitignore

@@ -0,0 +1,2 @@
+for_zenodo/
+sync_data_for_zenodo.sh

scripts/inner_ear/analysis/check_outliers.py

@@ -0,0 +1,47 @@
+import os
+from pathlib import Path
+
+import imageio.v3 as imageio
+import napari
+from synaptic_reconstruction.file_utils import get_data_path
+from elf.io import open_file
+
+
+tomos_with_outlier = [
+    "WT strong stim/Mouse 1/modiolar/11",
+    "WT strong stim/Mouse 1/pillar/3",
+    "WT strong stim/Mouse 2/modiolar/2",
+    "WT strong stim/Mouse 2/modiolar/2",
+    "WT strong stim/Mouse 2/modiolar/3",
+    "WT control/Mouse 1/modiolar/3",
+    "WT control/Mouse 1/modiolar/3",
+    "WT control/Mouse 2/modiolar/4",
+]
+root = "/home/pape/Work/data/moser/em-synapses/Electron-Microscopy-Susi/Analyse"
+
+for tomo in tomos_with_outlier:
+    folder = os.path.join(root, tomo)
+
+    result_folder = os.path.join(folder, "manuell")
+    if not os.path.exists(result_folder):
+        result_folder = os.path.join(folder, "Manuell")
+    if not os.path.exists(result_folder):
+        continue
+
+    raw_path = get_data_path(folder)
+    with open_file(raw_path, "r") as f:
+        raw = f["data"][:]
+
+    fname = Path(raw_path).stem
+    man_path = os.path.join(result_folder, "Vesikel.tif")
+    seg = imageio.imread(man_path)
+
+    az_path = os.path.join(result_folder, "Membrane.tif")
+    az_seg = imageio.imread(az_path)
+
+    v = napari.Viewer()
+    v.add_image(raw)
+    v.add_labels(seg)
+    v.add_labels(az_seg)
+    v.title = tomo
+    napari.run()

scripts/inner_ear/analysis/check_radii.py

@@ -0,0 +1,105 @@
+import os
+import sys
+from pathlib import Path
+
+import imageio
+import napari
+import pandas as pd
+from elf.io import open_file
+
+from tqdm import tqdm
+from synaptic_reconstruction.file_utils import get_data_path
+
+sys.path.append("../processing")
+
+
+def visualize_folder(folder):
+    result_folder = os.path.join(folder, "manuell")
+    if not os.path.exists(result_folder):
+        result_folder = os.path.join(folder, "Manuell")
+    if not os.path.exists(result_folder):
+        return
+
+    raw_path = get_data_path(folder)
+    with open_file(raw_path, "r") as f:
+        tomo = f["data"][:]
+
+    fname = Path(raw_path).stem
+    auto_path = os.path.join(folder, "automatisch", "v2", f"{fname}_vesicles.h5")
+    man_path = os.path.join(result_folder, "Vesikel.tif")
+
+    with open_file(auto_path, "r") as f:
+        automatic = f["segmentation"][:]
+    manual = imageio.imread(man_path)
+
+    v = napari.Viewer()
+    v.add_image(tomo)
+    v.add_labels(manual)
+    v.add_labels(automatic)
+    v.title = folder
+    napari.run()
+
+
+def visualize_all_radii():
+    from parse_table import parse_table, get_data_root
+
+    data_root = get_data_root()
+    table_path = os.path.join(data_root, "Electron-Microscopy-Susi", "Übersicht.xlsx")
+
+    table = parse_table(table_path, data_root)
+
+    for i, row in tqdm(table.iterrows(), total=len(table)):
+        folder = row["Local Path"]
+        if folder == "":
+            continue
+        visualize_folder(folder)
+
+
+def check_diameter_results():
+    diam_auto = pd.read_excel("./results/vesicle_diameters_tomos_with_manual_annotations.xlsx", sheet_name="Proofread")
+    diam_man = pd.read_excel("./results/vesicle_diameters_tomos_with_manual_annotations.xlsx", sheet_name="Manual")
+
+    print("Summary for the manual measurements:")
+    print(diam_man["diameter [nm]"].mean(), "+-", diam_man["diameter [nm]"].std())
+
+    print("Summary for the auto measurements:")
+    print(diam_auto["diameter [nm]"].mean(), "+-", diam_auto["diameter [nm]"].std())
+
+    print("Unique values")
+    print("Manual:", len(pd.unique(diam_man["diameter [nm]"])), "/", len(diam_man))
+    print("Automatic:", len(pd.unique(diam_auto["diameter [nm]"])), "/", len(diam_auto))
+
+    import matplotlib.pyplot as plt
+    import seaborn as sns
+
+    fig, axes = plt.subplots(2, sharex=True)
+    sns.histplot(data=diam_man, x="diameter [nm]", bins=16, kde=False, ax=axes[0])
+    sns.histplot(data=diam_auto, x="diameter [nm]", bins=16, kde=False, ax=axes[1])
+    plt.show()
+
+
+def test_export():
+    from synaptic_reconstruction.imod.to_imod import write_segmentation_to_imod_as_points
+    from subprocess import run
+
+    mrc_path = "/home/pape/Work/data/moser/em-synapses/Electron-Microscopy-Susi/Analyse/WT strong stim/Mouse 1/pillar/1/Emb71M1aGridA3sec3pil12.rec"  # noqa
+    seg_path = "/home/pape/Work/data/moser/em-synapses/Electron-Microscopy-Susi/Analyse/WT strong stim/Mouse 1/pillar/1/automatisch/v2/Emb71M1aGridA3sec3pil12_vesicles.h5"  # noqa
+    out_path = "exported_vesicles.mod"
+
+    with open_file(seg_path, "r") as f:
+        seg = f["segmentation"][:]
+
+    write_segmentation_to_imod_as_points(
+        mrc_path, seg, out_path, min_radius=10, radius_factor=0.85
+    )
+    run(["imod", mrc_path, out_path])
+
+
+def main():
+    # test_export()
+    check_diameter_results()
+    # visualize_all_radii()
+
+
+if __name__ == "__main__":
+    main()