Upload 2 files

clinical_metrics_v2.py

@@ -0,0 +1,382 @@
+"""
+Phase 3 — Clinical Logic & Verification Engine
+===============================================
+Implements deterministic clinical metric extraction from segmentation masks.
+No ML here — pure geometry and ophthalmology math.
+
+Mask convention (from Phase 2 training):
+    0 = Background
+    1 = Optic Disc
+    2 = Optic Cup
+"""
+
+import numpy as np
+import cv2
+from dataclasses import dataclass, field
+from typing import Tuple, Dict, Optional
+from enum import Enum
+
+# ── ISNT tolerance margin ──────────────────────────────────────────────────
+# Exact I > S > N > T fails on tiny numerical noise.
+# A margin of 0.2 rim-pixels absorbs rounding without masking real violations.
+_ISNT_MARGIN = 0.2
+
+
+class RiskLevel(Enum):
+    HEALTHY = "Healthy"
+    SUSPECT = "Glaucoma Suspect"
+    HIGH    = "High Risk"
+
+
+class SanityError(Exception):
+    """Raised when mask geometry violates anatomical constraints."""
+    pass
+
+
+@dataclass
+class ISNTResult:
+    inferior:       float = 0.0
+    superior:       float = 0.0
+    nasal:          float = 0.0
+    temporal:       float = 0.0
+    rule_satisfied: bool  = False
+
+    def to_dict(self) -> Dict[str, float]:
+        return {
+            'inferior':       round(self.inferior,  4),
+            'superior':       round(self.superior,  4),
+            'nasal':          round(self.nasal,     4),
+            'temporal':       round(self.temporal,  4),
+            'rule_satisfied': self.rule_satisfied,
+        }
+
+
+@dataclass
+class ClinicalResult:
+    vcdr:             float      = 0.0
+    isnt:             ISNTResult = field(default_factory=ISNTResult)
+    disc_area_px:     int        = 0
+    cup_area_px:      int        = 0
+    disc_center:      Tuple[int, int] = (0, 0)
+    cup_center:       Tuple[int, int] = (0, 0)
+    uncertainty:      float      = 0.0
+    high_uncertainty: bool       = False
+    risk_level:       RiskLevel  = RiskLevel.HEALTHY
+    sanity_passed:    bool       = False
+    warnings:         list       = field(default_factory=list)
+
+    def to_dict(self) -> dict:
+        return {
+            'vcdr':             round(self.vcdr, 4),
+            'isnt':             self.isnt.to_dict(),
+            'disc_area_px':     self.disc_area_px,
+            'cup_area_px':      self.cup_area_px,
+            'disc_center':      self.disc_center,
+            'cup_center':       self.cup_center,
+            'uncertainty':      round(self.uncertainty, 6),
+            'high_uncertainty': self.high_uncertainty,
+            'risk_level':       self.risk_level.value,
+            'sanity_passed':    self.sanity_passed,
+            'warnings':         self.warnings,
+        }
+
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Step 1 — Sanity checks
+# ─────────────────────────────────────────────────────────────────────────────
+
+def run_sanity_checks(disc_mask: np.ndarray, cup_mask: np.ndarray) -> None:
+    """
+    Enforce anatomical constraints.  Raises SanityError on hard violations.
+
+    Checks:
+        1. Masks are binary (0/1 values only).
+        2. Disc region is non-empty.
+        3. Cup is 100 % contained inside the disc (hard anatomical law).
+        4. Disconnected regions — warn only; upstream _clean_* handles them.
+    """
+    # 1. Binary check
+    for name, mask in [('disc', disc_mask), ('cup', cup_mask)]:
+        unique_vals = np.unique(mask)
+        if not set(unique_vals).issubset({0, 1}):
+            raise SanityError(
+                f"{name} mask contains non-binary values: {unique_vals}"
+            )
+
+    # 2. Non-empty disc
+    if int(disc_mask.sum()) == 0:
+        raise SanityError("Optic disc mask is empty — segmentation failure.")
+
+    # 3. Cup ⊂ Disc
+    cup_outside_disc = np.logical_and(cup_mask == 1, disc_mask == 0)
+    if cup_outside_disc.any():
+        n_violation = int(cup_outside_disc.sum())
+        raise SanityError(
+            f"Cup extends outside disc boundary ({n_violation} pixels). "
+            "Anatomically impossible — reject segmentation."
+        )
+
+    # 4. Single connected component — warn only, do not reject.
+    # Small segmentation gaps from noisy boundaries are handled upstream
+    # by _clean_disc / _clean_cup.
+    for name, mask in [('disc', disc_mask), ('cup', cup_mask)]:
+        if mask.sum() == 0:
+            continue
+        n_labels, _ = cv2.connectedComponents(mask.astype(np.uint8))
+        if n_labels > 2:
+            pass  # upstream cleanup handles; avoid hard rejection here
+
+
+# ──────────────────────────��──────────────────────────────────────────────────
+# Step 2 — vCDR
+# ─────────────────────────────────────────────────────────────────────────────
+
+def calculate_vcdr(
+    disc_mask: np.ndarray,
+    cup_mask:  np.ndarray,
+) -> Tuple[float, dict]:
+    """
+    Calculate vertical Cup-to-Disc Ratio using vertical extrema.
+
+    Clinical basis:
+        Horizontal disc expansion is less indicative of early glaucoma.
+        vCDR is the primary screening metric.
+
+    Returns:
+        vcdr    (float): ratio in [0, 1]
+        details (dict):  raw pixel measurements for transparency
+    """
+    disc_rows = np.where(disc_mask.any(axis=1))[0]
+    cup_rows  = np.where(cup_mask.any(axis=1))[0]
+
+    if disc_rows.size == 0:
+        return 0.0, {}
+
+    disc_v_diam = int(disc_rows.max() - disc_rows.min() + 1)
+    cup_v_diam  = int(cup_rows.max()  - cup_rows.min()  + 1) if cup_rows.size > 0 else 0
+
+    vcdr = cup_v_diam / disc_v_diam if disc_v_diam > 0 else 0.0
+
+    details = {
+        'disc_top_px':    int(disc_rows.min()),
+        'disc_bottom_px': int(disc_rows.max()),
+        'disc_v_diam_px': disc_v_diam,
+        'cup_top_px':     int(cup_rows.min())  if cup_rows.size > 0 else None,
+        'cup_bottom_px':  int(cup_rows.max())  if cup_rows.size > 0 else None,
+        'cup_v_diam_px':  cup_v_diam,
+    }
+    return round(vcdr, 4), details
+
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Step 3 — ISNT Rule
+# ─────────────────────────────────────────────────────────────────────────────
+
+def _disc_centroid(disc_mask: np.ndarray) -> Tuple[int, int]:
+    M = cv2.moments(disc_mask.astype(np.uint8))
+    if M['m00'] == 0:
+        h, w = disc_mask.shape
+        return h // 2, w // 2
+    cy = int(M['m01'] / M['m00'])
+    cx = int(M['m10'] / M['m00'])
+    return cy, cx
+
+
+def _rim_thickness_in_quadrant(
+    quadrant_mask: np.ndarray,
+    disc_mask:     np.ndarray,
+    cup_mask:      np.ndarray,
+) -> float:
+    """
+    Mean Euclidean distance from cup boundary to disc boundary,
+    measured inside a specific quadrant.
+
+    Uses a distance transform on the disc interior so that each pixel
+    inside the disc gets its distance to the disc edge.
+    We then sample only rim pixels (disc=1, cup=0) in this quadrant.
+    """
+    rim         = np.logical_and(disc_mask == 1, cup_mask == 0).astype(np.uint8)
+    rim_in_quad = np.logical_and(rim, quadrant_mask).astype(np.uint8)
+
+    if rim_in_quad.sum() == 0:
+        return 0.0
+
+    dist_to_disc_edge = cv2.distanceTransform(
+        disc_mask.astype(np.uint8), cv2.DIST_L2, 5
+    )
+
+    thicknesses = dist_to_disc_edge[rim_in_quad == 1]
+    return float(np.mean(thicknesses))
+
+
+def calculate_isnt(
+    disc_mask: np.ndarray,
+    cup_mask:  np.ndarray,
+) -> ISNTResult:
+    """
+    Calculate neuro-retinal rim thickness in four ISNT quadrants.
+
+    Quadrant definition (image coordinates):
+        Superior  — top half    (rows < cy)
+        Inferior  — bottom half (rows >= cy)
+        Nasal     — right half  (cols >= cx)   [standard right eye convention]
+        Temporal  — left half   (cols < cx)
+
+    ISNT rule (with tolerance margin):
+        Inferior > Superior - margin
+        Superior > Nasal    - margin
+        Nasal    > Temporal - margin
+
+    Using a strict exact ordering falsely triggers violations when
+    quadrant thicknesses differ by sub-pixel amounts due to rounding.
+    A margin of _ISNT_MARGIN (0.2 px) absorbs noise without masking
+    real rim thinning (which presents as differences of several pixels).
+    """
+    h, w   = disc_mask.shape
+    cy, cx = _disc_centroid(disc_mask)
+
+    superior_q = np.zeros((h, w), dtype=bool)
+    inferior_q = np.zeros((h, w), dtype=bool)
+    nasal_q    = np.zeros((h, w), dtype=bool)
+    temporal_q = np.zeros((h, w), dtype=bool)
+
+    superior_q[:cy, :]  = True
+    inferior_q[cy:, :]  = True
+    nasal_q[:,   cx:]   = True
+    temporal_q[:, :cx]  = True
+
+    I = _rim_thickness_in_quadrant(inferior_q,  disc_mask, cup_mask)
+    S = _rim_thickness_in_quadrant(superior_q,  disc_mask, cup_mask)
+    N = _rim_thickness_in_quadrant(nasal_q,     disc_mask, cup_mask)
+    T = _rim_thickness_in_quadrant(temporal_q,  disc_mask, cup_mask)
+
+    # Tolerated ISNT check — absorbs sub-pixel numerical noise
+    rule_ok = (
+        I > S - _ISNT_MARGIN and
+        S > N - _ISNT_MARGIN and
+        N > T - _ISNT_MARGIN
+    )
+
+    return ISNTResult(
+        inferior=I, superior=S, nasal=N, temporal=T,
+        rule_satisfied=rule_ok,
+    )
+
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Step 4 — Risk classification
+# ─────────────────────────────────────────────────────────────────────────────
+
+def classify_risk(
+    vcdr:                 float,
+    isnt:                 ISNTResult,
+    uncertainty:          float,
+    uncertainty_threshold: float = 0.05,
+) -> Tuple[RiskLevel, list]:
+    """
+    Rule-based risk stratification.
+
+    Thresholds from clinical literature:
+        vCDR < 0.65              → Healthy
+        vCDR 0.65–0.80           → Suspect
+        vCDR > 0.80              → High Risk
+        ISNT violation (any risk) → escalate to at least Suspect
+        High uncertainty         → override to Suspect regardless of vCDR
+    """
+    warnings = []
+
+    if uncertainty > uncertainty_threshold:
+        warnings.append(
+            f"High model uncertainty ({uncertainty:.4f}) — result may be unreliable."
+        )
+        return RiskLevel.SUSPECT, warnings
+
+    if vcdr > 0.80:
+        risk = RiskLevel.HIGH
+        warnings.append(
+            f"vCDR {vcdr:.2f} exceeds 0.80 — urgent referral recommended."
+        )
+    elif vcdr > 0.65:
+        risk = RiskLevel.SUSPECT
+        warnings.append(f"vCDR {vcdr:.2f} in borderline range (0.65–0.80).")
+    else:
+        risk = RiskLevel.HEALTHY
+
+    if not isnt.rule_satisfied:
+        warnings.append("ISNT rule violated — neuro-retinal rim thinning detected.")
+        if risk == RiskLevel.HEALTHY:
+            risk = RiskLevel.SUSPECT
+
+    return risk, warnings
+
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Main pipeline entry point
+# ─────────────────────────────────────────────────────────────────────────────
+
+def run_clinical_pipeline(
+    disc_mask:            np.ndarray,
+    cup_mask:             np.ndarray,
+    uncertainty:          float = 0.0,
+    uncertainty_threshold: float = 0.05,
+) -> ClinicalResult:
+    """
+    Execute complete Phase 3 clinical pipeline on binary masks.
+
+    Args:
+        disc_mask:             uint8 binary array (1 = disc, 0 = background)
+        cup_mask:              uint8 binary array (1 = cup,  0 = background)
+        uncertainty:           scalar from Phase 2 MC-Dropout (ROI-restricted)
+        uncertainty_threshold: flag above this value as high uncertainty
+
+    Returns:
+        ClinicalResult dataclass
+    """
+    result = ClinicalResult()
+    result.uncertainty      = float(uncertainty)
+    result.high_uncertainty = uncertainty > uncertainty_threshold
+
+    disc_mask = (disc_mask > 0).astype(np.uint8)
+    cup_mask  = (cup_mask  > 0).astype(np.uint8)
+
+    # ── Sanity checks ──────────────────────────────────────────────────
+    try:
+        run_sanity_checks(disc_mask, cup_mask)
+        result.sanity_passed = True
+    except SanityError as e:
+        result.warnings.append(f"SANITY FAIL: {e}")
+        result.sanity_passed = False
+        # Only abort if there is no disc at all — otherwise continue
+        # computing metrics on whatever geometry we have.
+        if disc_mask.sum() == 0:
+            result.risk_level = RiskLevel.SUSPECT
+            return result
+
+    # ── Structural measurements ────────────────────────────────────────
+    result.disc_area_px = int(disc_mask.sum())
+    result.cup_area_px  = int(cup_mask.sum())
+
+    dy, dx = _disc_centroid(disc_mask)
+    result.disc_center = (int(dx), int(dy))
+
+    if cup_mask.sum() > 0:
+        M = cv2.moments(cup_mask)
+        if M['m00'] > 0:
+            result.cup_center = (
+                int(M['m10'] / M['m00']),
+                int(M['m01'] / M['m00']),
+            )
+
+    # ── vCDR ──────────────────────────────────────────────────────────
+    result.vcdr, _ = calculate_vcdr(disc_mask, cup_mask)
+
+    # ── ISNT (with tolerance margin) ──────────────────────────────────
+    result.isnt = calculate_isnt(disc_mask, cup_mask)
+
+    # ── Risk classification ─────���─────────────────────────────────────
+    result.risk_level, warnings = classify_risk(
+        result.vcdr, result.isnt, uncertainty, uncertainty_threshold
+    )
+    result.warnings.extend(warnings)
+
+    return result

