{"cells":[{"cell_type":"markdown","metadata":{},"source":"# Knee MRI, End to End: EDA → DICOM Preprocessing → 2.5D Models → Honest Evaluation 🦵🧲\n\nHow I built a knee-MRI abnormality pipeline for the RSNA competition — every stage, with the real numbers and the real mistakes:\n\n1. **EDA** — what the data actually is: studies, series, planes, and a label set where only 58 studies are fully annotated\n2. **DICOM preprocessing** — windowing, plane slot-maps, and building 2.5D tensors\n3. **Four model architectures** — from a plain 2.5D ResNet to a plane-gated mixture of experts (full code)\n4. **Evaluation that doesn't lie** — why the 58 annotated studies misled me, measured with bootstrap CIs\n5. **Ensembling** — rank-mean across members, and what it bought\n\nMy four-member ensemble scores **0.919** on the leaderboard by itself; blended with the public community ensembles it reaches **0.937**. Everything here runs on CPU in minutes — the heavy training happened offline, and the member predictions needed for Part 4 are attached as a small public dataset."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"import os, glob\nimport numpy as np, pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pydicom, cv2\nfrom scipy.stats import rankdata\nfrom sklearn.metrics import roc_auc_score\nimport warnings; warnings.filterwarnings('ignore')\nsns.set_style('whitegrid')\n\nROOT = os.path.dirname(glob.glob('/kaggle/input/**/train_series.csv', recursive=True)[0])\ntrain = pd.read_csv(f'{ROOT}/train.csv')\ntser = pd.read_csv(f'{ROOT}/train_series.csv')\nLABELS = ['ACL', 'MCL', 'Medial Meniscus', 'Lateral Meniscus', 'Medial OA', 'Lateral OA',\n          'PF OA', 'Effusion', 'Synovitis', \"Baker's\", 'Contusion', 'Fracture']\nprint(f'{train.shape[0]:,} studies | {tser.shape[0]:,} series')"},{"cell_type":"markdown","metadata":{},"source":"## 1. EDA — the shape of a medical imaging dataset\n\nA *study* is one patient visit; each study contains several *series* (one MRI sweep each: a plane + a pulse sequence); each series is a stack of DICOM slices. The 12 targets are knee findings, scored by **macro AUC**.\n\nThe first thing to understand is the label situation — it defines everything downstream."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"ann = train[train[LABELS].notna().all(axis=1)]\nprint(f'studies fully annotated by radiologists: {len(ann)} of {len(train)} ({len(ann)/len(train):.1%})')\nfig, ax = plt.subplots(figsize=(10, 3))\navail = train[LABELS].notna().mean().sort_values()\navail.plot(kind='barh', ax=ax, color='#2a78d6')\nax.set_title('fraction of studies with a label, per finding'); plt.tight_layout(); plt.show()"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"prev = ann[LABELS].mean().sort_values()\nfig, axes = plt.subplots(1, 2, figsize=(13, 3.5))\nprev.plot(kind='barh', ax=axes[0], color='#2a78d6')\naxes[0].set_title('finding prevalence (58 annotated studies)')\nsns.heatmap(ann[LABELS].corr(), ax=axes[1], cmap='RdBu_r', center=0, cbar_kws={'shrink': .7})\naxes[1].set_title('finding co-occurrence'); plt.tight_layout(); plt.show()"},{"cell_type":"markdown","metadata":{},"source":"**Observations**\n- Only **58 studies** carry radiologist labels — the rest have labels *derivable from report text* (I measured which derivations to trust in a [separate notebook](https://www.kaggle.com/code/starkhushi/rsna-knee-which-report-labels-should-you-train-on)).\n- Findings co-occur (OA findings cluster; effusion travels with synovitis) — a multi-label head shares that signal.\n- With 58 labeled studies and findings having as few as 9 positives, **evaluation noise will dominate model differences** — measured properly in Part 4."},{"cell_type":"markdown","metadata":{},"source":"## 1b. What the images look like — one study, three planes"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"tser['plane'] = tser['Anatomical_Plane']\nprint(tser.groupby(['plane'])['SeriesInstanceUID'].count().rename('series count'))\nuid = ann['StudyInstanceUID'].iloc[0]\nsmap = tser[tser.StudyInstanceUID == uid]\nfig, axes = plt.subplots(1, 3, figsize=(13, 4.5))\nfor ax, plane in zip(axes, ['Sagittal', 'Coronal', 'Axial']):\n    row = smap[smap.plane == plane]\n    if len(row) == 0:\n        ax.axis('off'); continue\n    sd = f\"{ROOT}/train_series/{uid}/{row.SeriesInstanceUID.iloc[0]}\"\n    fs = sorted(glob.glob(sd + '/*.dcm'))\n    ds = pydicom.dcmread(fs[len(fs)//2])\n    ax.imshow(ds.pixel_array, cmap='gray'); ax.set_title(f'{plane} (middle slice)'); ax.axis('off')\nplt.tight_layout(); plt.show()"},{"cell_type":"markdown","metadata":{},"source":"## 2. Preprocessing — from DICOM to model tensor\n\nThree decisions carry all the weight:\n\n1. **Windowing.** Raw MRI intensities vary wildly between scanners. Percentile windowing (clip at the 1st/99th percentile, scale to [0,1]) makes every slice comparable — no fixed constants that break on outlier scanners.\n2. **Plane slot-maps.** Each study gets 3 fixed slots: *fluid-sensitive sagittal*, *coronal*, *axial*. Fluid-sensitive sequences light up exactly the pathology (effusion, edema) we're asked to find.\n3. **2.5D stacks.** 16 evenly-spaced slices per slot, resized to 224px → a `(3 slots, 16 slices, 224, 224)` tensor. 3D context at 2D compute cost."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def window(img):\n    lo, hi = np.percentile(img, [1, 99])\n    return np.clip((img - lo) / max(hi - lo, 1e-6), 0, 1).astype(np.float32)\n\nraw = pydicom.dcmread(fs[len(fs)//2]).pixel_array.astype(np.float32)\nfig, axes = plt.subplots(1, 3, figsize=(12, 3.6))\naxes[0].hist(raw.ravel(), bins=80, color='#2a78d6'); axes[0].set_title('raw intensity histogram')\naxes[1].imshow(raw, cmap='gray'); axes[1].set_title('raw'); axes[1].axis('off')\naxes[2].imshow(window(raw), cmap='gray'); axes[2].set_title('percentile-windowed'); axes[2].axis('off')\nplt.tight_layout(); plt.show()"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"SLOTS = [('sagittal fluid-sensitive', (tser.plane == 'Sagittal') & (tser.Fluid_Sensitive == 1)),\n         ('coronal', tser.plane == 'Coronal'),\n         ('axial', tser.plane == 'Axial')]\nfor name, mask in SLOTS:\n    cover = tser[mask].StudyInstanceUID.nunique() / train.shape[0]\n    print(f'slot {name:28s} available in {cover:.1%} of studies')\nprint('-> missing slots get a zero tensor and a mask flag; the models learn to route around them')"},{"cell_type":"markdown","metadata":{},"source":"## 3. The four architectures\n\nTrained offline (each ~2–4 GPU-hours). Full definitions below — the *differences between them* are what make the ensemble work later.\n\n| member | backbone | pooling idea |\n|---|---|---|\n| v2 | ResNet-34, single sagittal stack | mean over slices — the honest baseline |\n| v3 | ResNet-34, all 3 slots | **plane-gated mixture-of-experts**: each finding learns which plane to trust |\n| v4 | ConvNeXt-Tiny | gated **slice-attention** pooling, then plane gating |\n| v5 | ResNet-50 (RadImageNet init) | radiology-pretrained transfer + plane gating |"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"import torch, torch.nn as nn, torchvision\n\nclass NetV2(nn.Module):                       # plain 2.5D baseline\n    def __init__(self, n_out=12):\n        super().__init__()\n        bb = torchvision.models.resnet34(weights=None); bb.fc = nn.Identity()\n        self.bb = bb; self.head = nn.Linear(512, n_out)\n    def forward(self, x):                     # (B, S, H, W)\n        B, S, H, W = x.shape\n        f = self.bb(x.reshape(B*S, 1, H, W).repeat(1, 3, 1, 1)).reshape(B, S, -1).mean(1)\n        return self.head(f)\n\nclass PlaneMixNet(nn.Module):                 # v3: plane-gated mixture of experts\n    def __init__(self, n_slots=3, n_out=12):\n        super().__init__()\n        bb = torchvision.models.resnet34(weights=None); bb.fc = nn.Identity()\n        self.bb = bb\n        self.slot_emb = nn.Parameter(torch.zeros(n_slots, 512))\n        self.expert = nn.Linear(512, n_out); self.gate = nn.Linear(512, n_out)\n    def forward(self, x, mask):               # (B, slots, slices, H, W)\n        B, S, Csl, H, W = x.shape\n        f = self.bb(x.reshape(B*S*Csl, 1, H, W).repeat(1, 3, 1, 1)).reshape(B, S, Csl, -1).mean(2) + self.slot_emb\n        scores = self.gate(f).masked_fill(mask.unsqueeze(-1) == 0, -1e4)\n        return (torch.softmax(scores, 1) * self.expert(f)).sum(1)\n\nprint('v2 params: %.1fM | v3 params: %.1fM' % (\n    sum(p.numel() for p in NetV2().parameters()) / 1e6,\n    sum(p.numel() for p in PlaneMixNet().parameters()) / 1e6))"},{"cell_type":"markdown","metadata":{},"source":"**Why plane gating matters:** ACL tears live in the sagittal plane, MCL in the coronal, patellofemoral OA in the axial. A per-finding softmax over plane slots lets each finding *learn* its own reading order — exactly what a radiologist does. (v4 and v5 refine the same skeleton: attention over slices instead of a mean, and RadImageNet initialization — code in my inference kernel.)"},{"cell_type":"markdown","metadata":{},"source":"## 4. Evaluation that doesn't lie\n\nI selected checkpoints on the 58 annotated studies. **They ranked my models wrong** — my \"best\" model there was my worst on the leaderboard. Here's the measurement that explains it, using the saved member predictions (attached dataset)."},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"z = np.load(glob.glob('/kaggle/input/**/member_val_preds.npz', recursive=True)[0], allow_pickle=True)\nY = z['Y']; names = ['V2', 'V3', 'V4', 'V5', 'V6', 'V7']\ndef macro(P, Y):\n    aucs = [roc_auc_score(Y[:, j], P[:, j]) for j in range(Y.shape[1]) if len(np.unique(Y[:, j])) > 1]\n    return float(np.mean(aucs))\ntbl = pd.DataFrame({n: {'macro AUC (58 studies)': macro(z[n], Y)} for n in names}).T\ntbl.round(4)"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"# bootstrap the 58 studies: how much does macro AUC move from sampling alone?\nrng = np.random.RandomState(0)\nboots = {n: [] for n in names}\nfor _ in range(400):\n    idx = rng.choice(len(Y), len(Y), replace=True)\n    for n in names:\n        try: boots[n].append(macro(z[n][idx], Y[idx]))\n        except ValueError: pass\nfig, ax = plt.subplots(figsize=(9, 3.5))\nax.boxplot([boots[n] for n in names], tick_labels=names)\nax.set_title('bootstrap distribution of macro AUC on the 58 annotated studies')\nax.set_ylabel('macro AUC'); plt.tight_layout(); plt.show()\nspread = np.mean([np.percentile(boots[n], 97.5) - np.percentile(boots[n], 2.5) for n in names])\nprint(f'average 95% CI width: {spread:.3f} -- larger than every between-model gap above')"},{"cell_type":"markdown","metadata":{},"source":"**The lesson:** the CI on 58 studies is wider than the differences between models. Any checkpoint choice made on this set is a coin flip. The fix that worked: select on a ~250-study holdout built from *derived* labels (imperfect labels, but n crushes the variance), and keep the 58 as a final report only."},{"cell_type":"markdown","metadata":{},"source":"## 5. Ensembling — rank-mean across members"},{"cell_type":"code","metadata":{},"execution_count":null,"outputs":[],"source":"def rank_norm(P):\n    return np.stack([rankdata(P[:, j]) / len(P) for j in range(P.shape[1])], 1)\nR = {n: rank_norm(z[n]) for n in names}\ncombos = {'best single': max((macro(z[n], Y), n) for n in names)[0],\n          'prob-mean (V2..V5)': macro(np.mean([z[n] for n in names[:4]], 0), Y),\n          'rank-mean (V2..V5)': macro(np.mean([R[n] for n in names[:4]], 0), Y),\n          'rank-mean (all 6)': macro(np.mean([R[n] for n in names], 0), Y)}\npd.Series(combos).round(4)"},{"cell_type":"markdown","metadata":{},"source":"**Rank-mean beats probability-mean** because AUC reads only the ordering: rank-normalizing first stops an over-confident member from dominating the average. On the leaderboard this ensemble scores **0.919**; rank-blending it with the community's public ensembles reaches **0.937** — diversity is the whole game.\n\n## Takeaways 🎯\n1. **Understand your labels before your models** — 58 gold labels + thousands of derivable ones dictated the whole validation design.\n2. **Windowing + plane slots + 2.5D** turns messy multi-scanner DICOM into a fixed tensor without losing the radiologist's structure.\n3. **Give each finding a choice of planes** (gated experts) — anatomy is not plane-agnostic.\n4. **Bootstrap your validation set before trusting it** — mine was lying, measurably.\n5. **Rank-mean ensembling of diverse members** is the most reliable free lunch in the sport.\n\nIf this helped, an **upvote** is appreciated 🙌 Related: [my validation-noise notebook](https://www.kaggle.com/code/starkhushi/58-study-validation-set-is-lying-to-you) and [the report-label study](https://www.kaggle.com/code/starkhushi/rsna-knee-which-report-labels-should-you-train-on)."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11"}},"nbformat":4,"nbformat_minor":5}