{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd, cv2, numpy as np, matplotlib.pyplot as plt\nROOT = \"/kaggle/input/competitions/aptos2019-blindness-detection\"\ndf = pd.read_csv(f\"{ROOT}/train.csv\")\ndf[\"path\"] = ROOT + \"/train_images/\" + df.id_code + \".png\"\nprint(df.diagnosis.value_counts().sort_index())","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-08-30T08:54:35.039407Z","iopub.execute_input":"2026-08-30T08:54:35.039741Z","iopub.status.idle":"2026-08-30T08:54:36.880399Z","shell.execute_reply.started":"2026-08-30T08:54:35.039708Z","shell.execute_reply":"2026-08-30T08:54:36.879626Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import glob, os\nprint(os.listdir(ROOT))                      # what CSVs/dirs actually exist\nimgs = glob.glob(f\"{ROOT}/**/*.png\", recursive=True)[:20]\nprint(len(glob.glob(f\"{ROOT}/**/*.png\", recursive=True)), \"images total\")\nfor p in imgs[:8]:\n    print(cv2.imread(p).shape, os.path.basename(p))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T08:54:36.881899Z","iopub.execute_input":"2026-08-30T08:54:36.882235Z","iopub.status.idle":"2026-08-30T08:54:45.324699Z","shell.execute_reply.started":"2026-08-30T08:54:36.882204Z","shell.execute_reply":"2026-08-30T08:54:45.323969Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile fundus_quality.py\n\"\"\"\nfundus_quality.py\nImage Quality Assessment for fundus photographs.\n\nCovers requirement 1 of the problem statement: decide whether an image is\ngradeable, and if not, say *why* so the operator can recapture.\n\nDesign notes:\n  - All metrics are computed on a size-normalised copy (short side -> WORK_SIZE)\n    so that scores are comparable across APTOS's wildly varying resolutions.\n  - All metrics are computed *inside the FOV mask only*. The black surround\n    would otherwise dominate every statistic.\n  - Thresholds in QualityThresholds are PLACEHOLDERS. Calibrate them on real\n    data (see calibrate_thresholds) before quoting any numbers to a judge.\n\nNo deep learning. Runs in ~20ms per image on CPU.\n\"\"\"\n\nfrom dataclasses import dataclass, asdict, field\nimport cv2\nimport numpy as np\n\nWORK_SIZE = 512  # short side, in px, that metrics are computed at\n\n\n# --------------------------------------------------------------------------\n# Field of view\n# --------------------------------------------------------------------------\n\ndef fov_mask(bgr, thresh=12):\n    \"\"\"Binary mask of the illuminated retinal disc.\n\n    Uses the max across channels rather than a grey conversion: in very dark\n    or very red images the green channel alone can fall below threshold inside\n    genuine retina.\n    \"\"\"\n    m = bgr.max(axis=2)\n    m = cv2.medianBlur(m, 5)\n    _, mask = cv2.threshold(m, thresh, 255, cv2.THRESH_BINARY)\n\n    k = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (9, 9))\n    mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, k)\n    mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, k)\n\n    # keep only the largest blob (drops timestamp text, lens flare specks)\n    n, labels, stats, _ = cv2.connectedComponentsWithStats(mask, 8)\n    if n <= 1:\n        return np.zeros(mask.shape, np.uint8)\n    largest = 1 + np.argmax(stats[1:, cv2.CC_STAT_AREA])\n    return ((labels == largest) * 255).astype(np.uint8)\n\n\ndef crop_to_fov(bgr, mask=None, pad=0):\n    \"\"\"Tight-crop to the retinal disc. Returns (cropped_bgr, cropped_mask).\"\"\"\n    if mask is None:\n        mask = fov_mask(bgr)\n    ys, xs = np.where(mask > 0)\n    if len(xs) == 0:\n        return bgr, mask\n    y0, y1 = max(ys.min() - pad, 0), min(ys.max() + pad + 1, bgr.shape[0])\n    x0, x1 = max(xs.min() - pad, 0), min(xs.max() + pad + 1, bgr.shape[1])\n    return bgr[y0:y1, x0:x1], mask[y0:y1, x0:x1]\n\n\ndef _resize_short_side(img, target=WORK_SIZE, interp=cv2.INTER_AREA):\n    h, w = img.shape[:2]\n    s = target / min(h, w)\n    if s >= 1.0:\n        return img\n    return cv2.resize(img, (int(round(w * s)), int(round(h * s))), interpolation=interp)\n\n\n# --------------------------------------------------------------------------\n# Individual metrics\n# --------------------------------------------------------------------------\n\ndef focus_score(green, mask):\n    \"\"\"Variance of the Laplacian inside the FOV. Higher = sharper.\n\n    Green channel because vessels and haemorrhages have the most contrast there.\n    \"\"\"\n    lap = cv2.Laplacian(green, cv2.CV_64F, ksize=3)\n    vals = lap[mask > 0]\n    return float(vals.var()) if vals.size else 0.0\n\n\ndef illumination_scores(green, mask, grid=3):\n    \"\"\"Uniformity of illumination across a grid over the FOV.\n\n    Returns (cv_of_cell_means, frac_saturated, frac_crushed).\n    A high CV means one side of the image is much brighter than the other,\n    which is the classic signature of a misaligned handheld capture.\n    \"\"\"\n    h, w = green.shape\n    means = []\n    for i in range(grid):\n        for j in range(grid):\n            ys, ye = i * h // grid, (i + 1) * h // grid\n            xs, xe = j * w // grid, (j + 1) * w // grid\n            cell = green[ys:ye, xs:xe]\n            cm = mask[ys:ye, xs:xe]\n            if (cm > 0).sum() > 0.2 * cm.size:  # ignore mostly-black corners\n                means.append(cell[cm > 0].mean())\n\n    means = np.array(means, dtype=np.float64)\n    cv_ = float(means.std() / (means.mean() + 1e-6)) if means.size >= 2 else 0.0\n\n    inside = green[mask > 0]\n    if inside.size == 0:\n        return cv_, 0.0, 0.0\n    return cv_, float((inside >= 250).mean()), float((inside <= 10).mean())\n\n\ndef contrast_score(green, mask):\n    \"\"\"Robust contrast: 5th-to-95th percentile spread inside the FOV.\"\"\"\n    vals = green[mask > 0]\n    if vals.size == 0:\n        return 0.0\n    p5, p95 = np.percentile(vals, [5, 95])\n    return float(p95 - p5)\n\n\ndef framing_scores(mask):\n    \"\"\"Frame fill and edge clipping, measured separately per axis.\n\n    Vertical and horizontal clipping mean completely different things.\n    Standard fundus cameras output a wide frame (aspect ~1.4) in which the\n    retinal disc is cropped flat at top and bottom with black bars left and\n    right. That is the normal format, not a defect, so clip_v is expected to\n    be high on most clinical images.\n\n    Horizontal clipping is the real signal: it means the disc ran off the\n    left or right edge, which happens when the operator was off-centre.\n\n    Returns (coverage, clip_v, clip_h).\n    \"\"\"\n    h, w = mask.shape\n    coverage = float((mask > 0).mean())\n\n    top = (mask[0, :] > 0).sum()\n    bottom = (mask[-1, :] > 0).sum()\n    left = (mask[:, 0] > 0).sum()\n    right = (mask[:, -1] > 0).sum()\n\n    clip_v = float((top + bottom) / (2.0 * w))\n    clip_h = float((left + right) / (2.0 * h))\n    return coverage, clip_v, clip_h\n\n\n# --------------------------------------------------------------------------\n# Verdict\n# --------------------------------------------------------------------------\n\n@dataclass\nclass QualityThresholds:\n    \"\"\"PLACEHOLDER values. Run calibrate_thresholds() on APTOS before trusting.\"\"\"\n    focus_reject: float = 40.0\n    focus_usable: float = 120.0\n    illum_cv_reject: float = 0.45\n    illum_cv_usable: float = 0.28\n    saturated_reject: float = 0.06\n    crushed_reject: float = 0.55\n    contrast_reject: float = 26.0\n    contrast_usable: float = 36.0\n    coverage_reject: float = 0.18\n    # Physical judgements, deliberately NOT percentile-calibrated: a retina\n    # running off the side of the frame is bad regardless of how common it\n    # happens to be in whatever dataset you calibrated on.\n    clip_h_reject: float = 0.40     # disc runs off left/right edge\n    clip_h_usable: float = 0.15\n    clip_v_usable: float = 0.90     # near-total vertical fill; peripheral loss\n\n\n@dataclass\nclass QualityReport:\n    verdict: str                     # \"good\" | \"usable\" | \"reject\"\n    reasons: list = field(default_factory=list)   # operator-facing messages\n    metrics: dict = field(default_factory=dict)   # raw numbers, for the report/audit\n\n    def as_dict(self):\n        return asdict(self)\n\n\ndef assess_quality(bgr, thresholds=None):\n    \"\"\"Full quality assessment. Input: BGR uint8 image as read by cv2.imread.\"\"\"\n    t = thresholds or QualityThresholds()\n\n    small = _resize_short_side(bgr)\n    mask = fov_mask(small)\n\n    if (mask > 0).sum() < 0.02 * mask.size:\n        return QualityReport(\n            \"reject\",\n            [\"No retina detected. Check that the lens is aligned with the pupil.\"],\n            {},\n        )\n\n    green = small[:, :, 1]\n    focus = focus_score(green, mask)\n    illum_cv, sat, crushed = illumination_scores(green, mask)\n    contrast = contrast_score(green, mask)\n    coverage, clip_v, clip_h = framing_scores(mask)\n\n    metrics = dict(\n        focus=round(focus, 2),\n        illum_cv=round(illum_cv, 4),\n        saturated_frac=round(sat, 4),\n        crushed_frac=round(crushed, 4),\n        contrast=round(contrast, 2),\n        fov_coverage=round(coverage, 4),\n        clip_v=round(clip_v, 4),\n        clip_h=round(clip_h, 4),\n    )\n\n    reject, warn = [], []\n\n    if focus < t.focus_reject:\n        reject.append(\"Image is out of focus. Hold steady and refocus before recapturing.\")\n    elif focus < t.focus_usable:\n        warn.append(\"Slightly soft focus.\")\n\n    if sat > t.saturated_reject:\n        reject.append(\"Overexposed with glare. Reduce illumination or change the angle.\")\n    if crushed > t.crushed_reject:\n        reject.append(\"Image is too dark. Increase illumination or dilate the pupil.\")\n\n    if illum_cv > t.illum_cv_reject:\n        reject.append(\"Uneven illumination across the retina. Re-centre the lens on the pupil.\")\n    elif illum_cv > t.illum_cv_usable:\n        warn.append(\"Uneven illumination.\")\n\n    if contrast < t.contrast_reject:\n        reject.append(\"Very low contrast; retinal detail is not visible.\")\n    elif contrast < t.contrast_usable:\n        warn.append(\"Low contrast.\")\n\n    if coverage < t.coverage_reject:\n        reject.append(\"Retina fills too little of the frame. Move closer to the eye.\")\n    if clip_h > t.clip_h_reject:\n        reject.append(\"Retina runs off the side of the frame. Re-centre the camera horizontally.\")\n    elif clip_h > t.clip_h_usable:\n        warn.append(\"Retina slightly off-centre horizontally.\")\n\n    if clip_v > t.clip_v_usable:\n        warn.append(\"Frame is cropped tight; peripheral retina may be missing.\")\n\n    if reject:\n        return QualityReport(\"reject\", reject, metrics)\n    if warn:\n        return QualityReport(\"usable\", warn, metrics)\n    return QualityReport(\"good\", [\"Gradeable.\"], metrics)\n\n\n# --------------------------------------------------------------------------\n# Threshold calibration\n# --------------------------------------------------------------------------\n\ndef load_thresholds(path):\n    \"\"\"Read calibrate_thresholds() output back into a QualityThresholds.\n\n    Only the distribution-relative fields are overridden. The exposure and\n    framing limits (saturated, crushed, coverage, clipping) stay at their\n    defaults on purpose: those are physical facts about an image, not\n    properties of a particular dataset, so calibrating them to APTOS would\n    make them wrong for your handheld rig.\n    \"\"\"\n    import json\n    with open(path) as f:\n        cal = json.load(f)\n    t = QualityThresholds()\n    for k, v in cal.items():\n        if hasattr(t, k):\n            setattr(t, k, float(v))\n    return t\n\n\ndef calibrate_thresholds(metric_rows, reject_pct=5, usable_pct=25):\n    \"\"\"Derive thresholds from the observed distribution instead of guessing.\n\n    metric_rows: list of the `metrics` dicts returned by assess_quality.\n\n    Rationale to give a judge: 'we set the reject threshold at the 5th\n    percentile of the APTOS distribution, then hand-reviewed the flagged\n    images to confirm the cut point.' That is defensible. Magic numbers\n    are not.\n    \"\"\"\n    import pandas as pd\n    df = pd.DataFrame(metric_rows)\n    return {\n        \"focus_reject\":    float(df.focus.quantile(reject_pct / 100)),\n        \"focus_usable\":    float(df.focus.quantile(usable_pct / 100)),\n        \"contrast_reject\": float(df.contrast.quantile(reject_pct / 100)),\n        \"contrast_usable\": float(df.contrast.quantile(usable_pct / 100)),\n        \"illum_cv_reject\": float(df.illum_cv.quantile(1 - reject_pct / 100)),\n        \"illum_cv_usable\": float(df.illum_cv.quantile(1 - usable_pct / 100)),\n    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T08:54:45.325685Z","iopub.execute_input":"2026-08-30T08:54:45.325913Z","iopub.status.idle":"2026-08-30T08:54:45.337575Z","shell.execute_reply.started":"2026-08-30T08:54:45.325890Z","shell.execute_reply":"2026-08-30T08:54:45.336853Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile fundus_dedup.py\n\"\"\"\nfundus_dedup.py\nNear-duplicate detection for fundus images.\n\nWhy plain dHash fails here\n--------------------------\nA 9x8 dHash of a fundus photo encodes the circular FOV boundary and nothing\nelse -- every retina is a bright disc on black, so every hash is nearly\nidentical and union-find chains the whole dataset into one group.\n\nWhat this does instead\n----------------------\n1. Crop away the black surround, so the frame stops dominating.\n2. Take the green channel and divide out a heavy blur. This removes the\n   global orange cast and per-camera illumination, leaving vessel structure,\n   which is what actually distinguishes two retinas.\n3. Hash at 33x32 -> 1024 bits instead of 64. Two orders of magnitude more\n   discriminative.\n4. VERIFY every candidate pair by direct correlation before merging.\n   This is the step that stops chaining: a false candidate is rejected\n   before union-find ever sees it.\n\"\"\"\n\nfrom concurrent.futures import ProcessPoolExecutor\n\nimport cv2\nimport numpy as np\nimport pandas as pd\n\nfrom fundus_quality import fov_mask, crop_to_fov\n\nHASH_W, HASH_H = 33, 32      # -> 32*32 = 1024 bits\nTHUMB = 64                   # verification thumbnail size\nHAMMING_FRAC = 0.10          # candidate if < 10% of bits differ\nNCC_MIN = 0.90               # verified duplicate if correlation >= this\n\n\ndef _normalised(path):\n    \"\"\"FOV-cropped, illumination-flattened green channel.\"\"\"\n    img = cv2.imread(path, cv2.IMREAD_REDUCED_COLOR_4)\n    if img is None:\n        img = cv2.imread(path)\n    if img is None:\n        return None\n    img, _ = crop_to_fov(img, fov_mask(img))\n    if min(img.shape[:2]) < 16:\n        return None\n\n    g = img[:, :, 1].astype(np.float32)\n    sigma = max(g.shape) / 20.0\n    bg = cv2.GaussianBlur(g, (0, 0), sigmaX=sigma)\n    flat = g / (bg + 1e-3)                      # retinex-style flattening\n    return np.clip(flat, 0.5, 1.5)\n\n\ndef fingerprint(path):\n    \"\"\"Returns (1024-bit hash as uint8 array, 64x64 z-scored thumbnail).\"\"\"\n    flat = _normalised(path)\n    if flat is None:\n        return np.zeros(HASH_H * (HASH_W - 1), np.uint8), np.zeros((THUMB, THUMB), np.float32)\n\n    small = cv2.resize(flat, (HASH_W, HASH_H), interpolation=cv2.INTER_AREA)\n    bits = (small[:, 1:] > small[:, :-1]).flatten().astype(np.uint8)\n\n    thumb = cv2.resize(flat, (THUMB, THUMB), interpolation=cv2.INTER_AREA)\n    thumb = (thumb - thumb.mean()) / (thumb.std() + 1e-6)\n    return bits, thumb.astype(np.float32)\n\n\ndef find_duplicate_groups(df, hamming_frac=HAMMING_FRAC, ncc_min=NCC_MIN,\n                          workers=4, verbose=True):\n    with ProcessPoolExecutor(workers) as ex:\n        fps = list(ex.map(fingerprint, df.path.tolist(), chunksize=16))\n\n    H = np.array([f[0] for f in fps], np.float32)          # (N, 1024)\n    T = np.array([f[1] for f in fps], np.float32).reshape(len(fps), -1)  # (N, 4096)\n    n_bits = H.shape[1]\n\n    # Hamming distance without an (N, N, bits) tensor:\n    #   d(a,b) = |a| + |b| - 2 a.b   for binary vectors\n    s = H.sum(1)\n    d = s[:, None] + s[None, :] - 2.0 * (H @ H.T)\n    np.fill_diagonal(d, n_bits)\n\n    cand_i, cand_j = np.where(d <= hamming_frac * n_bits)\n    keep = cand_i < cand_j\n    cand_i, cand_j = cand_i[keep], cand_j[keep]\n    if verbose:\n        print(f\"  {len(cand_i)} candidate pairs from hashing\")\n\n    # verification: normalised cross-correlation of the z-scored thumbnails\n    if len(cand_i):\n        ncc = (T[cand_i] * T[cand_j]).sum(1) / T.shape[1]\n        ok = ncc >= ncc_min\n        if verbose:\n            print(f\"  {int(ok.sum())} survive correlation check \"\n                  f\"(rejected {int((~ok).sum())} as false matches)\")\n        cand_i, cand_j = cand_i[ok], cand_j[ok]\n\n    parent = list(range(len(df)))\n\n    def find(x):\n        while parent[x] != x:\n            parent[x] = parent[parent[x]]\n            x = parent[x]\n        return x\n\n    for i, j in zip(cand_i, cand_j):\n        ri, rj = find(int(i)), find(int(j))\n        if ri != rj:\n            parent[ri] = rj\n\n    out = df[[\"id_code\", \"diagnosis\"]].copy()\n    out[\"group_id\"] = [find(i) for i in range(len(df))]\n\n    sizes = out.group_id.value_counts()\n    dup = sizes[sizes > 1]\n    if verbose:\n        print(f\"  {len(dup)} duplicate groups, {int(dup.sum())} images, \"\n              f\"largest group {int(sizes.max())}\")\n        conflict = int((out.groupby('group_id').diagnosis.nunique() > 1).sum())\n        if conflict:\n            print(f\"  {conflict} groups have conflicting labels\")\n        if sizes.max() > 10:\n            print(\"  !! still chaining -- raise ncc_min toward 0.95\")\n    return out\n\n\ndef inspect_group(df, dupes, group_id, out_png=\"dup_group.png\"):\n    \"\"\"Save a strip of every image in a group so you can eyeball the call.\"\"\"\n    ids = dupes.loc[dupes.group_id == group_id, \"id_code\"].tolist()\n    paths = df.set_index(\"id_code\").loc[ids, \"path\"].tolist()\n    tiles = []\n    for p in paths[:8]:\n        im = cv2.imread(p, cv2.IMREAD_REDUCED_COLOR_8)\n        im, _ = crop_to_fov(im, fov_mask(im))\n        tiles.append(cv2.resize(im, (200, 200)))\n    cv2.imwrite(out_png, np.hstack(tiles))\n    return ids","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T08:54:45.339225Z","iopub.execute_input":"2026-08-30T08:54:45.339551Z","iopub.status.idle":"2026-08-30T08:54:45.359004Z","shell.execute_reply.started":"2026-08-30T08:54:45.339516Z","shell.execute_reply":"2026-08-30T08:54:45.358152Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nstep2_prepare.py\nSurvey -> deduplicate -> split -> quality sweep.\n\nRun this once. Everything downstream (training, calibration, the demo app)\nreads the CSVs it writes. Nothing here needs a GPU.\n\nOutputs:\n  splits.csv           id_code, diagnosis, split      <- the only split you ever use\n  duplicate_groups.csv id_code, group_id              <- audit trail\n  quality_metrics.csv  id_code, verdict, focus, ...   <- feeds the quality gate\n  thresholds.json      calibrated quality cut points\n\"\"\"\n\nimport json\nimport os\nfrom concurrent.futures import ProcessPoolExecutor\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nfrom sklearn.model_selection import train_test_split\n\nfrom fundus_quality import assess_quality, calibrate_thresholds\nfrom fundus_dedup import find_duplicate_groups\n\nROOT = \"/kaggle/input/competitions/aptos2019-blindness-detection\"\nOUT = \"/kaggle/working\"\nSEED = 42\nN_WORKERS = max(1, (os.cpu_count() or 2) - 1)\n\n\ndef load_index():\n    df = pd.read_csv(f\"{ROOT}/train.csv\")\n    df[\"path\"] = ROOT + \"/train_images/\" + df.id_code + \".png\"\n    return df\n\n\n# ---------------------------------------------------------------------------\n# 1. Resolution survey  (reads headers only -- no image decoding, so it is instant)\n# ---------------------------------------------------------------------------\n\ndef survey_resolutions(df):\n    dims = [Image.open(p).size for p in df.path]          # (w, h)\n    s = pd.DataFrame(dims, columns=[\"w\", \"h\"])\n    s[\"aspect\"] = (s.w / s.h).round(2)\n    s[\"megapixels\"] = (s.w * s.h / 1e6).round(2)\n\n    print(\"distinct (w,h) pairs:\", len(s.groupby(['w', 'h'])))\n    print(\"\\nmost common resolutions:\")\n    print(s.groupby([\"w\", \"h\"]).size().sort_values(ascending=False).head(8))\n    print(\"\\nshort side: min %d  median %d  max %d\"\n          % (s[[\"w\", \"h\"]].min(axis=1).min(),\n             s[[\"w\", \"h\"]].min(axis=1).median(),\n             s[[\"w\", \"h\"]].min(axis=1).max()))\n    return s\n\n\n# ---------------------------------------------------------------------------\n# 2. Near-duplicate detection via difference hash\n# ---------------------------------------------------------------------------\n\n# ---------------------------------------------------------------------------\n# 3. Group-aware stratified split\n# ---------------------------------------------------------------------------\n\ndef make_split(dupes, val_frac=0.15, test_frac=0.15, seed=SEED):\n    \"\"\"Split at the GROUP level so no image and its near-twin land on\n    opposite sides of the split. Stratify on the group's label.\"\"\"\n    g = (dupes.groupby(\"group_id\")\n               .agg(diagnosis=(\"diagnosis\", \"first\"), n=(\"id_code\", \"size\"))\n               .reset_index())\n\n    train_g, hold_g = train_test_split(\n        g, test_size=val_frac + test_frac,\n        stratify=g.diagnosis, random_state=seed)\n    val_g, test_g = train_test_split(\n        hold_g, test_size=test_frac / (val_frac + test_frac),\n        stratify=hold_g.diagnosis, random_state=seed)\n\n    assign = {}\n    for name, part in [(\"train\", train_g), (\"val\", val_g), (\"test\", test_g)]:\n        for gid in part.group_id:\n            assign[gid] = name\n\n    out = dupes.copy()\n    out[\"split\"] = out.group_id.map(assign)\n\n    print(\"\\nsplit sizes:\\n\", out.split.value_counts())\n    print(\"\\nclass balance per split (%):\")\n    print((pd.crosstab(out.split, out.diagnosis, normalize=\"index\") * 100).round(1))\n    return out[[\"id_code\", \"diagnosis\", \"split\"]]\n\n\n# ---------------------------------------------------------------------------\n# 4. Quality sweep\n# ---------------------------------------------------------------------------\n\ndef _quality_row(args):\n    id_code, path = args\n    img = cv2.imread(path)\n    if img is None:\n        return {\"id_code\": id_code, \"verdict\": \"reject\", \"reason\": \"unreadable\"}\n    r = assess_quality(img)\n    return {\"id_code\": id_code, \"verdict\": r.verdict,\n            \"reason\": r.reasons[0] if r.reasons else \"\", **r.metrics}\n\n\ndef quality_sweep(df):\n    args = list(zip(df.id_code, df.path))\n    with ProcessPoolExecutor(N_WORKERS) as ex:\n        rows = list(ex.map(_quality_row, args, chunksize=16))\n    q = pd.DataFrame(rows)\n    print(\"\\nverdicts with placeholder thresholds:\\n\", q.verdict.value_counts())\n    return q\n\n\ndef main():\n    df = load_index()\n    print(f\"{len(df)} labelled images\\n\")\n\n    print(\"=== 1. resolution survey ===\")\n    survey_resolutions(df)\n\n    print(\"\\n=== 2. duplicate detection ===\")\n    dupes = find_duplicate_groups(df, workers=N_WORKERS)\n    dupes.to_csv(f\"{OUT}/duplicate_groups.csv\", index=False)\n\n    print(\"\\n=== 3. split ===\")\n    splits = make_split(dupes)\n    splits.to_csv(f\"{OUT}/splits.csv\", index=False)\n\n    print(\"\\n=== 4. quality sweep ===\")\n    q = quality_sweep(df)\n    q.to_csv(f\"{OUT}/quality_metrics.csv\", index=False)\n\n    # calibrate on TRAIN ONLY -- calibrating on val/test leaks information\n    train_ids = set(splits.loc[splits.split == \"train\", \"id_code\"])\n    train_metrics = q[q.id_code.isin(train_ids)].dropna(subset=[\"focus\"])\n    cal = calibrate_thresholds(train_metrics.to_dict(\"records\"))\n    with open(f\"{OUT}/thresholds.json\", \"w\") as f:\n        json.dump(cal, f, indent=2)\n    print(\"\\ncalibrated thresholds (from train split only):\")\n    print(json.dumps(cal, indent=2))\n\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T08:54:45.360137Z","iopub.execute_input":"2026-08-30T08:54:45.360540Z","iopub.status.idle":"2026-08-30T09:02:37.092793Z","shell.execute_reply.started":"2026-08-30T08:54:45.360511Z","shell.execute_reply":"2026-08-30T09:02:37.091854Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd, inspect, fundus_quality\nprint(fundus_quality.__file__)\nprint(\"clip_h\" in inspect.getsource(fundus_quality.framing_scores))\nprint(pd.read_csv(\"/kaggle/working/quality_metrics.csv\").columns.tolist())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:02:37.093924Z","iopub.execute_input":"2026-08-30T09:02:37.094212Z","iopub.status.idle":"2026-08-30T09:02:37.109169Z","shell.execute_reply.started":"2026-08-30T09:02:37.094172Z","shell.execute_reply":"2026-08-30T09:02:37.108347Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"q = pd.read_csv(f\"{OUT}/quality_metrics.csv\")\nprint(q[q.verdict == \"reject\"].reason.value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:02:37.110468Z","iopub.execute_input":"2026-08-30T09:02:37.111217Z","iopub.status.idle":"2026-08-30T09:02:37.125923Z","shell.execute_reply.started":"2026-08-30T09:02:37.111180Z","shell.execute_reply":"2026-08-30T09:02:37.125078Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"q = pd.read_csv(\"/kaggle/working/quality_metrics.csv\")\nprint(q[q.verdict==\"reject\"].reason.value_counts())\n\n# 12 worst by focus, and 12 sitting right at the boundary\nimport matplotlib.pyplot as plt\nfor name, sub in [(\"worst\", q.nsmallest(12,\"focus\")),\n                  (\"boundary\", q.iloc[(q.focus-85.8).abs().argsort()[:12]])]:\n    fig, ax = plt.subplots(3,4, figsize=(16,12))\n    for a, (_, r) in zip(ax.flat, sub.iterrows()):\n        a.imshow(cv2.cvtColor(cv2.imread(f\"{ROOT}/train_images/{r.id_code}.png\",\n                 cv2.IMREAD_REDUCED_COLOR_4), cv2.COLOR_BGR2RGB))\n        a.set_title(f\"{r.verdict} f={r.focus:.0f}\", fontsize=9); a.axis(\"off\")\n    fig.suptitle(name); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:18:28.862193Z","iopub.execute_input":"2026-08-30T09:18:28.863127Z","iopub.status.idle":"2026-08-30T09:18:34.000143Z","shell.execute_reply.started":"2026-08-30T09:18:28.863035Z","shell.execute_reply":"2026-08-30T09:18:33.998986Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nWider boundary visual check.\n\nPulls a larger, focus-stratified sample around the focus_reject threshold\nso you can eyeball whether the cut point is in the right place, not just\nwhether individual borderline images look ambiguous (they will, by\ndefinition, near a boundary).\n\nUsage: paste into a notebook cell after quality_metrics.csv has been\nwritten (i.e. after step2_prepare.py has run).\n\"\"\"\n\nimport pandas as pd\nimport cv2\nimport matplotlib.pyplot as plt\n\nOUT = \"/kaggle/working\"\nROOT = \"/kaggle/input/competitions/aptos2019-blindness-detection\"\n\nq = pd.read_csv(f\"{OUT}/quality_metrics.csv\")\ndf = pd.read_csv(f\"{ROOT}/train.csv\")\ndf[\"path\"] = ROOT + \"/train_images/\" + df.id_code + \".png\"\ndf_path = df.set_index(\"id_code\")\n\n# widen the window around the current focus_reject threshold (85.78)\n# so you get real spread on both sides, not just images already near it\nLOW, HIGH = 60, 110\nband = q[(q.focus >= LOW) & (q.focus <= HIGH)].copy()\nband = band.sort_values(\"focus\").reset_index(drop=True)\n\nprint(f\"{len(band)} images with focus in [{LOW}, {HIGH}]\")\nprint(band.verdict.value_counts())\n\n# stratified sample: take images spread evenly across the sorted focus\n# range rather than randomly, so you see the full gradient low->high\nN = 24\nif len(band) > N:\n    # N evenly-spaced positions across [0, len(band)-1]\n    idx = (pd.Series(range(N)) * (len(band) - 1) / (N - 1)).round().astype(int).unique()\n    sample = band.iloc[idx].reset_index(drop=True)\nelse:\n    sample = band\n\nn = len(sample)\ncols = 4\nrows = (n + cols - 1) // cols\nfig, axes = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows))\naxes = axes.flatten()\n\nfor ax, (_, row) in zip(axes, sample.iterrows()):\n    img = cv2.cvtColor(cv2.imread(df_path.loc[row.id_code, \"path\"]), cv2.COLOR_BGR2RGB)\n    ax.imshow(img)\n    ax.set_title(f\"{row.verdict}  focus={row.focus:.0f}\", fontsize=10)\n    ax.axis(\"off\")\n\nfor ax in axes[n:]:\n    ax.axis(\"off\")\n\nplt.suptitle(f\"focus band [{LOW}-{HIGH}], sorted low -> high, current focus_reject=85.78\", y=1.0)\nplt.tight_layout()\nplt.savefig(f\"{OUT}/boundary_wide.png\", dpi=110, bbox_inches=\"tight\")\nplt.show()\n\n# tally sheet to fill in by eye as you look at the grid above\ntally = sample[[\"id_code\", \"verdict\", \"focus\"]].copy()\ntally[\"my_call_gradeable\"] = \"\"  # fill in: y / n / borderline, for each row as you review\ntally.to_csv(f\"{OUT}/boundary_tally.csv\", index=False)\nprint(f\"\\nFill in 'my_call_gradeable' in {OUT}/boundary_tally.csv (y/n/borderline) as you review,\")\nprint(\"then compare against the algorithm's verdict column.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:34:51.044111Z","iopub.execute_input":"2026-08-30T09:34:51.044474Z","iopub.status.idle":"2026-08-30T09:35:13.874815Z","shell.execute_reply.started":"2026-08-30T09:34:51.044442Z","shell.execute_reply":"2026-08-30T09:35:13.873864Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile dr_ordinal.py\n\"\"\"\ndr_ordinal.py\nOrdinal grading head + referable-DR metrics.\n\nWhy ordinal instead of 5-way softmax\n------------------------------------\nA softmax over {0,1,2,3,4} treats the grades as unrelated categories, so\nconfusing grade 0 with grade 4 costs exactly what confusing grade 3 with\ngrade 4 costs. That is wrong for DR: the scale is ordered, and the clinically\nimportant boundary is a threshold on that order, not a category.\n\nInstead we predict K-1 = 4 independent cumulative probabilities:\n\n    y1 = P(grade >= 1)   any DR\n    y2 = P(grade >= 2)   REFERABLE DR   <- the number the PS grades you on\n    y3 = P(grade >= 3)   severe\n    y4 = P(grade >= 4)   proliferative\n\nThis buys three things:\n  1. Ordinality is built into the target encoding.\n  2. Referable DR is a single output you can threshold directly, so hitting\n     \">90% sensitivity, >85% specificity\" becomes an explicit operating-point\n     choice on one ROC curve rather than a side effect of argmax.\n  3. That same output is the thing you temperature-calibrate for the\n     confidence score the explainability module needs.\n\"\"\"\n\nimport numpy as np\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nNUM_GRADES = 5\nREFERABLE_IDX = 1          # y2 above: P(grade >= 2)\n\n\n# ---------------------------------------------------------------------------\n# Target encoding / decoding\n# ---------------------------------------------------------------------------\n\ndef encode_ordinal(grades, num_grades=NUM_GRADES):\n    \"\"\"grade -> cumulative binary vector.  2 -> [1, 1, 0, 0]\"\"\"\n    g = torch.as_tensor(grades).long().view(-1, 1)\n    ks = torch.arange(1, num_grades, device=g.device).view(1, -1)\n    return (g >= ks).float()\n\n\ndef decode_ordinal(probs, threshold=0.5):\n    \"\"\"Cumulative probabilities -> integer grade.\n\n    Uses the count of exceeded thresholds rather than a first-crossing scan.\n    The head can emit non-monotonic outputs (e.g. [0.9, 0.3, 0.6, 0.1]); a\n    scan would stop at the first dip and silently under-grade, while counting\n    degrades gracefully.\n    \"\"\"\n    p = torch.as_tensor(probs)\n    return (p > threshold).sum(dim=1).clamp(0, NUM_GRADES - 1)\n\n\nclass OrdinalHead(nn.Module):\n    def __init__(self, in_features, num_grades=NUM_GRADES):\n        super().__init__()\n        self.fc = nn.Linear(in_features, num_grades - 1)\n\n    def forward(self, x):\n        return self.fc(x)                      # logits, (B, 4)\n\n\ndef ordinal_loss(logits, grades, pos_weight=None):\n    \"\"\"Binary cross-entropy over the cumulative targets.\n\n    pos_weight (shape (4,)) counteracts the imbalance at each threshold.\n    The grade>=4 threshold has only ~8% positives, so without weighting the\n    model learns to always say no there.\n    \"\"\"\n    target = encode_ordinal(grades).to(logits.device)\n    return F.binary_cross_entropy_with_logits(\n        logits, target, pos_weight=pos_weight)\n\n\ndef ordinal_pos_weight(train_grades, num_grades=NUM_GRADES):\n    \"\"\"neg/pos ratio at each cumulative threshold, from the training split.\"\"\"\n    g = np.asarray(train_grades)\n    w = []\n    for k in range(1, num_grades):\n        pos = float((g >= k).sum())\n        neg = float((g < k).sum())\n        w.append(neg / max(pos, 1.0))\n    return torch.tensor(w, dtype=torch.float32)\n\n\n# ---------------------------------------------------------------------------\n# Metrics\n# ---------------------------------------------------------------------------\n\ndef quadratic_weighted_kappa(y_true, y_pred, num_grades=NUM_GRADES):\n    y_true = np.asarray(y_true, int)\n    y_pred = np.asarray(y_pred, int)\n    O = np.zeros((num_grades, num_grades))\n    for t, p in zip(y_true, y_pred):\n        O[t, p] += 1\n\n    w = np.zeros((num_grades, num_grades))\n    for i in range(num_grades):\n        for j in range(num_grades):\n            w[i, j] = ((i - j) ** 2) / ((num_grades - 1) ** 2)\n\n    hist_t = np.bincount(y_true, minlength=num_grades)\n    hist_p = np.bincount(y_pred, minlength=num_grades)\n    E = np.outer(hist_t, hist_p).astype(float)\n    E *= O.sum() / E.sum()\n\n    denom = (w * E).sum()\n    return 1.0 - (w * O).sum() / denom if denom > 0 else 0.0\n\n\ndef referable_operating_point(y_true_grade, referable_prob,\n                              min_sensitivity=0.90):\n    \"\"\"Pick the threshold meeting the sensitivity floor, then maximise specificity.\n\n    This is the crux of the problem statement. Sensitivity and specificity are\n    NOT properties of the model -- they are a point you choose on its ROC\n    curve. Reporting whatever argmax happens to give is leaving the PS target\n    to chance.\n    \"\"\"\n    y = (np.asarray(y_true_grade) >= 2).astype(int)\n    p = np.asarray(referable_prob, dtype=float)\n\n    best = None\n    for thr in np.unique(np.round(p, 4)):\n        pred = (p >= thr).astype(int)\n        tp = int(((pred == 1) & (y == 1)).sum())\n        fn = int(((pred == 0) & (y == 1)).sum())\n        tn = int(((pred == 0) & (y == 0)).sum())\n        fp = int(((pred == 1) & (y == 0)).sum())\n        sens = tp / max(tp + fn, 1)\n        spec = tn / max(tn + fp, 1)\n        if sens >= min_sensitivity and (best is None or spec > best[\"specificity\"]):\n            best = dict(threshold=float(thr), sensitivity=sens, specificity=spec,\n                        tp=tp, fp=fp, tn=tn, fn=fn)\n\n    if best is None:\n        return dict(threshold=float(p.min()), sensitivity=1.0, specificity=0.0,\n                    note=\"sensitivity floor unreachable at any threshold\")\n    return best\n\n\ndef roc_auc(y_true_binary, scores):\n    \"\"\"AUC via the rank identity (equivalent to Mann-Whitney U).\"\"\"\n    y = np.asarray(y_true_binary, int)\n    s = np.asarray(scores, float)\n    n_pos, n_neg = int(y.sum()), int((1 - y).sum())\n    if n_pos == 0 or n_neg == 0:\n        return float(\"nan\")\n    order = s.argsort()\n    ranks = np.empty(len(s), float)\n    ranks[order] = np.arange(1, len(s) + 1)\n    # average ranks within ties\n    _, inv, counts = np.unique(s, return_inverse=True, return_counts=True)\n    sums = np.zeros(len(counts))\n    np.add.at(sums, inv, ranks)\n    ranks = (sums / counts)[inv]\n    return (ranks[y == 1].sum() - n_pos * (n_pos + 1) / 2) / (n_pos * n_neg)\n\n\n# ---------------------------------------------------------------------------\n# Confidence calibration\n# ---------------------------------------------------------------------------\n\nclass TemperatureScaler(nn.Module):\n    \"\"\"One scalar per output, fitted on the VALIDATION split after training.\n\n    Raw network probabilities are systematically overconfident. The PS asks\n    for 'calibrated confidence scores', and this is the standard, cheap,\n    defensible way to get them. Fit on val, report ECE on test.\n    \"\"\"\n\n    def __init__(self, n_outputs=NUM_GRADES - 1):\n        super().__init__()\n        self.log_t = nn.Parameter(torch.zeros(n_outputs))\n\n    def forward(self, logits):\n        return logits / torch.exp(self.log_t)\n\n    def fit(self, logits, grades, iters=300, lr=0.02):\n        logits = torch.as_tensor(logits, dtype=torch.float32)\n        target = encode_ordinal(grades)\n        opt = torch.optim.LBFGS([self.log_t], lr=lr, max_iter=iters)\n\n        def closure():\n            opt.zero_grad()\n            loss = F.binary_cross_entropy_with_logits(self(logits), target)\n            loss.backward()\n            return loss\n\n        opt.step(closure)\n        return self\n\n\ndef expected_calibration_error(probs, labels, n_bins=15):\n    \"\"\"ECE for one binary output. Report this next to a reliability diagram.\"\"\"\n    p = np.asarray(probs, float)\n    y = np.asarray(labels, int)\n    edges = np.linspace(0, 1, n_bins + 1)\n    ece = 0.0\n    for lo, hi in zip(edges[:-1], edges[1:]):\n        m = (p > lo) & (p <= hi)\n        if m.sum():\n            ece += (m.sum() / len(p)) * abs(y[m].mean() - p[m].mean())\n    return float(ece)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:51:55.353871Z","iopub.execute_input":"2026-08-30T09:51:55.354224Z","iopub.status.idle":"2026-08-30T09:51:55.363219Z","shell.execute_reply.started":"2026-08-30T09:51:55.354193Z","shell.execute_reply":"2026-08-30T09:51:55.362326Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"q = pd.read_csv(\"/kaggle/working/quality_metrics.csv\")\ns = pd.read_csv(\"/kaggle/working/splits.csv\")\nm = q.merge(s, on=\"id_code\")\nprint(m.groupby(\"diagnosis\").focus.describe()[[\"25%\",\"50%\",\"75%\"]])\nprint(pd.crosstab(m.diagnosis, m.verdict, normalize=\"index\").round(3))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:53:30.530358Z","iopub.execute_input":"2026-08-30T09:53:30.531345Z","iopub.status.idle":"2026-08-30T09:53:30.590677Z","shell.execute_reply.started":"2026-08-30T09:53:30.531297Z","shell.execute_reply":"2026-08-30T09:53:30.589873Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"m[\"res\"] = m.id_code.map(lambda i: Image.open(f\"{ROOT}/train_images/{i}.png\").size)\nprint(m.groupby([\"res\",\"diagnosis\"]).focus.median().unstack().round(0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T10:09:27.616709Z","iopub.execute_input":"2026-08-30T10:09:27.619359Z","iopub.status.idle":"2026-08-30T10:09:34.139729Z","shell.execute_reply.started":"2026-08-30T10:09:27.619160Z","shell.execute_reply":"2026-08-30T10:09:34.138717Z"}},"outputs":[],"execution_count":null}]}