Merge branch 'main' of https://github.com/computational-cell-analytics/synaptic-reconstruction into more-inner-ear-analysis

Author: constantinpape

environment.yaml

@@ -8,5 +8,9 @@ dependencies:
     - pip
     - pyqt
     - magicgui
+    - pytorch
+    - bioimageio.core
+    - kornia
+    - tensorboard
     - pip:
         - napari-skimage-regionprops

run_correction.sh

@@ -0,0 +1,66 @@
+import os
+from glob import glob
+
+import h5py
+import numpy as np
+
+from scipy.ndimage import binary_closing
+from skimage.measure import label
+from synaptic_reconstruction.ground_truth.shape_refinement import edge_filter
+from tqdm import tqdm
+
+ROOT = "/mnt/lustre-emmy-hdd/projects/nim00007/data/synaptic-reconstruction/cooper/20241102_TOMO_DATA_Imig2014/final_Imig2014_seg_autoComp"  # noqa
+
+OUTPUT_AZ = "./boundary_az"
+
+
+def filter_az(path):
+    # Check if we have the output already.
+    ds, fname = os.path.split(path)
+    ds = os.path.basename(ds)
+    out_path = os.path.join(OUTPUT_AZ, ds, fname)
+    os.makedirs(os.path.join(OUTPUT_AZ, ds), exist_ok=True)
+    if os.path.exists(out_path):
+        return
+
+    with h5py.File(path, "r") as f:
+        raw = f["raw"][:]
+        az = f["AZ/segment_from_AZmodel_v3"][:]
+        vesicles = f["/vesicles/segment_from_combined_vesicles"][:]
+
+    # Compute the sato filter of the raw data, smooth it afterwards.
+    # This will highlight dark ridge-like structures, and so
+    # will yield high values for the plasma membrane.
+    hmap = edge_filter(raw, sigma=1.0, method="sato", per_slice=True, n_threads=8)
+
+    # Filter the active zone by combining a bunch of things:
+    # 1. Find a mask with high values in the ridge filter.
+    threshold_hmap = 0.5
+    az_filtered = hmap > threshold_hmap
+    # 2. Intersect it with the active zone predictions.
+    az_filtered = np.logical_and(az_filtered, az)
+    # 3. Intersect it with the negative vesicle mask.
+    az_filtered = np.logical_and(az_filtered, vesicles == 0)
+
+    # Postprocessing of the filtered active zone:
+    # 1. Apply connected components and only keep the largest component.
+    az_filtered = label(az_filtered)
+    ids, sizes = np.unique(az_filtered, return_counts=True)
+    ids, sizes = ids[1:], sizes[1:]
+    az_filtered = (az_filtered == ids[np.argmax(sizes)]).astype("uint8")
+    # 2. Apply binary closing.
+    az_filtered = np.logical_or(az_filtered, binary_closing(az_filtered, iterations=4)).astype("uint8")
+
+    # Save the result.
+    with h5py.File(out_path, "a") as f:
+        f.create_dataset("filtered_az", data=az_filtered, compression="gzip")
+
+
+def main():
+    files = sorted(glob(os.path.join(ROOT, "**/*.h5"), recursive=True))
+    for ff in tqdm(files):
+        filter_az(ff)
+
+
+if __name__ == "__main__":
+    main()

scripts/aggregate_data_information.py

@@ -8,7 +8,7 @@
 import numpy as np
 import pandas
 
-from synaptic_reconstruction.distance_measurements import create_object_distance_lines
+from synaptic_reconstruction.distance_measurements import create_object_distance_lines, load_distances
 from synaptic_reconstruction.file_utils import get_data_path
 from synaptic_reconstruction.tools.distance_measurement import _downsample
 
@@ -21,11 +21,12 @@
 def get_distance_visualization(
     tomo, segmentations, distance_paths, vesicle_ids, scale, return_mem_props=False
 ):
