Merge branch 'synapse-analysis' of https://github.com/computational-cell-analytics/lightsheet-moser into synapse-analysis

constantinpape · constantinpape · commit 180329209397 · 2025-06-30T15:50:28.000+02:00
diff --git a/scripts/validation/SGNs/consensus_step2.py b/scripts/validation/SGNs/consensus_step2.py
@@ -0,0 +1,115 @@
+import os
+from glob import glob
+
+import pandas as pd
+from flamingo_tools.validation import create_consensus_annotations
+
+ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationSGNs/for_consensus_annotation"  # noqa
+
+ANNOTATION_FOLDERS = ["AnnotationsAMD", "AnnotationsEK"]
+COLOR = ["blue", "yellow"]
+CONSENSUS_ANNOTATIONS = "consensus_annotations"
+OUTPUT_FOLDER = os.path.join(ROOT, "final_consensus_annotations")
+
+
+def match_annotations(image_path, annotation_folders):
+    annotations = {}
+    prefix = os.path.basename(image_path).split("_")[:3]
+    prefix = "_".join(prefix)
+
+    annotations = {}
+    for annotation_folder in annotation_folders:
+        all_annotations = glob(os.path.join(ROOT, annotation_folder, "*.csv"))
+        matches = [ann for ann in all_annotations if os.path.basename(ann).startswith(prefix)]
+        if len(matches) != 1:
+            breakpoint()
+        assert len(matches) == 1
+        annotation_path = matches[0]
+        annotations[annotation_folder] = annotation_path
+
+    assert len(annotations) == len(annotation_folders)
+    return annotations
+
+
+def create_consensus_step2(image_path, check):
+    print("Compute consensus annotations for", image_path)
+    annotation_paths = match_annotations(image_path, ANNOTATION_FOLDERS)
+    matching_distance = 8  # TODO
+    consensus_annotations, unmatched_annotations = create_consensus_annotations(
+        annotation_paths, matching_distance=matching_distance, min_matches_for_consensus=2,
+    )
+    fname = os.path.basename(image_path)
+
+    # We exclude the two images for which we don't have previous consensus annotations yet.
+    # (This was due to resizing problems.)
+    if fname in (
+        "MLR169R_PV_z1913_base_full_rescaled.tif",
+        "MLR169R_PV_z2594_mid_full_rescaled.tif",
+    ):
+        prev_consensus = None
+    else:
+        prev_consensus = match_annotations(image_path, [CONSENSUS_ANNOTATIONS])[CONSENSUS_ANNOTATIONS]
+        prev_consensus = pd.read_csv(prev_consensus)[["axis-0", "axis-1", "axis-2"]]
+
+    if check:
+        import napari
+        import tifffile
+
+        consensus_annotations = consensus_annotations[["axis-0", "axis-1", "axis-2"]].values
+        unmatched_annotators = unmatched_annotations.annotator.values
+        unmatched_annotations = unmatched_annotations[["axis-0", "axis-1", "axis-2"]].values
+
+        image = tifffile.imread(image_path)
+        v = napari.Viewer()
+        v.add_image(image)
+        if prev_consensus is not None:
+            v.add_points(prev_consensus.values, face_color="gray", name="previous-consensus-annotations")
+        v.add_points(consensus_annotations, face_color="green")
+        v.add_points(
+            unmatched_annotations,
+            properties={"annotator": unmatched_annotators},
+            face_color="annotator",
+            face_color_cycle=COLOR,  # TODO reorder
+        )
+        v.title = os.path.basename(fname)
+        napari.run()
+
+    else:
+        # Combine consensus and previous annotations.
+        consensus_annotations = consensus_annotations[["axis-0", "axis-1", "axis-2"]]
+        n_consensus = len(consensus_annotations)
+        n_unmatched = len(unmatched_annotations)
+        if prev_consensus is not None:
+            n_prev = len(prev_consensus)
+            print("Number of previous consensus annotations:", n_prev)
+
+        print("Number of new consensus annotations:", n_consensus)
+        print("Number of unmatched annotations:", n_unmatched)
+
+        consensus_annotations = pd.concat([consensus_annotations, prev_consensus])
+
+        out_name = fname.replace(".tif", ".csv")
+        out_path = os.path.join(OUTPUT_FOLDER, out_name)
+
+        os.makedirs(OUTPUT_FOLDER, exist_ok=True)
+        consensus_annotations.to_csv(out_path, index=False)
+
+
+def main():
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--images", nargs="+")
+    parser.add_argument("--check", action="store_true")
+    args = parser.parse_args()
+
+    if args.images is None:
+        image_paths = sorted(glob(os.path.join(ROOT, "*.tif")))
+    else:
+        image_paths = args.images
+
+    for image_path in image_paths:
+        create_consensus_step2(image_path, args.check)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/validation/SGNs/rescale_annotations.py b/scripts/validation/SGNs/rescale_annotations.py
@@ -1,13 +1,15 @@
+import json
 import os
 import shutil
 from glob import glob
 
