Update evaluation for IHCs

Author: constantinpape

flamingo_tools/validation.py

@@ -43,6 +43,8 @@ def fetch_data_for_evaluation(
     seg_name: str = "SGN_v2",
     z_extent: int = 0,
     components_for_postprocessing: Optional[List[int]] = None,
+    cochlea: Optional[str] = None,
+    extra_data: Optional[str] = None,
 ) -> Tuple[np.ndarray, pd.DataFrame]:
     """Fetch segmentation from S3 matching the annotation path for evaluation.
 
@@ -53,28 +55,31 @@ def fetch_data_for_evaluation(
         z_extent: Additional z-slices to load from the segmentation.
         components_for_postprocessing: The component ids for restricting the segmentation to.
             Choose [1] for the default componentn containing the helix.
+        cochlea: Optional name of the cochlea.
+        extra_data: Extra data to fetch.
 
     Returns:
         The segmentation downloaded from the S3 bucket.
         The annotations loaded from pandas and matching the segmentation.
     """
     # Load the annotations and normalize them for the given z-extent.
     annotations = pd.read_csv(annotation_path)
-    annotations = annotations.drop(columns="index")
+    if "index" in annotations.columns:
+        annotations = annotations.drop(columns="index")
     if z_extent == 0:  # If we don't have a z-extent then we just drop the first axis and rename the other two.
         annotations = annotations.drop(columns="axis-0")
         annotations = annotations.rename(columns={"axis-1": "axis-0", "axis-2": "axis-1"})
-    else:  # Otherwise we have to center the first axis.
-        # TODO
-        raise NotImplementedError
 
     # Load the segmentaiton from cache path if it is given and if it is already cached.
     if cache_path is not None and os.path.exists(cache_path):
         segmentation = imageio.imread(cache_path)
         return segmentation, annotations
 
     # Parse which ID and which cochlea from the name.
-    cochlea, slice_id = _parse_annotation_path(annotation_path)
+    if cochlea is None:
+        cochlea, slice_id = _parse_annotation_path(annotation_path)
+    else:
+        _, slice_id = _parse_annotation_path(annotation_path)
 
     # Open the S3 connection, get the path to the SGN segmentation in S3.
     internal_path = os.path.join(cochlea, "images",  "ome-zarr", f"{seg_name}.ome.zarr")