phase3pipeline_v2.py

@@ -0,0 +1,351 @@
+"""
+Phase 3 — Inference Pipeline
+"""
+
+import os
+import numpy as np
+import cv2
+import torch
+import torch.nn.functional as F
+from typing import Optional, Dict, Any, Tuple
+
+from model import UNet
+from checkpoint_loader import load_model_for_inference
+from clinical_metrics import run_clinical_pipeline, ClinicalResult
+
+# ── Minimum plausible disc area in pixels (512×512 image).
+# Anything smaller is almost certainly not a real fundus image.
+# A disc typically covers ~3–5 % of the 512×512 canvas ≈ 7,000–13,000 px.
+# We gate at 1 % (≈ 2,600 px) to be conservative.
+_MIN_DISC_AREA_PX  = 2_600
+_MIN_CUP_AREA_PX   = 100       # below this the cup reading is meaningless
+_MIN_CUP_DISC_RATIO = 0.01     # cup/disc area fraction — below = unreliable
+
+
+class Phase3Pipeline:
+
+    def __init__(
+        self,
+        repo_id: str = "Nj-1111/EyeeSEE",
+        epoch: Optional[int] = None,
+        mc_passes: int = 20,
+        uncertainty_threshold: float = 0.05,
+        device: Optional[torch.device] = None,
+        token: Optional[str] = None,
+        debug: bool = False,          # gate all debug I/O behind this flag
+    ):
+        self.repo_id               = repo_id
+        self.mc_passes             = mc_passes
+        self.uncertainty_threshold = uncertainty_threshold
+        self.debug                 = debug
+
+        self.device = device or torch.device(
+            'cuda' if torch.cuda.is_available() else 'cpu'
+        )
+
+        self.token = (
+            token or
+            os.getenv('HF_TOKEN_2') or
+            os.getenv('HF_TOKEN')
+        )
+
+        # IMP: lower dropout improves stability
+        self.model = UNet(
+            in_channels=1,
+            n_classes=3,
+            base_filters=64,
+            dropout=0.1
+        )
+
+        load_model_for_inference(
+            model=self.model,
+            repo_id=repo_id,
+            epoch=epoch,
+            device=self.device,
+            token=self.token
+        )
+
+    # ──────────────────────────────────────────────────────────────────────
+    # preprocessing
+    # ──────────────────────────────────────────────────────────────────────
+
+    def _preprocess(self, image: np.ndarray) -> torch.Tensor:
+
+        if image.ndim == 3:
+            image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
+
+        image = cv2.resize(
+            image,
+            (512, 512),
+            interpolation=cv2.INTER_AREA
+        )
+
+        image = image.astype(np.float32) / 255.0
+
+        return (
+            torch.from_numpy(image)
+            .unsqueeze(0)
+            .unsqueeze(0)
+            .to(self.device)
+        )
+
+    # ──────────────────────────────────────────────────────────────────────
+    # mask cleanup  (separate functions — disc and cup have different scales)
+    # ──────────────────────────────────────────────────────────────────────
+
+    def _clean_disc(self, binary_mask: np.ndarray) -> np.ndarray:
+        """
+        Morphological cleanup for the disc mask.
+        Disc is large enough that a 5×5 kernel is safe.
+        """
+        binary_mask = binary_mask.astype(np.uint8)
+        if binary_mask.sum() == 0:
+            return binary_mask
+
+        kernel = np.ones((5, 5), np.uint8)
+        binary_mask = cv2.morphologyEx(binary_mask, cv2.MORPH_OPEN,  kernel)
+        binary_mask = cv2.morphologyEx(binary_mask, cv2.MORPH_CLOSE, kernel)
+
+        n, labels, stats, _ = cv2.connectedComponentsWithStats(binary_mask)
+        if n <= 1:
+            return binary_mask
+        largest = 1 + np.argmax(stats[1:, cv2.CC_STAT_AREA])
+        cleaned = np.zeros_like(binary_mask)
+        cleaned[labels == largest] = 1
+        return cleaned
+
+    def _clean_cup(self, binary_mask: np.ndarray) -> np.ndarray:
+        """
+        Morphological cleanup for the cup mask.
+        Uses an adaptive kernel: after anatomical clipping the remaining cup
+        may be small — a 5×5 open would erase it.  We drop to 3×3 for
+        small remnants.
+        """
+        binary_mask = binary_mask.astype(np.uint8)
+        if binary_mask.sum() == 0:
+            return binary_mask
+
+        k = 3 if binary_mask.sum() < 3_000 else 5
+        kernel = np.ones((k, k), np.uint8)
+        binary_mask = cv2.morphologyEx(binary_mask, cv2.MORPH_OPEN,  kernel)
+        binary_mask = cv2.morphologyEx(binary_mask, cv2.MORPH_CLOSE, kernel)
+
+        n, labels, stats, _ = cv2.connectedComponentsWithStats(binary_mask)
+        if n <= 1:
+            return binary_mask
+        largest = 1 + np.argmax(stats[1:, cv2.CC_STAT_AREA])
+        cleaned = np.zeros_like(binary_mask)
+        cleaned[labels == largest] = 1
+        return cleaned
+
+    # ──────────────────────────────────────────────────────────────────────
+    # anatomical enforcement
+    # ──────────────────────────────────────────────────────────────────────
+
+    def _enforce_cup_in_disc(
+        self,
+        disc_mask: np.ndarray,
+        cup_mask: np.ndarray,
+    ) -> Tuple[np.ndarray, np.ndarray, int]:
+        """
+        Hard-clip cup to disc boundary.
+        Returns (disc_mask, corrected_cup, n_pixels_removed).
+        Pure logical AND — no morphology here.
+        """
+        corrected_cup = np.logical_and(cup_mask == 1, disc_mask == 1).astype(np.uint8)
+        violations    = int(cup_mask.sum()) - int(corrected_cup.sum())
+        return disc_mask, corrected_cup, violations
+
+    # ──────────────────────────────────────────────────────────────────────
+    # MC-Dropout segmentation
+    # ──────────────────────────────────────────────────────────────────────
+
+    def _mc_segment(
+        self,
+        tensor: torch.Tensor,
+    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """
+        Run MC-Dropout inference and return:
+            disc_mask (raw, thresholded, NOT yet cleaned)
+            cup_mask  (raw, thresholded, NOT yet cleaned)
+            var_probs (3, H, W)  — full variance map, used for ROI uncertainty later
+
+        Masks are intentionally returned uncleaned so that the caller can
+        enforce topology FIRST and clean AFTER (see run()).
+        """
+        # BatchNorm stays frozen; only dropout layers go to train mode
+        self.model.eval()
+        for m in self.model.modules():
+            if isinstance(m, (torch.nn.Dropout, torch.nn.Dropout2d, torch.nn.Dropout3d)):
+                m.train()
+
+        all_probs = []
+        with torch.no_grad():
+            for _ in range(self.mc_passes):
+                logits = self.model(tensor)
+                probs  = F.softmax(logits, dim=1)
+                all_probs.append(probs.cpu().numpy())
+
+        self.model.eval()
+
+        all_probs  = np.stack(all_probs, axis=0)           # (T, 1, 3, H, W)
+        mean_probs = all_probs.mean(axis=0)[0]             # (3, H, W)
+        var_probs  = all_probs.var(axis=0)[0]              # (3, H, W)
+
+        # ── Separate thresholds — critical fix ─────────────────────────
+        # Cup probability maps are more diffuse than disc maps.
+        # Using 0.25 for both causes the cup to spill far outside the disc;
+        # the anatomical clipping then wipes it out entirely.
+        # A higher cup threshold trims diffuse edges back inside the disc.
+        disc_mask = (mean_probs[1] > 0.35).astype(np.uint8)
+        cup_mask  = (mean_probs[2] > 0.55).astype(np.uint8)
+
+        if self.debug:
+            cv2.imwrite("disc_prob_debug.png", (mean_probs[1] * 255).astype(np.uint8))
+            cv2.imwrite("cup_prob_debug.png",  (mean_probs[2] * 255).astype(np.uint8))
+            print(f"[debug] raw disc px={disc_mask.sum()}  raw cup px={cup_mask.sum()}")
+
+        return disc_mask, cup_mask, var_probs
+
+    # ──────────────────────────────────────────────────────────────────────
+    # uncertainty over anatomical ROI
+    # ──────────────────────────────────────────────────────────────────────
+
+    @staticmethod
+    def _roi_uncertainty(
+        var_probs: np.ndarray,      # (3, H, W)
+        disc_mask: np.ndarray,      # (H, W)
+        cup_mask:  np.ndarray,      # (H, W)
+    ) -> float:
+        """
+        Mean MC-Dropout variance restricted to the disc+cup region.
+
+        Averaging variance over the WHOLE image (including background)
+        suppresses clinically important local uncertainty because the model
+        is very confident about the large background class.
+        Restricting to the anatomical ROI makes the score meaningful.
+        """
+        roi = (disc_mask == 1) | (cup_mask == 1)
+        if not roi.any():
+            # No anatomical region found — fall back to global (will be high)
+            return float(var_probs[1:].mean())
+        # Channels 1 (disc) and 2 (cup) variance, sampled at ROI pixels
+        return float(var_probs[1:, roi].mean())
+
+    # ──────────────────────────────────────────────────────────────────────
+    # public API
+    # ──────────────────────────────────────────────────────────────────────
+
+    def run(self, image: np.ndarray) -> Dict[str, Any]:
+        """
+        Full pipeline.  Processing order (corrected):
+
+            1. threshold (separate disc / cup thresholds)
+            2. enforce topology  ← BEFORE any morphology
+            3. clean disc + clean cup  (separate kernels)
+            4. enforce topology again  ← morphology can reintroduce violations
+            5. compute ROI uncertainty
+            6. minimum-cup sanity gate
+            7. clinical metrics
+        """
+        tensor = self._preprocess(image)
+
+        # Step 1+2: raw threshold → first topology enforcement
+        # Cup is intentionally NOT cleaned yet so morphology doesn't expand
+        # it past the disc before the AND clip.
+        disc_mask, cup_mask, var_probs = self._mc_segment(tensor)
+
+        disc_mask, cup_mask, violations = self._enforce_cup_in_disc(disc_mask, cup_mask)
+
+        # Step 3: clean both masks (cup with adaptive small kernel)
+        disc_mask = self._clean_disc(disc_mask)
+        cup_mask  = self._clean_cup(cup_mask)
+
+        # Step 4: second enforcement — morphology can re-expand cup slightly
+        # outside a cleaned (potentially slightly shrunk) disc
+        disc_mask, cup_mask, extra_violations = self._enforce_cup_in_disc(disc_mask, cup_mask)
+        violations += extra_violations
+
+        # Step 5: uncertainty over anatomical ROI (not the whole image)
+        uncertainty = self._roi_uncertainty(var_probs, disc_mask, cup_mask)
+
+        if self.debug:
+            print(
+                f"[debug] post-clean disc px={disc_mask.sum()}  "
+                f"cup px={cup_mask.sum()}  "
+                f"violations={violations}  "
+                f"roi_uncertainty={uncertainty:.6f}"
+            )
+
+        # ── Step 6: minimum-cup sanity gate ───────────────────────────
+        # A tiny cup remnant (e.g. a few dozen pixels surviving morphology)
+        # produces a meaningless vCDR.  Zero it out and warn instead.
+        disc_area = int(disc_mask.sum())
+        cup_area  = int(cup_mask.sum())
+        extra_warnings: list[str] = []
+
+        if disc_area < _MIN_DISC_AREA_PX:
+            extra_warnings.append(
+                f"Disc too small ({disc_area} px) — image may not be a fundus photo. "
+                "Segmentation result is unreliable."
+            )
+
+        if cup_area > 0 and cup_area < _MIN_CUP_AREA_PX:
+            extra_warnings.append(
+                f"Cup remnant too small ({cup_area} px) — suppressed. "
+                "Cup segmentation is unreliable for this image."
+            )
+            cup_mask = np.zeros_like(cup_mask)
+            cup_area = 0
+
+        if cup_area > 0 and disc_area > 0:
+            if (cup_area / disc_area) < _MIN_CUP_DISC_RATIO:
+                extra_warnings.append(
+                    f"Cup/disc area ratio ({cup_area}/{disc_area} = "
+                    f"{cup_area/disc_area:.3f}) is below minimum — "
+                    "cup reading may be unreliable."
+                )
+
+        # Step 7: clinical pipeline
+        clinical = run_clinical_pipeline(
+            disc_mask=disc_mask,
+            cup_mask=cup_mask,
+            uncertainty=uncertainty,
+            uncertainty_threshold=self.uncertainty_threshold
+        )
+
+        # Inject all accumulated warnings
+        if violations > 0:
+            clinical.warnings.append(
+                f"Mask corrected: {violations} cup pixels clipped to disc boundary."
+            )
+        clinical.warnings.extend(extra_warnings)
+
+        report = {
+            'vcdr':             clinical.vcdr,
+            'isnt':             clinical.isnt.to_dict(),
+            'risk_level':       clinical.risk_level.value,
+            'uncertainty':      round(uncertainty, 6),
+            'high_uncertainty': clinical.high_uncertainty,
+            'disc_area_px':     clinical.disc_area_px,
+            'cup_area_px':      clinical.cup_area_px,
+            'disc_center':      clinical.disc_center,
+            'cup_center':       clinical.cup_center,
+            'sanity_passed':    clinical.sanity_passed,
+            'warnings':         clinical.warnings,
+        }
+
+        return {
+            'disc_mask':   disc_mask,
+            'cup_mask':    cup_mask,
+            'uncertainty': uncertainty,
+            'clinical':    clinical,
+            'report':      report,
+        }
+
+    def run_from_path(self, image_path: str) -> Dict[str, Any]:
+
+        image = cv2.imread(image_path)
+        if image is None:
+            raise FileNotFoundError(f"Cannot load image: {image_path}")
+        return self.run(image)