intra_annotator_agreement.py

#!/usr/bin/env python3
"""Intra-annotator agreement analysis for MISGRA annotations.

Reads Intra_annotator.xlsx to pair findings across two annotation sessions
by the same annotator. Computes per-pair clinical metrics (Euclidean distance,
Dice, volume/length/angle differences) and aggregate ICC(3,1) statistics.

Multi-annotation findings are matched via Hungarian algorithm.
Zero-size ROIs are excluded; mixed-type pairs are reported but not compared.
"""

import argparse
import csv
import itertools
import json
import math
import os
import sys
from collections import defaultdict

import numpy as np
import openpyxl


# ── Annotator directory mapping ──────────────────────────────────────────────
# Maps spreadsheet annotator name -> (directory name, annotations subfolder)
ANNOTATOR_MAP = {
    "Annotator_1": ("Annotator_1", "A1_annotations"),
    "Annotator_2": ("Annotator_2", "A2_annotations"),
    "Annotator_3": ("Annotator_3", "A3_annotations"),
    "Annotator_4": ("Annotator_4", "A4_annotations"),
    "Annotator_5": ("Annotator_5", "A5_annotations"),
}

OUTPUT_COLS = [
    "annotator", "path1", "path2",
    "annotation1", "annotation2",
    "type1", "type2",
    "status", "match_method", "match_group", "match_cost",
    "description1", "description2",
    # Point metrics
    "euclidean_distance_mm",
    # ROI metrics
    "centre_distance_mm", "volume_1_mm3", "volume_2_mm3",
    "relative_volume_diff_pct", "dice_coefficient",
    # Line metrics
    "length_1_mm", "length_2_mm", "length_diff_mm", "angle_diff_deg",
]


# ── ICC computation ─────────────────────────────────────────────────────────

def calculate_icc_2rater(measurements1, measurements2):
    """ICC(3,1) for two raters — two-way mixed, consistency, single measures.

    Returns ICC value, or None if < 2 cases.
    """
    m1 = np.asarray(measurements1, dtype=np.float64)
    m2 = np.asarray(measurements2, dtype=np.float64)
    n = len(m1)
    if n < 2:
        return None
    k = 2

    data = np.column_stack([m1, m2])
    grand_mean = np.mean(data)
    row_means = np.mean(data, axis=1)
    col_means = np.mean(data, axis=0)

    ss_cases = k * np.sum((row_means - grand_mean) ** 2)
    ss_raters = n * np.sum((col_means - grand_mean) ** 2)
    ss_total = np.sum((data - grand_mean) ** 2)
    ss_error = ss_total - ss_cases - ss_raters

    ms_cases = ss_cases / (n - 1)
    ms_error = ss_error / ((n - 1) * (k - 1))

    denom = ms_cases + (k - 1) * ms_error
    if denom == 0:
        return None
    return float((ms_cases - ms_error) / denom)


# ── JSON loaders ─────────────────────────────────────────────────────────────

def annotation_type_from_filename(filename):
    """Extract annotation type from filename prefix (e.g. 'point_1.json' -> 'point')."""
    return filename.split("_")[0]


def load_annotation(json_path):
    """Load a 3D Slicer Markups JSON file and extract type-specific fields.

    Supports Fiducial (point), ROI (bounding box), and Line markups.
    Returns None if file missing, ROI is zero-size, or type unrecognised.
    """
    if not os.path.isfile(json_path):
        return None
    with open(json_path) as f:
        data = json.load(f)
    markup = data["markups"][0]
    mtype = markup["type"]

    if mtype == "Fiducial":
        cp = markup["controlPoints"][0]
        return {
            "type": "point",
            "position": np.array(cp["position"], dtype=np.float64),
        }
    elif mtype == "ROI":
        size = markup.get("size", [0, 0, 0])
        if all(s == 0 for s in size):
            return None
        centre = np.array(markup["center"], dtype=np.float64)
        size_arr = np.array(size, dtype=np.float64)
        return {
            "type": "roi",
            "centre": centre,
            "size": size_arr,
            "volume": float(np.prod(size_arr)),
        }
    elif mtype == "Line":
        cps = markup["controlPoints"]
        p1 = np.array(cps[0]["position"], dtype=np.float64)
        p2 = np.array(cps[1]["position"], dtype=np.float64)
        # Prefer stored measurement length, fall back to Euclidean distance
        length = None
        for m in markup.get("measurements", []):
            if m["name"] == "length" and m.get("enabled", False):
                length = m["value"]
        if length is None:
            length = float(np.linalg.norm(p2 - p1))
        direction = p2 - p1
        norm_val = np.linalg.norm(direction)
        if norm_val > 0:
            direction = direction / norm_val
        return {
            "type": "line",
            "p1": p1, "p2": p2,
            "length": length,
            "direction": direction,
        }
    return None