@@ -111,6 +116,14 @@ def fetch_data_for_evaluation(
     if cache_path is not None:
         imageio.imwrite(cache_path, segmentation, compression="zlib")
 
+    if extra_data is not None:
+        internal_path = os.path.join(cochlea, "images",  "ome-zarr", f"{extra_data}.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:
+            extra_im_data = f[input_key][roi]
+        return segmentation, annotations, extra_im_data
+
     return segmentation, annotations
 
 
@@ -347,6 +360,62 @@ def union(a, b):
     return consensus_df, unmatched_df
 
 
+def match_detections(
+    detections: np.ndarray,
+    annotations: np.ndarray,
+    max_dist: float
+):
+    """One-to-one matching between 3-D detections and ground-truth points.
+
+    Args:
+        detections: N x 3 candidate detections.
+        annotations: M x 3 ground-truth annotations for the reference points.
+        max_dist: Maximum Euclidean distance allowed for a match.
+
+    Returns:
+        Indices in `detections` that were matched (true positives).
+        Indices in `annotations` that were matched (true positives).
+        Unmatched detection indices (false positives).
+        Unmatched annotation indices (false negatives).
+    """
+    det = np.asarray(detections, dtype=float)
+    ann = np.asarray(annotations, dtype=float)
+    N, M = len(det), len(ann)
+
+    # trivial corner cases --------------------------------------------------------
+    if N == 0:
+        return np.empty(0, int), np.empty(0, int), np.empty(0, int), np.arange(M)
+    if M == 0:
+        return np.empty(0, int), np.empty(0, int), np.arange(N), np.empty(0, int)
+
+    # 1. build sparse radius-filtered distance matrix -----------------------------
+    tree_det = cKDTree(det)
+    tree_ann = cKDTree(ann)
+    coo = tree_det.sparse_distance_matrix(tree_ann, max_dist, output_type="coo_matrix")
+
+    if coo.nnz == 0:                       # nothing is close enough
+        return np.empty(0, int), np.empty(0, int), np.arange(N), np.arange(M)
+
+    cost = np.full((N, M), 5 * max_dist, dtype=float)
+    cost[coo.row, coo.col] = coo.data      # fill only existing edges
+
+    # 2. optimal one-to-one assignment (Hungarian) --------------------------------
+    row_ind, col_ind = linear_sum_assignment(cost)
+
+    # Filter assignments that were padded with +∞ cost for non-existent edges
+    # (linear_sum_assignment automatically does that padding internally).
+    valid_mask = cost[row_ind, col_ind] <= max_dist
+    tp_det_ids = row_ind[valid_mask]
+    tp_ann_ids = col_ind[valid_mask]
+    assert len(tp_det_ids) == len(tp_ann_ids)
+
+    # 3. derive FP / FN -----------------------------------------------------------
+    fp_det_ids = np.setdiff1d(np.arange(N), tp_det_ids, assume_unique=True)
+    fn_ann_ids = np.setdiff1d(np.arange(M), tp_ann_ids, assume_unique=True)
+
+    return tp_det_ids, tp_ann_ids, fp_det_ids, fn_ann_ids
+
+
 def for_visualization(segmentation, annotations, matches):
     green_red = ["#00FF00", "#FF0000"]
 

scripts/validation/IHCs/consensus_annotations.py

@@ -0,0 +1,91 @@
+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/F1ValidationIHCs"
+ANNOTATION_FOLDERS = ["Annotations_AMD", "Annotations_EK", "Annotations_LR"]
+
+OUTPUT_FOLDER = os.path.join(ROOT, "consensus_annotation")
+COLOR = ["blue", "yellow", "orange"]
+
+
+def match_annotations(image_path):
+    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)]
+        assert len(matches) == 1
+        annotation_path = matches[0]
+        annotations[annotation_folder] = annotation_path
+
+    return annotations
+
+
+def consensus_annotations(image_path, check):
+    annotation_paths = match_annotations(image_path)
+    assert len(annotation_paths) == len(ANNOTATION_FOLDERS)
+
+    # I tried first with a matching distnce of 8, but that is too conservative.
+    # A mathing distance of 16 seems better, but might still need to refine this.
+    matching_distance = 16.0
+    consensus_annotations, unmatched_annotations = create_consensus_annotations(
+        annotation_paths, matching_distance=matching_distance, min_matches_for_consensus=2,
+    )
+    fname = os.path.basename(image_path)
+
+    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)
+        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:
+        os.makedirs(OUTPUT_FOLDER, exist_ok=True)
+        consensus_annotations = consensus_annotations[["axis-0", "axis-1", "axis-2"]]
+        consensus_annotations.insert(0, "annotator", ["consensus"] * len(consensus_annotations))
+        unmatched_annotations = unmatched_annotations[["axis-0", "axis-1", "axis-2", "annotator"]]
+        annotations = pd.concat([consensus_annotations, unmatched_annotations])
+        output_path = os.path.join(OUTPUT_FOLDER, fname.replace(".tif", ".csv"))
+        annotations.to_csv(output_path, index=False)
+        print("Saved to", output_path)
+
+
+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:
+        consensus_annotations(image_path, args.check)
+
+
+if __name__ == "__main__":
+    main()

scripts/validation/IHCs/consensus_evaluation.csv