+import imageio.v3 as imageio
 import numpy as np
 import pandas as pd
 import tifffile
 import zarr
 
-from flamingo_tools.s3_utils import get_s3_path, BUCKET_NAME, SERVICE_ENDPOINT
+from flamingo_tools.s3_utils import get_s3_path, BUCKET_NAME, SERVICE_ENDPOINT, create_s3_target
 
 
 def get_scale_factor():
@@ -28,54 +30,146 @@ def get_scale_factor():
     return scale_factor
 
 
-def rescale_annotations(input_path, scale_factor, bkp_folder):
-    annotations = pd.read_csv(input_path)
+def get_shape():
+    cochlea = "M_LR_000169_R_fused"
+    s3 = create_s3_target()
+    content = s3.open(f"{BUCKET_NAME}/{cochlea}/dataset.json", mode="r", encoding="utf-8")
+    info = json.loads(content.read())
+    print("Available sources:")
+    for source in info["sources"].keys():
+        print(source)
 
+    internal_path = os.path.join(cochlea, "images",  "ome-zarr", "PV.ome.zarr")
+    s3_store, fs = get_s3_path(internal_path, bucket_name=BUCKET_NAME, service_endpoint=SERVICE_ENDPOINT)
+
+    input_key = "s0"
+    with zarr.open(s3_store, mode="r") as f:
+        new_shape = f[input_key].shape
+    return new_shape
+
+
+def rescale_annotations(annotation_file, output_folder, new_shape, original_shape):
+    # 0.) Split the name into its parts.
+    fname = os.path.basename(annotation_file)
+    name_components = fname.split("_")
+    z = int(name_components[2][1:])
+
+    # 1.) Find the matching raw file and get its shape.
+    root = os.path.split(os.path.split(annotation_file)[0])[0]
+    tif_name = "_".join(name_components[:-1])
+    image_file = os.path.join(root, f"{tif_name}.tif")
+    assert os.path.exists(image_file), image_file
+    this_shape = tifffile.memmap(image_file).shape
+
+    # 2.) Determine if the annotations have to be reshaped,
+    if this_shape[1:] == new_shape[1:]:  # No, they don't have to be reshaped.
+        # In this case we copy the annotations and that's it.
+        print(annotation_file, "does not need to be rescaled")
+        output_path = os.path.join(output_folder, fname)
+        shutil.copyfile(annotation_file, output_path)
+        return
+    elif this_shape[1:] == original_shape[1:]:  # Yes, they have to be reshaped
+        pass
+    else:
+        raise ValueError(f"Unexpected shape: {this_shape}")
+
+    # 3.) Rescale the annotations.
+    scale_factor = [float(ns) / float(os) for ns, os in zip(new_shape, original_shape)]
+
+    annotations = pd.read_csv(annotation_file)
     annotations_rescaled = annotations.copy()
     annotations_rescaled["axis-1"] = annotations["axis-1"] * scale_factor[1]
     annotations_rescaled["axis-2"] = annotations["axis-2"] * scale_factor[2]
 
-    fname = os.path.basename(input_path)
-    name_components = fname.split("_")
-    z = int(name_components[2][1:])
     new_z = int(np.round(z * scale_factor[0]))
-
     name_components[2] = f"z{new_z}"
     name_components = name_components[:-1] + ["rescaled"] + name_components[-1:]
     new_fname = "_".join(name_components)
 
-    input_folder = os.path.split(input_path)[0]
-    out_path = os.path.join(input_folder, new_fname)
-    bkp_path = os.path.join(bkp_folder, fname)
+    output_path = os.path.join(output_folder, new_fname)
+    annotations_rescaled.to_csv(output_path, index=False)
 
-    # print(input_path)
-    # print(out_path)
-    # print(bkp_path)
-    # print()
-    # return
 
-    shutil.move(input_path, bkp_path)
-    annotations_rescaled.to_csv(out_path, index=False)
-
-
-def main():
-    # scale_factor = get_scale_factor()
-    # print(scale_factor)
-    scale_factor = (2.6314,) * 3
+def rescale_all_annotations():
+    prefix = "MLR169R_PV"
+    # shape = get_shape()
+    original_shape = (1921, 1479, 2157)
+    new_shape = (5089, 3915, 5665)
 
     root = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationSGNs"
     annotation_folders = ["AnnotationsEK", "AnnotationsAMD", "AnnotationsLR"]
     for folder in annotation_folders:
-        bkp_folder = os.path.join(root, folder, "rescaled_bkp")
-        os.makedirs(bkp_folder, exist_ok=True)
+        output_folder = os.path.join(root, folder, "rescaled")
+        os.makedirs(output_folder, exist_ok=True)
 
         files = glob(os.path.join(root, folder, "*.csv"))
         for annotation_file in files:
             fname = os.path.basename(annotation_file)