-    ribbon_lines, _ = create_object_distance_lines(distance_paths["ribbon"], seg_ids=vesicle_ids, scale=scale)
-    pd_lines, _ = create_object_distance_lines(distance_paths["PD"], seg_ids=vesicle_ids, scale=scale)
-    membrane_lines, mem_props = create_object_distance_lines(
-        distance_paths["membrane"], seg_ids=vesicle_ids, scale=scale
-    )
+    d, e1, e2, ids = load_distances(distance_paths["ribbon"])
+    ribbon_lines, _ = create_object_distance_lines(d, e1, e2, ids, filter_seg_ids=vesicle_ids, scale=scale)
+    d, e1, e2, ids = load_distances(distance_paths["PD"])
+    pd_lines, _ = create_object_distance_lines(d, e1, e2, ids, filter_seg_ids=vesicle_ids, scale=scale)
+    d, e1, e2, ids = load_distances(distance_paths["membrane"])
+    membrane_lines, mem_props = create_object_distance_lines(d, e1, e2, ids, filter_seg_ids=vesicle_ids, scale=scale)
 
     distance_lines = {
         "ribbon_distances": ribbon_lines,

scripts/cooper/analysis/active_zone_analysis.py

@@ -10,7 +10,22 @@
 from tqdm import tqdm
 
 ROOT = "/mnt/lustre-emmy-hdd/projects/nim00007/data/synaptic-reconstruction/moser/inner_ear_data"
+ROOT_OTHER_TOMOS = "/mnt/lustre-emmy-hdd/projects/nim00007/data/synaptic-reconstruction/moser/other_tomograms/"
+
 LABEL_KEY = "labels/inner_ear_structures"
+OTHER_NAMES = ["vesicle_pools", "tether", "rat"]
+
+
+def get_other_paths(name):
+    assert name in OTHER_NAMES, f"Invalid name {name}"
+    if name == "vesicle_pools":
+        folder = "01_vesicle_pools"
+    elif name == "tether":
+        folder = "02_tether"
+    else:
+        folder = "03_ratten_tomos"
+    paths = sorted(glob(os.path.join(ROOT, folder, "*.h5")))
+    return paths
 
 
 def get_train_val_test_split(root):
@@ -34,8 +49,7 @@ def get_train_val_test_split(root):
     return train_tomos, val_tomos, test_tomos
 
 
-def preprocess_labels(tomograms):
-    structure_keys = ("ribbon", "PD", "membrane")
+def preprocess_labels(tomograms, structure_keys=("ribbon", "PD", "membrane")):
     nc = len(structure_keys)
 
     for tomo in tqdm(tomograms, desc="Preprocess labels"):
@@ -60,11 +74,11 @@ def noop(x):
     return x
 
 
-def train_inner_ear_structures(train_tomograms, val_tomograms):
+def train_inner_ear_structures(train_tomograms, val_tomograms, name):
     patch_shape = (64, 512, 512)
     sampler = MinForegroundSampler(min_fraction=0.05, p_reject=0.95)
     supervised_training(
-        name="inner_ear_structure_model",
+        name=name,
         train_paths=train_tomograms, val_paths=val_tomograms,
         label_key=LABEL_KEY, patch_shape=patch_shape, save_root=".",
         sampler=sampler, label_transform=noop, out_channels=3,
@@ -73,11 +87,27 @@ def train_inner_ear_structures(train_tomograms, val_tomograms):
     )
 
 
-def main():
+def training_v1():
+    train_tomograms, val_tomograms, _ = get_train_val_test_split(ROOT)
+    preprocess_labels(train_tomograms)
+    preprocess_labels(val_tomograms)
+    train_inner_ear_structures(train_tomograms, val_tomograms, name="inner_ear_structure_model")
+
+
+def training_v2():
     train_tomograms, val_tomograms, _ = get_train_val_test_split(ROOT)
     preprocess_labels(train_tomograms)
+    for name in OTHER_NAMES:
+        other_tomograms = get_other_paths(name)
+        preprocess_labels(other_tomograms, structure_keys=("ribbons", "presynapse", "membrane"))
+        train_tomograms.extend(other_tomograms)
     preprocess_labels(val_tomograms)
-    train_inner_ear_structures(train_tomograms, val_tomograms)
+    train_inner_ear_structures(train_tomograms, val_tomograms, name="inner_ear_structure_model_v2")
+
+
+def main():
+    # training_v1()
+    training_v2()
 
 
 if __name__ == "__main__":

scripts/inner_ear/check_results.py

@@ -8,13 +8,15 @@
     name="synaptic_reconstruction",
     packages=find_packages(exclude=["test"]),
     version=__version__,
-    author="Constantin Pape; Sarah Muth",
+    author="Constantin Pape; Sarah Muth; Luca Freckmann",
     url="https://github.com/computational-cell-analytics/synaptic_reconstruction",
     license="MIT",
     entry_points={
         "console_scripts": [
-            "sr_tools.correct_segmentation = synaptic_reconstruction.tools.segmentation_correction:main",
-            "sr_tools.measure_distances = synaptic_reconstruction.tools.distance_measurement:main",
-        ]
+            "synapse_net.run_segmentation = synaptic_reconstruction.tools.cli:segmentation_cli"
+        ],
+        "napari.manifest": [
+            "synaptic_reconstruction = synaptic_reconstruction:napari.yaml",
+        ],
     },
 )