# ── Comparison functions ─────────────────────────────────────────────────────

def compare_points(a1, a2):
    """Euclidean distance (mm) between two point annotations."""
    dist = float(np.linalg.norm(a1["position"] - a2["position"]))
    return {"euclidean_distance_mm": dist}


def dice_coefficient_3d(centre1, size1, centre2, size2):
    """3D Dice coefficient between two axis-aligned bounding boxes (0-1)."""
    min1 = centre1 - size1 / 2.0
    max1 = centre1 + size1 / 2.0
    min2 = centre2 - size2 / 2.0
    max2 = centre2 + size2 / 2.0

    inter_min = np.maximum(min1, min2)
    inter_max = np.minimum(max1, max2)
    inter_size = np.maximum(inter_max - inter_min, 0.0)
    intersection = float(np.prod(inter_size))

    vol1 = float(np.prod(size1))
    vol2 = float(np.prod(size2))
    if (vol1 + vol2) == 0:
        return 0.0
    return (2.0 * intersection) / (vol1 + vol2)


def compare_rois(a1, a2):
    """Compare two ROIs by centre distance, relative volume difference, and Dice."""
    centre_dist = float(np.linalg.norm(a1["centre"] - a2["centre"]))
    v1, v2 = a1["volume"], a2["volume"]
    mean_vol = (v1 + v2) / 2.0
    rel_vol_diff = abs(v1 - v2) / mean_vol * 100.0 if mean_vol > 0 else ""
    dice = dice_coefficient_3d(a1["centre"], a1["size"], a2["centre"], a2["size"])
    return {
        "centre_distance_mm": centre_dist,
        "volume_1_mm3": v1,
        "volume_2_mm3": v2,
        "relative_volume_diff_pct": rel_vol_diff,
        "dice_coefficient": dice,
    }


def angle_between_lines(d1, d2):
    """Angle (0-90 deg) between two line directions. Uses abs(dot) so antiparallel = 0."""
    cos_angle = abs(float(np.dot(d1, d2)))
    cos_angle = min(cos_angle, 1.0)
    return math.degrees(math.acos(cos_angle))


def compare_lines(a1, a2):
    """Compare two lines by length difference (mm) and angular difference (deg)."""
    length_diff = abs(a1["length"] - a2["length"])
    angle_diff = angle_between_lines(a1["direction"], a2["direction"])
    return {
        "length_1_mm": a1["length"],
        "length_2_mm": a2["length"],
        "length_diff_mm": length_diff,
        "angle_diff_deg": angle_diff,
    }


# ── Hungarian matching ──────────────────────────────────────────────────────
# Optimal assignment for comma-separated multi-annotation findings.

def compute_cost(a1, a2, atype):
    """Matching cost between two annotations (Euclidean dist for points/ROIs,
    length diff + angle diff for lines)."""
    if atype == "point":
        return float(np.linalg.norm(a1["position"] - a2["position"]))
    elif atype == "roi":
        return float(np.linalg.norm(a1["centre"] - a2["centre"]))
    elif atype == "line":
        length_cost = abs(a1["length"] - a2["length"])
        angle_cost = angle_between_lines(a1["direction"], a2["direction"])
        return length_cost + angle_cost
    return float("inf")