-            if not fname.startswith(("MLR169R_PV_z722", "MLR169R_PV_z979")):
+            if not fname.startswith(prefix):
                 continue
             print("Rescaling", annotation_file)
-            rescale_annotations(annotation_file, scale_factor, bkp_folder)
+            rescale_annotations(annotation_file, output_folder, new_shape, original_shape)
+
+
+# Download the two new slices.
+def download_new_data():
+    from flamingo_tools.validation import fetch_data_for_evaluation
+
+    root = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationSGNs"
+    output_folder = os.path.join(root, "for_consensus_annotation")
+    os.makedirs(output_folder, exist_ok=True)
+
+    cochlea = "M_LR_000169_R_fused"
+    files = [
+        "AnnotationsEK/rescaled/MLR169R_PV_z1913_base_full_rescaled_annotationsEK.csv",
+        "AnnotationsEK/rescaled/MLR169R_PV_z2594_mid_full_rescaled_annotationsEK.csv"
+    ]
+    for ff in files:
+        annotation_path = os.path.join(root, ff)
+
+        fname = os.path.basename(annotation_path)
+        name_components = fname.split("_")
+
+        tif_name = "_".join(name_components[:-1])
+        image_file = os.path.join(output_folder, f"{tif_name}.tif")
+
+        _, _, image = fetch_data_for_evaluation(
+            annotation_path, cache_path=None, cochlea=cochlea, extra_data="PV", z_extent=10
+        )
+        print(image.shape)
+        print("Writing to:", image_file)
+        imageio.imwrite(image_file, image)
+
+
+def check_rescaled_annotations():
+    import napari
+    from flamingo_tools.validation import fetch_data_for_evaluation
+
+    root = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationSGNs"
+    annotation_folders = ["AnnotationsEK/rescaled", "AnnotationsAMD/rescaled", "AnnotationsLR/rescaled"]
+    cochlea = "M_LR_000169_R_fused"
+
+    for folder in annotation_folders:
+        annotation_paths = sorted(glob(os.path.join(root, folder, "*.csv")))
+        for annotation_path in annotation_paths:
+            segmentation, annotations, image = fetch_data_for_evaluation(
+                annotation_path, cache_path=None, components_for_postprocessing=[1], cochlea=cochlea, extra_data="PV",
+            )
+            v = napari.Viewer()
+            v.add_image(image)
+            v.add_labels(segmentation)
+            v.add_points(annotations)
+            v.title = annotation_path
+            napari.run()
+
+
+def main():
+    # rescale_all_annotations()
+    # check_rescaled_annotations()
+
+    # MLR169R_PV_z1913_base_full_rescaled.tif
+    # MLR169R_PV_z2594_mid_full_rescaled.tif
+    download_new_data()
 
 
 # Rescale the point annotations for the cochlea MLR169R, which was
diff --git a/scripts/validation/SGNs/run_evaluation.py b/scripts/validation/SGNs/run_evaluation.py
@@ -3,7 +3,7 @@
 
 import pandas as pd
 from flamingo_tools.validation import (
-    fetch_data_for_evaluation, parse_annotation_path, compute_scores_for_annotated_slice
+    fetch_data_for_evaluation, _parse_annotation_path, compute_scores_for_annotated_slice
 )
 
 ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationSGNs"
@@ -27,11 +27,16 @@ def run_evaluation(root, annotation_folders, result_file, cache_folder):
         annotator = folder[len("Annotations"):]
         annotations = sorted(glob(os.path.join(root, folder, "*.csv")))
         for annotation_path in annotations:
-            cochlea, slice_id = parse_annotation_path(annotation_path)
+            cochlea, slice_id = _parse_annotation_path(annotation_path)
 
             print("Run evaluation for", annotator, cochlea, "z=", slice_id)
+            if cochlea == "M_LR_000169_R":
+                mobie_name = f"{cochlea}_fused"
+            else:
+                mobie_name = None
+
             segmentation, annotations = fetch_data_for_evaluation(
-                annotation_path, components_for_postprocessing=[1],
+                annotation_path, components_for_postprocessing=[1], cochlea=mobie_name,
                 cache_path=None if cache_folder is None else os.path.join(cache_folder, f"{cochlea}_{slice_id}.tif")
             )
             scores = compute_scores_for_annotated_slice(segmentation, annotations, matching_tolerance=5)
@@ -51,7 +56,7 @@ def main():
     import argparse
     parser = argparse.ArgumentParser()
     parser.add_argument("-i", "--input", default=ROOT)
-    parser.add_argument("--folders", default=ANNOTATION_FOLDERS)
+    parser.add_argument("--folders", default=ANNOTATION_FOLDERS, nargs="+")
     parser.add_argument("--result_file", default="results.csv")
     parser.add_argument("--cache_folder")
     args = parser.parse_args()