@@ -0,0 +1,37 @@
+annotator,tps,fps,fns
+Annotations_AMD,24,0,4
+Annotations_EK,27,1,1
+Annotations_LR,24,6,4
+Annotations_AMD,5,0,0
+Annotations_EK,4,0,1
+Annotations_LR,4,2,1
+Annotations_AMD,31,1,1
+Annotations_EK,29,2,3
+Annotations_LR,31,1,1
+Annotations_AMD,5,0,2
+Annotations_EK,6,0,1
+Annotations_LR,5,2,2
+Annotations_AMD,26,0,1
+Annotations_EK,27,0,0
+Annotations_LR,27,4,0
+Annotations_AMD,31,0,0
+Annotations_EK,30,0,1
+Annotations_LR,31,2,0
+Annotations_AMD,7,0,2
+Annotations_EK,8,0,1
+Annotations_LR,8,1,1
+Annotations_AMD,28,1,2
+Annotations_EK,30,1,0
+Annotations_LR,27,6,3
+Annotations_AMD,43,0,0
+Annotations_EK,42,0,1
+Annotations_LR,43,4,0
+Annotations_AMD,34,0,2
+Annotations_EK,35,0,1
+Annotations_LR,36,1,0
+Annotations_AMD,5,0,0
+Annotations_EK,5,0,0
+Annotations_LR,5,0,0
+Annotations_AMD,31,0,2
+Annotations_EK,33,0,0
+Annotations_LR,33,1,0

scripts/validation/IHCs/evaluate_consensus.py

@@ -0,0 +1,56 @@
+import os
+from glob import glob
+
+import pandas as pd
+from flamingo_tools.validation import match_detections
+
+ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationIHCs"
+ANNOTATION_FOLDERS = ["Annotations_AMD", "Annotations_EK", "Annotations_LR"]
+CONSENSUS_FOLDER = "consensus_annotation"
+
+
+def match_annotations(consensus_path):
+    annotations = {}
+    prefix = os.path.basename(consensus_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)]
+        assert len(matches) == 1
+        annotation_path = matches[0]
+        annotations[annotation_folder] = annotation_path
+
+    return annotations
+
+
+def main():
+    consensus_files = sorted(glob(os.path.join(ROOT, CONSENSUS_FOLDER, "*.csv")))
+
+    results = {
+        "annotator": [],
+        "tps": [],
+        "fps": [],
+        "fns": [],
+    }
+    for consensus_file in consensus_files:
+        consensus = pd.read_csv(consensus_file)
+        consensus = consensus[consensus.annotator == "consensus"][["axis-0", "axis-1", "axis-2"]]
+
+        annotations = match_annotations(consensus_file)
+        for name, annotation_path in annotations.items():
+            annotation = pd.read_csv(annotation_path)[["axis-0", "axis-1", "axis-2"]]
+            tp, _, fp, fn = match_detections(annotation, consensus, max_dist=12.0)
+            results["annotator"].append(name)
+            results["tps"].append(len(tp))
+            results["fps"].append(len(fp))
+            results["fns"].append(len(fn))
+
+    results = pd.DataFrame(results)
+    print(results)
+    results.to_csv("consensus_evaluation.csv", index=False)
+
+
+if __name__ == "__main__":
+    main()

scripts/validation/IHCs/run_evaluation.py

@@ -8,7 +8,8 @@
 
 ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/AnnotatedImageCrops/F1ValidationIHCs"
 # ANNOTATION_FOLDERS = ["AnnotationsEK", "AnnotationsAMD", "AnnotationsLR"]
-ANNOTATION_FOLDERS = ["Annotations_AMD", "Annotations_LR"]
+# ANNOTATION_FOLDERS = ["Annotations_AMD", "Annotations_LR"]
+ANNOTATION_FOLDERS = ["consensus_annotation"]
 
 
 def run_evaluation(root, annotation_folders, result_file, cache_folder):
@@ -25,7 +26,7 @@ def run_evaluation(root, annotation_folders, result_file, cache_folder):
         os.makedirs(cache_folder, exist_ok=True)
 
     for folder in annotation_folders:
-        annotator = folder[len("Annotations"):]
+        annotator = "consensus" if folder == "consensus_annotation" else folder[len("Annotations"):]
         annotations = sorted(glob(os.path.join(root, folder, "*.csv")))
         for annotation_path in annotations:
             print(annotation_path)