def hungarian_match(list1, list2, atype):
    """Optimal one-to-one matching between two annotation lists.

    Brute-force permutation search (n is typically 2-5).
    Returns (matched, unmatched_1, unmatched_2).
    """
    n1, n2 = len(list1), len(list2)
    if n1 == 0 or n2 == 0:
        return [], list(range(n1)), list(range(n2))

    cost_matrix = np.full((n1, n2), float("inf"))
    for i in range(n1):
        for j in range(n2):
            if list1[i] is not None and list2[j] is not None:
                cost_matrix[i, j] = compute_cost(list1[i], list2[j], atype)

    n_match = min(n1, n2)
    best_cost = float("inf")
    best_assignment = None

    if n1 <= n2:
        for perm in itertools.permutations(range(n2), n_match):
            total = sum(cost_matrix[i, perm[i]] for i in range(n_match))
            if total < best_cost:
                best_cost = total
                best_assignment = [(i, perm[i]) for i in range(n_match)]
    else:
        for perm in itertools.permutations(range(n1), n_match):
            total = sum(cost_matrix[perm[j], j] for j in range(n_match))
            if total < best_cost:
                best_cost = total
                best_assignment = [(perm[j], j) for j in range(n_match)]

    matched = [(i, j, cost_matrix[i, j]) for i, j in best_assignment]
    matched_i = {i for i, _, _ in matched}
    matched_j = {j for _, j, _ in matched}
    return (matched,
            [i for i in range(n1) if i not in matched_i],
            [j for j in range(n2) if j not in matched_j])


# ── Main ─────────────────────────────────────────────────────────────────────

def main():
    parser = argparse.ArgumentParser(description=__doc__,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument(
        "--spreadsheet", required=True,
        help="Path to Intra_annotator.xlsx spreadsheet defining matched findings",
    )
    parser.add_argument(
        "--annotations-root", required=True,
        help="Root folder containing annotator directories with annotation data",
    )
    parser.add_argument(
        "--output-csv",
        default="intra_annotator_results.csv",
    )
    args = parser.parse_args()

    ann_root = args.annotations_root
    if not os.path.isdir(ann_root):
        print(f"ERROR: annotations root not found: {ann_root}", file=sys.stderr)
        sys.exit(1)

    wb = openpyxl.load_workbook(args.spreadsheet)
    ws = wb.active
    rows_data = list(ws.iter_rows(min_row=2, values_only=True))

    results = []

    # Accumulate paired measurements for ICC computation
    icc_points_data = []   # (annotator, position1, position2)
    icc_rois_data = []     # (annotator, roi_dict_1, roi_dict_2)
    icc_lines_data = []    # (annotator, length1, length2)

    # Track annotation counts per session pair for finding-level agreement
    finding_counts = defaultdict(lambda: [0, 0])

    for row_idx, row in enumerate(rows_data, start=2):
        annotator_raw, _, path1_raw, ann1_raw, desc1_raw, path2_raw, ann2_raw, desc2_raw = row

        if annotator_raw is None:
            continue
        annotator = str(annotator_raw).strip()
        if annotator not in ANNOTATOR_MAP:
            print(f"WARNING R{row_idx}: Unknown annotator '{annotator}', skipping",
                  file=sys.stderr)
            continue

        full_name, sub = ANNOTATOR_MAP[annotator]
        base_dir = os.path.join(ann_root, full_name, sub)

        path1 = str(path1_raw).strip() if path1_raw else ""
        path2 = str(path2_raw).strip() if path2_raw else ""
        ann1_str = str(ann1_raw).strip() if ann1_raw else ""
        ann2_str = str(ann2_raw).strip() if ann2_raw else ""
        desc1 = str(desc1_raw).strip() if desc1_raw else ""
        desc2 = str(desc2_raw).strip() if desc2_raw else ""

        session_key = (annotator, path1, path2)
        if ann1_str:
            finding_counts[session_key][0] += len(ann1_str.split(","))
        if ann2_str:
            finding_counts[session_key][1] += len(ann2_str.split(","))

        # ── Case 1: Unmatched findings ──
        if not ann1_str and not ann2_str:
            continue
        if not ann1_str or not ann2_str:
            present_side = "1" if ann1_str else "2"
            present_anns = ann1_str if ann1_str else ann2_str
            present_desc = desc1 if ann1_str else desc2
            for af in present_anns.split(","):
                af = af.strip()
                r = {c: "" for c in OUTPUT_COLS}
                r["annotator"] = annotator
                r["path1"] = path1
                r["path2"] = path2
                r["status"] = "unmatched"
                r["match_method"] = "unmatched"
                if present_side == "1":
                    r["annotation1"] = af
                    r["type1"] = annotation_type_from_filename(af)
                    r["description1"] = present_desc
                else:
                    r["annotation2"] = af
                    r["type2"] = annotation_type_from_filename(af)
                    r["description2"] = present_desc
                results.append(r)
            continue

        ann1_files = [a.strip() for a in ann1_str.split(",")]
        ann2_files = [a.strip() for a in ann2_str.split(",")]

        type1_set = set(annotation_type_from_filename(f) for f in ann1_files)
        type2_set = set(annotation_type_from_filename(f) for f in ann2_files)

        # ── Case 2: Mixed-type pair ──
        if len(ann1_files) == 1 and len(ann2_files) == 1 and type1_set != type2_set:
            r = {c: "" for c in OUTPUT_COLS}
            r["annotator"] = annotator
            r["path1"] = path1
            r["path2"] = path2
            r["annotation1"] = ann1_files[0]
            r["annotation2"] = ann2_files[0]
            r["type1"] = annotation_type_from_filename(ann1_files[0])
            r["type2"] = annotation_type_from_filename(ann2_files[0])
            r["status"] = "mixed_type"
            r["match_method"] = "mixed_type"
            r["description1"] = desc1
            r["description2"] = desc2
            results.append(r)
            continue

        def process_match(af1, af2, atype, method, group_id="", cost=""):
            file1 = os.path.join(base_dir, path1.lstrip("./"), af1)
            file2 = os.path.join(base_dir, path2.lstrip("./"), af2)
            a1 = load_annotation(file1)
            a2 = load_annotation(file2)

            r = {c: "" for c in OUTPUT_COLS}
            r["annotator"] = annotator
            r["path1"] = path1
            r["path2"] = path2
            r["annotation1"] = af1
            r["annotation2"] = af2
            r["type1"] = atype
            r["type2"] = atype
            r["match_method"] = method
            r["match_group"] = group_id
            r["match_cost"] = cost
            r["description1"] = desc1
            r["description2"] = desc2

            if a1 is None or a2 is None:
                r["status"] = "excluded_zero_roi" if atype == "roi" else "file_error"
                results.append(r)
                return

            r["status"] = "compared"
            if atype == "point":
                r.update(compare_points(a1, a2))
                icc_points_data.append((annotator, a1["position"], a2["position"]))
            elif atype == "roi":
                r.update(compare_rois(a1, a2))
                icc_rois_data.append((annotator, a1, a2))
            elif atype == "line":
                r.update(compare_lines(a1, a2))
                icc_lines_data.append((annotator, a1["length"], a2["length"]))
            results.append(r)

        # ── Case 3: Single 1:1 same-type pair ──
        if len(ann1_files) == 1 and len(ann2_files) == 1:
            atype = annotation_type_from_filename(ann1_files[0])
            process_match(ann1_files[0], ann2_files[0], atype, "direct")
            continue

        # ── Case 4: Multi-annotation matching ──
        group_id = f"grp_R{row_idx}"
        all_types = type1_set | type2_set
        if len(all_types) > 1:
            for af in ann1_files:
                r = {c: "" for c in OUTPUT_COLS}
                r["annotator"] = annotator
                r["path1"] = path1
                r["path2"] = path2
                r["annotation1"] = af
                r["type1"] = annotation_type_from_filename(af)
                r["status"] = "mixed_type"
                r["match_method"] = "mixed_type"
                r["match_group"] = group_id
                r["description1"] = desc1
                results.append(r)
            continue

        atype = all_types.pop()
        loaded1 = [load_annotation(os.path.join(base_dir, path1.lstrip("./"), af))
                    for af in ann1_files]
        loaded2 = [load_annotation(os.path.join(base_dir, path2.lstrip("./"), af))
                    for af in ann2_files]

        valid1 = [(i, a) for i, a in enumerate(loaded1) if a is not None]
        valid2 = [(i, a) for i, a in enumerate(loaded2) if a is not None]

        if not valid1 or not valid2:
            for af in ann1_files:
                r = {c: "" for c in OUTPUT_COLS}
                r["annotator"] = annotator
                r["path1"] = path1
                r["path2"] = path2
                r["annotation1"] = af
                r["type1"] = atype
                r["status"] = "excluded_zero_roi" if atype == "roi" else "file_error"
                r["match_method"] = "auto_unmatched"
                r["match_group"] = group_id
                r["description1"] = desc1
                results.append(r)
            continue

        anns1 = [a for _, a in valid1]
        anns2 = [a for _, a in valid2]
        matched, unmatched_1, unmatched_2 = hungarian_match(anns1, anns2, atype)

        for i, j, cost_val in matched:
            orig_i = valid1[i][0]
            orig_j = valid2[j][0]
            process_match(ann1_files[orig_i], ann2_files[orig_j], atype,
                          "auto_matched", group_id, f"{cost_val:.4f}")

        for i in unmatched_1:
            orig_i = valid1[i][0]
            r = {c: "" for c in OUTPUT_COLS}
            r["annotator"] = annotator
            r["path1"] = path1
            r["path2"] = path2
            r["annotation1"] = ann1_files[orig_i]
            r["type1"] = atype
            r["status"] = "auto_unmatched"
            r["match_method"] = "auto_unmatched"
            r["match_group"] = group_id
            r["description1"] = desc1
            results.append(r)

        for j in unmatched_2:
            orig_j = valid2[j][0]
            r = {c: "" for c in OUTPUT_COLS}
            r["annotator"] = annotator
            r["path1"] = path1
            r["path2"] = path2
            r["annotation2"] = ann2_files[orig_j]
            r["type2"] = atype
            r["status"] = "auto_unmatched"
            r["match_method"] = "auto_unmatched"
            r["match_group"] = group_id
            r["description2"] = desc2
            results.append(r)

    # ── Write output CSV ─────────────────────────────────────────────────
    with open(args.output_csv, "w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=OUTPUT_COLS)
        writer.writeheader()
        for r in results:
            writer.writerow(r)
    print(f"Detailed results written to: {args.output_csv}\n")

    # ── Summary statistics ───────────────────────────────────────────────
    # Filter to only successfully compared pairs for aggregate reporting
    compared = [r for r in results if r["status"] == "compared"]

    print("=" * 70)
    print("INTRA-ANNOTATOR AGREEMENT SUMMARY")
    print("=" * 70)

    # ── Points summary ──────────────────────────────────────────────────
    # Report descriptive statistics of Euclidean distances across all matched
    # point pairs, then ICC(3,1) per spatial axis (L/P/S), both overall and
    # per annotator.
    points = [r for r in compared if r["type1"] == "point"]
    if points:
        dists = np.array([r["euclidean_distance_mm"] for r in points])
        print(f"\n--- Points (n={len(points)}) ---")
        print(f"  Euclidean distance (mm):  mean={dists.mean():.2f}  "
              f"std={dists.std():.2f}  median={np.median(dists):.2f}  "
              f"min={dists.min():.2f}  max={dists.max():.2f}")

        if icc_points_data:
            pos1_all = np.array([p[1] for p in icc_points_data])
            pos2_all = np.array([p[2] for p in icc_points_data])
            for axis, label in enumerate(["L (left-right)", "P (post-ant)", "S (sup-inf)"]):
                icc_val = calculate_icc_2rater(pos1_all[:, axis], pos2_all[:, axis])
                icc_str = f"{icc_val:.4f}" if icc_val is not None else "n/a"
                print(f"  ICC(3,1) {label}-axis: {icc_str}")

        print("  Per annotator:")
        for ann in sorted(set(r["annotator"] for r in points)):
            d = np.array([r["euclidean_distance_mm"]
                          for r in points if r["annotator"] == ann])
            print(f"    {ann:12s}  n={len(d):3d}  "
                  f"mean={d.mean():.2f}  std={d.std():.2f}")

            ann_icc = [(p[1], p[2]) for p in icc_points_data if p[0] == ann]
            if len(ann_icc) >= 2:
                p1a = np.array([x[0] for x in ann_icc])
                p2a = np.array([x[1] for x in ann_icc])
                iccs = []
                for axis in range(3):
                    v = calculate_icc_2rater(p1a[:, axis], p2a[:, axis])
                    iccs.append(f"{v:.3f}" if v is not None else "n/a")
                print(f"    {'':12s}        ICC(L,P,S): {', '.join(iccs)}")
    else:
        print("\n--- Points: no matched pairs ---")

    # ── ROIs summary ──
    rois = [r for r in compared if r["type1"] == "roi"]
    if rois:
        cd = np.array([r["centre_distance_mm"] for r in rois])
        vd = np.array([r["relative_volume_diff_pct"]
                       for r in rois if r["relative_volume_diff_pct"] != ""])
        dice_vals = np.array([r["dice_coefficient"]
                              for r in rois if r["dice_coefficient"] != ""])
        print(f"\n--- ROIs (n={len(rois)}) ---")
        print(f"  Centre distance (mm):       mean={cd.mean():.2f}  "
              f"std={cd.std():.2f}  median={np.median(cd):.2f}  "
              f"min={cd.min():.2f}  max={cd.max():.2f}")
        if len(vd) > 0:
            print(f"  Relative volume diff (%):   mean={vd.mean():.2f}  "
                  f"std={vd.std():.2f}  median={np.median(vd):.2f}  "
                  f"min={vd.min():.2f}  max={vd.max():.2f}")
        if len(dice_vals) > 0:
            print(f"  Dice coefficient:           mean={dice_vals.mean():.4f}  "
                  f"std={dice_vals.std():.4f}  median={np.median(dice_vals):.4f}  "
                  f"min={dice_vals.min():.4f}  max={dice_vals.max():.4f}")

        if icc_rois_data:
            vol1 = np.array([d[1]["volume"] for d in icc_rois_data])
            vol2 = np.array([d[2]["volume"] for d in icc_rois_data])
            icc_vol = calculate_icc_2rater(vol1, vol2)
            print(f"  ICC(3,1) volume: {icc_vol:.4f}" if icc_vol is not None
                  else "  ICC(3,1) volume: n/a")

            size1 = np.array([d[1]["size"] for d in icc_rois_data])
            size2 = np.array([d[2]["size"] for d in icc_rois_data])
            for axis, label in enumerate(["LR (left-right)", "PA (post-ant)", "SI (sup-inf)"]):
                icc_dim = calculate_icc_2rater(size1[:, axis], size2[:, axis])
                icc_str = f"{icc_dim:.4f}" if icc_dim is not None else "n/a"
                print(f"  ICC(3,1) {label}: {icc_str}")

        print("  Per annotator:")
        for ann in sorted(set(r["annotator"] for r in rois)):
            c = np.array([r["centre_distance_mm"]
                          for r in rois if r["annotator"] == ann])
            v = np.array([r["relative_volume_diff_pct"]
                          for r in rois
                          if r["annotator"] == ann and r["relative_volume_diff_pct"] != ""])
            dc = np.array([r["dice_coefficient"]
                           for r in rois
                           if r["annotator"] == ann and r["dice_coefficient"] != ""])
            vstr = f"vol_diff: mean={v.mean():.2f}%" if len(v) > 0 else "vol_diff: n/a"
            dstr = f"dice: mean={dc.mean():.3f}" if len(dc) > 0 else "dice: n/a"
            print(f"    {ann:12s}  n={len(c):3d}  "
                  f"centre: mean={c.mean():.2f}  {vstr}  {dstr}")
    else:
        print("\n--- ROIs: no matched pairs ---")

    # ── Lines summary ──
    lines = [r for r in compared if r["type1"] == "line"]
    if lines:
        ld = np.array([r["length_diff_mm"] for r in lines])
        ad = np.array([r["angle_diff_deg"] for r in lines])
        print(f"\n--- Lines (n={len(lines)}) ---")
        print(f"  Length difference (mm):   mean={ld.mean():.2f}  "
              f"std={ld.std():.2f}  median={np.median(ld):.2f}  "
              f"min={ld.min():.2f}  max={ld.max():.2f}")
        print(f"  Angular difference (deg): mean={ad.mean():.2f}  "
              f"std={ad.std():.2f}  median={np.median(ad):.2f}  "
              f"min={ad.min():.2f}  max={ad.max():.2f}")

        if icc_lines_data:
            len1 = np.array([d[1] for d in icc_lines_data])
            len2 = np.array([d[2] for d in icc_lines_data])
            icc_len = calculate_icc_2rater(len1, len2)
            print(f"  ICC(3,1) length: {icc_len:.4f}" if icc_len is not None
                  else "  ICC(3,1) length: n/a")

        print("  Per annotator:")
        for ann in sorted(set(r["annotator"] for r in lines)):
            l_arr = np.array([r["length_diff_mm"]
                              for r in lines if r["annotator"] == ann])
            a_arr = np.array([r["angle_diff_deg"]
                              for r in lines if r["annotator"] == ann])
            print(f"    {ann:12s}  n={len(l_arr):3d}  "
                  f"length: mean={l_arr.mean():.2f} std={l_arr.std():.2f}  "
                  f"angle: mean={a_arr.mean():.2f} std={a_arr.std():.2f}")
    else:
        print("\n--- Lines: no matched pairs ---")

    # ── Finding-level agreement ──
    print(f"\n--- Finding-level agreement ---")
    session_diffs = []
    for (ann, p1, p2), (n1, n2) in finding_counts.items():
        session_diffs.append((ann, n1, n2, abs(n1 - n2)))

    total_sessions = len(session_diffs)
    matching = sum(1 for _, _, _, d in session_diffs if d == 0)
    abs_diffs = np.array([d for _, _, _, d in session_diffs])
    print(f"  Total session pairs: {total_sessions}")
    print(f"  Sessions with identical finding count: {matching} "
          f"({matching / total_sessions * 100:.1f}%)")
    print(f"  Mean abs difference in annotation count: "
          f"{abs_diffs.mean():.2f} +/- {abs_diffs.std():.2f}")

    print("  Per annotator:")
    ann_sessions = defaultdict(list)
    for ann, n1, n2, d in session_diffs:
        ann_sessions[ann].append(d)
    for ann in sorted(ann_sessions):
        diffs = np.array(ann_sessions[ann])
        m = sum(1 for d in diffs if d == 0)
        print(f"    {ann:12s}  sessions={len(diffs):3d}  "
              f"matching={m}/{len(diffs)}  "
              f"mean_diff={diffs.mean():.2f} +/- {diffs.std():.2f}")

    # ── Match summary ──
    n_unmatched = sum(1 for r in results if r["status"] == "unmatched")
    n_auto_unmatched = sum(1 for r in results if r["status"] == "auto_unmatched")
    n_mixed = sum(1 for r in results if r["status"] == "mixed_type")
    n_excluded = sum(1 for r in results if r["status"] == "excluded_zero_roi")
    n_auto_matched = sum(1 for r in results if r["match_method"] == "auto_matched")

    print(f"\n--- Match summary ---")
    print(f"  Direct 1:1 pairs compared: "
          f"{sum(1 for r in results if r['match_method'] == 'direct' and r['status'] == 'compared')}")
    print(f"  Auto-matched pairs (Hungarian): {n_auto_matched}")
    print(f"  Unmatched findings (one session only): {n_unmatched}")
    print(f"  Auto-unmatched (unequal list length): {n_auto_unmatched}")
    print(f"  Mixed-type pairs: {n_mixed}")
    print(f"  Excluded zero-size ROIs: {n_excluded}")

    if n_mixed > 0:
        print(f"\n  Mixed-type pair details:")
        for r in results:
            if r["status"] == "mixed_type":
                print(f"    {r['annotator']}: {r['annotation1']} ({r['type1']}) <> "
                      f"{r['annotation2']} ({r['type2']})")

    print("\n" + "=" * 70)


if __name__ == "__main__":
    main()