{"cells":[{"cell_type":"markdown","metadata":{},"source":"# Knee MRI: EDA to baseline, decision by decision\n\nRSNA Knee Abnormality Detection asks, for one knee MRI exam, for the probability\nthat each of **12 findings** is present. The score is the mean AUC over the 12\nfindings (macro AUC).\n\nThis notebook walks through what the data actually looks like, and how each\nobservation pins down a piece of the baseline. In this competition a surprising\nnumber of things get decided before you touch a model.\n\nWhat you'll get out of it:\n\n- where the real difficulty is (it isn't the image processing)\n- which provided columns you cannot use as-is\n- how to split validation so it means something\n- how all of the above fixes the baseline configuration\n\nAt the end I walk through the path to **public LB 0.898**, one change at a time.\n\n## 0. The shape of the problem\n\n| | |\n|---|---|\n| Input | One exam (study) = several imaging series (sagittal, coronal, axial, ...) |\n| Output | Probability for each of 12 findings |\n| Metric | Mean AUC over the 12 findings |\n| Format | Code Competition (submit a notebook, under 9 hours) |\n\nThe 12 findings:\n\n`ACL` `MCL` `Medial Meniscus` `Lateral Meniscus` `Medial OA` `Lateral OA`\n`PF OA` `Effusion` `Synovitis` `Baker's` `Contusion` `Fracture`\n\n```\n[ one exam (study) ]\n\n  images -+- Sagittal series -+\n          +- Coronal  series -+- each series = tens of slices\n          +- Axial    series -+   (cross-sections through the knee)\n\n  text ---- one radiology report written by a radiologist\n            (not necessarily English: 9 languages are mixed in -> section 2)\n\n                          |\n                          v  predict\n\n  +--------------------------------------------+\n  |  probability of each of the 12 findings     |\n  |  ACL / MCL / Medial Meniscus / ...          |\n  +--------------------------------------------+\n```\n\nThe same knee is imaged from three directions. Which finding shows up in which\ndirection is not arbitrary, and I measure it later (section 6).\n\nThe names are medical but they all describe \"which part of the knee is damaged,\nand how\". One example: the **`ACL` (anterior cruciate ligament) sits in the\nmiddle of the joint and connects the femur to the tibia**; it tears during sudden\nchanges of direction in sport. The label is whether that ligament is torn. The\nother eleven are the same idea - a **structure** (ligament, meniscus, cartilage,\nbone, joint fluid) crossed with a **kind of damage**. You don't need more anatomy\nthan that to follow the rest.\n\n### What the metric decides for you\n\nAUC only looks at **ranking**. \"0.9 vs 0.8\" and \"0.6 vs 0.5\" score the same as\nlong as the order is the same. Two consequences:\n\n1. **No calibration needed.** Any monotone transform of your predictions leaves\n   the score unchanged.\n2. **You can blend on ranks.** Models whose outputs live on different scales can\n   still be averaged fairly once converted to ranks.\n"},{"cell_type":"code","metadata":{},"source":"from pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nROOT = next(p for p in [Path(\"/kaggle/input/rsna-knee-abnormality-detection\"),\n                        Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")]\n            if p.exists())\n\ntrain = pd.read_csv(ROOT / \"train.csv\")\nseries = pd.read_csv(ROOT / \"train_series.csv\")\n\nTARGETS = [\"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\", \"Medial OA\",\n           \"Lateral OA\", \"PF OA\", \"Effusion\", \"Synovitis\", \"Baker's\",\n           \"Contusion\", \"Fracture\"]\n\nprint(f\"train   : {len(train):,} exams\")\nprint(f\"series  : {len(series):,} series\")\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 1. The thing to know first: only 58 exams have labels\n\nOpen `train.csv` and it's immediately visible.\n\n**train has 4,407 exams, but only 58 of them carry the 12 findings as labels.**\nThe other 4,349 come with the **radiologist's free-text report** and nothing else.\n\n```\ntrain 4,407 exams\n  |\n  +-    58 : images + report + [12 findings labelled]\n  |            -> ground truth, so useful for checking labels you build\n  |\n  +- 4,349 : images + report + (no labels)\n               |\n               -> unusable for training as-is.\n                  Read the report and they become usable.\n```\n\nTo train on those 4,349 exams you have to **read the report text and decide the\npresence/absence of each of the 12 findings yourself**. Skip that and your\ntraining set stays at 58 exams, which is not enough to train an image model.\n\nSo before touching \"predict findings from images\", this competition makes you\nsolve **\"build training labels from reports\"** first. Three steps:\n\n1. Read the reports and **build your own labels** for all 4,407 exams\n2. Train an image model on those labels\n3. Check label quality against the 58 ground-truth exams\n\nThe 58 exams are few, but they are the only honest check you have on step 1.\n"},{"cell_type":"code","metadata":{},"source":"# ---- how many exams actually carry labels ---------------------------------\nlabeled = train[TARGETS].notna().all(axis=1)\nprint(f\"exams with all 12 findings labelled: {labeled.sum()} / {len(train):,}\")\nprint(f\"exams with a report                : {train['Report'].notna().sum():,}\")\n\nprint(\"\\n--- example report (first 200 chars) ---\")\nprint(train.loc[train[\"Report\"].notna(), \"Report\"].iloc[0][:200])\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 2. The reports are in nine languages\n\nNot just English. English, Spanish, French, Dutch, German, Turkish, Croatian,\nGreek and Russian are all present.\n\nThe design decision here is **not to detect the language and branch on it**.\n\nThe reason is simple: **if the language detection is wrong, that report becomes\nunreadable entirely**. Cheap heuristics like \"contains 'the', so English\" fire\nhappily on Dutch and German text too.\n\nInstead I build **one dictionary holding the cues for all nine languages** and\napply all of it to every report. There is no language to get wrong, so that\nwhole failure mode disappears.\n\nReading a report needs three kinds of expression:\n\n| kind | examples |\n|---|---|\n| the finding itself | `anterior cruciate`, `cruzado anterior`, `vorderes kreuzband` |\n| negation | `no`, `sin`, `ohne`, `unremarkable` |\n| hedging | `possible`, `suspected`, `cannot exclude` |\n\nMiss the negation and hedging and \"no ACL tear\" becomes \"ACL tear present\".\n"},{"cell_type":"code","metadata":{},"source":"# ---- which languages show up ----------------------------------------------\n# No language detection (section 2). This is only to confirm it isn't all English.\nCUE = {\n    \"English\": r\"\\bmeniscus\\b|\\bligament\\b\",\n    \"Spanish\": r\"\\bmenisco\\b|\\bligamento\\b\",\n    \"French\":  r\"\\bm[eé]nisque\\b|ligament crois\",\n    \"Dutch\":   r\"\\bkruisband\\b\",\n    \"German\":  r\"\\bmeniskus\\b|\\bkreuzband\\b\",\n    \"Turkish\": r\"men[iı]sk[uü]s|çapraz bağ\",\n    \"Russian\": r\"мениск|связк\",\n    \"Greek\":   r\"μηνίσκ|σύνδεσμ\",\n}\nrep = train[\"Report\"].fillna(\"\").str.lower()\nfor name, pat in CUE.items():\n    n = rep.str.contains(pat, regex=True).sum()\n    print(f\"  reports containing {name:8s} cues: {n:5,}\")\nprint(\"\\n-> Not English-only. Apply all nine dictionaries at once, no detection.\")\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 3. The image data, in numbers\n\n| | |\n|---|---|\n| DICOM files | **819,078** |\n| Total size | **~490 GB** |\n| Median file size | 0.526 MB |\n\n0.526 MB matches exactly what a 512x512 16-bit image takes **uncompressed**. And\nindeed, all 22,859 files I sampled were uncompressed (Explicit VR Little Endian).\n\nSo decoding is cheap - there is no compression to undo. The expensive parts are\n**opening a very large number of files** and **resizing**.\n\nWhich means: re-reading 490 GB every epoch is out of the question. Run the\npreprocessing **once** and cache small arrays; training then only ever touches\nthe cache.\n\n## 4. Check what the provided columns actually contain\n\n`train_series.csv` describes each series with `Anatomical_Plane`,\n`Fluid_Sensitive` and `Fat_Suppression`.\n\nBy their names, `Fluid_Sensitive` (does this sequence highlight fluid?) and\n`Fat_Suppression` (does it null out fat?) look like **different properties**.\nPhysically they are independent. But count them:\n\n**The two columns agree on all 24,371 series (100%).** One is a copy of the\nother, and neither of them tells you the contrast (T1 / PD / T2).\n\nContrast matters - it's what makes a finding visible at all:\n\n| | contrast | what it shows |\n|---|---|---|\n| TR < 800 | T1 | fat bright. anatomy |\n| TE < 50 | PD | intermediate. meniscus |\n| otherwise | T2 | fluid bright. effusion, inflammation |\n\nYou can also read it off the series description text, but **11.8% are anonymised\nto `DummySeriesDesc!`**, so text alone classifies only 77%. TE/TR is present on\n95.1% of series, so **using both gets you to 100%**.\n\n> The 11.8% / 95.1% figures come from scanning the headers of all 24,371 series\n> once. The cell below samples 300 series to keep the runtime down, so the\n> percentages it prints will differ somewhat.\n\nA provided column does not necessarily contain what its name suggests. Counting\ntwo of them against each other is what surfaced both facts here: that one\nduplicates the other, and that neither encodes T1/PD/T2. Used as features on the\nstrength of their names, neither would have come up.\n"},{"cell_type":"code","metadata":{},"source":"# ---- do the provided columns describe different properties? ---------------\nct = pd.crosstab(series[\"Fluid_Sensitive\"], series[\"Fat_Suppression\"])\nprint(\"Fluid_Sensitive x Fat_Suppression:\")\nprint(ct)\n\nagree = (series[\"Fluid_Sensitive\"] == series[\"Fat_Suppression\"]).mean()\nprint(f\"\\nfraction where the two agree: {agree:.1%}\")\nprint(\"-> They say the same thing. Neither carries contrast (T1/PD/T2).\")\nprint(\"   Recover it from DICOM TR/TE.\")\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- peek at DICOM headers (physical scale and contrast) ------------------\n# Opening all 819,078 files is expensive, so take one file per series for a\n# random 300 series. That is plenty to see the distribution.\nimport pydicom\n\nsample = series.sample(min(300, len(series)), random_state=0)\n\nrows = []\nfor r in sample.itertuples():\n    d = ROOT / \"train_series\" / r.StudyInstanceUID / r.SeriesInstanceUID\n    files = list(d.glob(\"*.dcm\"))[:1] if d.exists() else []\n    if not files:\n        continue\n    try:\n        ds = pydicom.dcmread(files[0], stop_before_pixels=True)\n    except Exception:\n        continue\n    ps = getattr(ds, \"PixelSpacing\", None)\n    rows.append(dict(\n        spacing=float(ps[0]) if ps is not None else np.nan,\n        rows_px=float(getattr(ds, \"Rows\", np.nan)),\n        TR=float(getattr(ds, \"RepetitionTime\", np.nan) or np.nan),\n        TE=float(getattr(ds, \"EchoTime\", np.nan) or np.nan),\n        desc=str(getattr(ds, \"SeriesDescription\", \"\")),\n    ))\n\nhdr = pd.DataFrame(rows)\nprint(f\"series read: {len(hdr)}\")\nprint(\"\\nmillimetres per pixel:\")\nprint(hdr[\"spacing\"].describe()[[\"min\", \"25%\", \"50%\", \"75%\", \"max\"]].round(3))\nprint(\"\\n-> It varies by site and scanner. Resizing alone will not align sizes.\")\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- recover contrast from TR/TE ------------------------------------------\ndef contrast_of(tr, te):\n    \"\"\"Map a TR/TE pair to T1 / PD / T2. A simple rule is good enough.\"\"\"\n    if not np.isfinite(tr) or not np.isfinite(te):\n        return \"unknown\"\n    if tr < 800:\n        return \"T1\"\n    if te < 50:\n        return \"PD\"\n    return \"T2\"\n\n\nhdr[\"contrast\"] = [contrast_of(t, e) for t, e in zip(hdr[\"TR\"], hdr[\"TE\"])]\nprint(\"contrast from TR/TE:\")\nprint(hdr[\"contrast\"].value_counts().to_string())\n\ndummy = hdr[\"desc\"].str.contains(\"Dummy\", case=False, na=False).mean()\nprint(f\"\\nfraction of series with an anonymised description: {dummy:.1%}\")\nprint(\"-> Names alone cannot classify them. Combine with TR/TE.\")\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 5. Aligning physical scale - at 160mm it was doing nothing 28.9% of the time\n\nIn MRI, **how many millimetres one pixel covers differs between sites and\nscanners**. DICOM stores it in `PixelSpacing`.\n\nIgnore it and \"resize everything to 224x224\" means **the same meniscus appears at\nmore than twice the size** in one image versus another. To the model those are\ndifferent objects.\n\nSo: fix a **physical field of view in millimetres** around the joint, crop that,\nand resize the crop.\n\n```\ncrop width (pixels) = field of view (mm) / PixelSpacing (mm/pixel)\n```\n\nI started with a 160mm field of view. It **was not working**.\n\nSome acquisitions are narrower than 160mm to begin with. When that happens the\ncrop width saturates at the image size, and **the physical-scale alignment\nsilently does nothing** - you are just using the whole image.\n\nCounted: at 160mm this happened on **28.9% of series**. At 130mm it drops to\n0.5%. (Both from scanning every series; the cell below samples 300, so its\nnumbers differ slightly.)\n\nFixing only this one thing moved the LB from 0.857 to 0.861.\n\nWorth stating plainly: \"implemented\" and \"having an effect\" are different\nthings. The way to tell them apart is to **count the cases where it didn't\napply**.\n"},{"cell_type":"code","metadata":{},"source":"# ---- what field of view should we use? ------------------------------------\n# crop width = fov / spacing. Above the image size it saturates and the\n# physical-scale alignment stops doing anything. Count how often, per fov.\nok = hdr.dropna(subset=[\"spacing\", \"rows_px\"])\nprint(f\"{'fov':>6s}  {'series where it no-ops':>24s}\")\nfor fov in (160, 150, 140, 130, 120):\n    need = fov / ok[\"spacing\"]\n    clipped = (need > ok[\"rows_px\"]).mean()\n    print(f\"{fov:4d}mm  {clipped:23.1%}\")\nprint(\"\\n-> 160mm no-ops on a non-trivial fraction. I went with 130mm.\")\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 6. Which plane shows which finding\n\nA knee MRI is acquired from several directions: sagittal (S), coronal (C),\naxial (A).\n\nRather than assume \"feed it everything\", I built versions with planes removed and\n**submitted them**, to measure which findings depend on which plane.\n\n| planes used | LB |\n|---|---|\n| coronal only (C) | 0.718 |\n| sagittal only (S) | 0.772 |\n| sagittal + axial (S+A) | 0.802 |\n| **all three** | **0.820** |\n\nBroken down per finding, the ordering lands where anatomy says it should. These\nare **local validation AUCs per finding**, a different quantity from the LB\nnumbers above.\n\n| finding | best single plane (validation AUC) | anatomical expectation |\n|---|---|---|\n| MCL | **coronal 0.705** (sagittal is worst at 0.477) | the MCL is read on coronal |\n| Medial Meniscus | **sagittal 0.729** | menisci are read on sagittal |\n| PF OA (patellofemoral) | **axial 0.629** (sagittal 0.590 / coronal 0.528) | cartilage behind the patella is axial |\n| Synovitis | **axial 0.804** | joint fluid and synovium are axial |\n\nThat is evidence the model is looking somewhere anatomically sensible, rather\nthan latching onto a shortcut such as a site-specific imaging quirk.\n\nIf you only add one plane to sagittal, add axial.\n\n## 7. How to split validation\n\nImages come from many sites and scanners. If the same scanner appears on both\nsides of a split, the model can learn \"this is scanner X, and scanner X exams\ntend to look like this\", and validation flatters it.\n\nSo I build a **scanner fingerprint** from the DICOM headers\n(`Manufacturer | ModelName | SoftwareVersions | MagneticFieldStrength`), which\ngives 112 groups, and split with `GroupKFold` so **no group straddles the\nsplit**. Validation is then always \"a scanner the model has not seen\".\n"},{"cell_type":"code","metadata":{},"source":"# ---- how to split validation ----------------------------------------------\n# Group by scanner fingerprint and split so no group straddles the boundary.\n# The real fingerprint comes from one pass over every series header; here we\n# just demonstrate the splitting discipline.\nfrom sklearn.model_selection import GroupKFold\n\ndemo = train[[\"StudyInstanceUID\"]].copy()\ndemo[\"group\"] = np.random.default_rng(0).integers(0, 112, len(demo))  # 112 real groups\n\ngkf = GroupKFold(n_splits=5)\nfor f, (tr_i, va_i) in enumerate(gkf.split(demo, groups=demo[\"group\"])):\n    a = set(demo.iloc[tr_i][\"group\"])\n    b = set(demo.iloc[va_i][\"group\"])\n    assert not (a & b), \"a group straddles train and validation\"\n    print(f\"  fold {f}: train {len(tr_i):,} / valid {len(va_i):,}  overlap {len(a & b)}\")\nprint(\"\\n-> Validation is always an unseen scanner.\")\n","execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{},"source":"## 8. All of that fixes the baseline\n\nThe EDA findings turn directly into settings.\n\n| decision | because |\n|---|---|\n| Training target = labels built from the reports | only 58 exams are labelled (1) |\n| Crop a fixed 130mm physical field of view | 160mm no-ops on 28.9% (5) |\n| Use all three planes | different planes carry different findings (6) |\n| GroupKFold on scanner fingerprint, 5 folds | stop the model learning the scanner (7) |\n| Average folds; blend on ranks | AUC only sees ranking (0) |\n\nWhat's left is the image model, and there I use a conventional setup.\n\n### Building the input (2.5D)\n\nMRI is 3D, but 3D CNNs are heavy and there isn't much training data. So: **2.5D**.\n\n1. Pick one representative series per plane\n2. Sample **24 slices** at even spacing\n3. Take **3 consecutive slices as one group**, treated as the RGB channels (8 groups)\n4. 3 planes x 8 groups = 24 \"images\" make up one exam\n\nFor a single plane:\n\n```\noriginal series (tens of slices)\n  ...............................\n\n        | pick 24 at even spacing\n        v\n  .  .  .  .  .  .  .  .   ... (24)\n\n        | consecutive 3 -> RGB\n        v\n  [1,2,3] [4,5,6] [7,8,9] ...    8 groups\n\n        | encoder, per group\n        v\n  one feature vector per group\n        |\n        | learn which groups matter, per finding\n        v\n  probability for each of the 12 findings\n```\n\nPutting neighbouring slices in the channels lets a plain 2D model see\n**continuity perpendicular to the plane**. A structure like the meniscus, which\nonly appears on certain slices, ends up concentrated in the groups containing\nthose slices.\n\n### Model\n\n- **Encoder**: DINOv2 ViT-S/14 (self-supervised pretraining)\n- **Pooling**: learn which of the 8 groups to weight **separately for each\n  finding** (the best slices for the meniscus are not the best slices for the ACL)\n- **Learning rate**: 1e-3 on the head, 8e-6 on the encoder, and only the **last\n  6 blocks** unfrozen\n\nThat last point matters. Fine-tuning the whole pretrained encoder at a normal\nlearning rate **overwrites the pretrained weights and destroys the\nrepresentation**. The version that trained all layers at 1e-4 lost to a far\nsmaller resnet18. Adapt the encoder; don't retrain it.\n\n## 9. What each change was worth\n\nNone of this was decided up front. I submitted one change at a time and measured\neach contribution.\n\n| change | public LB |\n|---|---|\n| resnet18 2.5D + weak labels from my own rules | 0.820 |\n| **target replaced with LLM-read labels from the reports** | **0.857** (+0.037) |\n| crop 160mm -> 130mm | 0.861 (+0.004) |\n| encoder -> DINOv2 ViT-S/14 with discriminative LRs | 0.882 (+0.021) |\n| pooling made per-finding | 0.884 (+0.002) |\n| input resolution 224px -> **336px** | **0.898** (+0.014) |\n\n**The biggest single move was the label swap (+0.037)** - larger than every image\nmodelling improvement put together. The claim in section 1, that this is a\nlabel-building problem, shows up directly in the results.\n\nResolution was second (+0.014). 130mm at 336px works out to **0.387mm per\npixel**, and I read that as roughly what it takes to see a meniscal tear.\n\nBecause the changes were separated, each contribution is a number. Bundled\ntogether, there would be no way to tell whether the gain came from the labels or\nthe resolution - and no basis for deciding where to spend effort next.\n\n## 10. Where to go from here\n\n- **Improve label quality.** It was the largest axis by a wide margin, and\n  better report reading still has room.\n- **Train several models and blend.** Though: blending in a model that is weak on\n  its own does not help. What helps is a model that is **about as strong as your\n  baseline and wrong in different places**.\n- **Blend with public models.** Someone else's model trained with a different\n  design differs from yours far more than diversity you manufacture in-house.\n"},{"cell_type":"markdown","metadata":{},"source":"## Appendix: running this baseline end to end\n\nEverything above is analysis. The cells below actually run the baseline described\nin sections 8-9 and write `submission.csv`, so the notebook produces the **public\nLB 0.898** entry rather than just describing it.\n\nThe weights are the 5 folds of that configuration, published as a CC0 dataset:\n`kitopl/rsna-knee-eda-baseline-weights`. They were trained on the public\nreport labels from `stevenleehans/rsna-knee-llm-report-labels`, so nothing here\ndepends on private data.\n\nTwo things worth repeating because they are easy to get wrong:\n\n- **Preprocessing at inference must match training exactly.** The 130mm crop, the\n  24-slice sampling, the 1st/99th percentile normalisation - if any of it drifts,\n  the score drops and nothing tells you why. The code below is a verbatim copy of\n  the training-time preprocessing.\n- **The encoder is built with `pretrained=False`.** The backbone weights are\n  already inside the checkpoint, so no internet access is needed.\n"},{"cell_type":"code","metadata":{},"source":"# ---- inference config ------------------------------------------------------\nimport os\nimport time\nfrom concurrent.futures import ProcessPoolExecutor\n\nimport cv2\nimport pydicom\nimport timm\nimport torch\nimport torch.nn as nn\nfrom scipy.stats import rankdata\n\nOUT = \"/kaggle/working\"\nINPUT = str(ROOT)\n\nPLANES_ORDER = [\"Sagittal\", \"Coronal\", \"Axial\"]\nPLANES, N_SLICES, HW = 3, 24, 336\nTRIPLETS = N_SLICES // 3\nFOV_MM = 130.0          # must match training (section 5)\nBACKBONE = \"vit_small_patch14_dinov2.lvd142m\"\n\n\ndef find_weights(fname=\"model_fold0.pt\"):\n    \"\"\"Locate the weights by searching /kaggle/input.\n\n    Hardcoding the mount path is brittle - it depends on how the dataset was\n    attached. Search for the file instead, and print what is actually there so\n    a miss is diagnosable rather than mysterious.\n    \"\"\"\n    for root, _, files in os.walk(\"/kaggle/input\"):\n        if fname in files:\n            return root\n    return None\n\nBATCH = 8\nTIME_BUDGET = 7.5 * 3600   # margin against the 9h limit; fill the rest with 0.5\n\ndev = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nprint(f\"device={dev}  torch={torch.__version__}  timm={timm.__version__}\")\nprint(\"mounted under /kaggle/input:\", sorted(os.listdir(\"/kaggle/input\")))\nWEIGHTS = find_weights()\nprint(\"weights dir:\", WEIGHTS)\nassert WEIGHTS, \"model_fold0.pt not found anywhere under /kaggle/input\"\nprint(\"weights files:\", sorted(f for f in os.listdir(WEIGHTS) if f.endswith(\".pt\")))\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- preprocessing: verbatim from training ---------------------------------\n# Any divergence between this and the training-time preprocessing costs score\n# silently, so it is copied rather than re-derived.\n\ndef pick_series(series_df):\n    \"\"\"One series per (study, plane), preferring fluid-sensitive.\"\"\"\n    df = series_df.copy()\n    df[\"_pref\"] = -df[\"Fluid_Sensitive\"].fillna(0).astype(int)\n    return (df.sort_values([\"StudyInstanceUID\", \"Anatomical_Plane\",\n                            \"_pref\", \"SeriesInstanceUID\"])\n              .groupby([\"StudyInstanceUID\", \"Anatomical_Plane\"], as_index=False)\n              .head(1))\n\n\ndef slice_position(ds):\n    \"\"\"Projection of the image position onto the slice normal.\n\n    Do not use the z component directly: for sagittal the slice axis is x and\n    for coronal it is y.\n    \"\"\"\n    iop = getattr(ds, \"ImageOrientationPatient\", None)\n    ipp = getattr(ds, \"ImagePositionPatient\", None)\n    if iop is None or ipp is None:\n        return None\n    n = np.cross([float(v) for v in iop[:3]], [float(v) for v in iop[3:]])\n    return float(np.dot([float(v) for v in ipp], n))\n\n\ndef load_volume(split, study, series):\n    d = os.path.join(INPUT, f\"{split}_series\", study, series)\n    try:\n        paths = [f.path for f in os.scandir(d) if f.name.endswith(\".dcm\")]\n    except OSError as e:\n        return None, f\"scandir: {e}\"\n    if not paths:\n        return None, \"no dcm\"\n\n    items = []\n    for p in paths:\n        try:\n            ds = pydicom.dcmread(p, force=True)\n            arr = ds.pixel_array\n        except Exception:\n            continue\n        ps = getattr(ds, \"PixelSpacing\", None)\n        items.append((slice_position(ds), arr,\n                      float(ps[0]) if ps is not None else None))\n    if not items:\n        return None, \"all slices failed to decode\"\n\n    if all(it[0] is None for it in items):\n        items = [(i, a, s) for i, (_, a, s) in enumerate(items)]\n    else:\n        items = [it for it in items if it[0] is not None]\n        items.sort(key=lambda x: x[0])\n\n    out = np.zeros((N_SLICES, HW, HW), dtype=np.uint8)\n    idx = np.linspace(0, len(items) - 1, N_SLICES).round().astype(int)\n    for k, i in enumerate(idx):\n        _, a, spacing = items[i]\n        a = a.astype(np.float32)\n        if a.ndim != 2:\n            a = a[..., 0] if a.ndim == 3 else a.reshape(a.shape[-2:])\n        if spacing and spacing > 0:\n            # crop a fixed physical field of view, then resize (section 5)\n            crop_px = int(round(FOV_MM / spacing))\n            crop_px = max(16, min(crop_px, min(a.shape)))\n            cy, cx = a.shape[0] // 2, a.shape[1] // 2\n            h = crop_px // 2\n            a = a[max(0, cy - h):cy + h, max(0, cx - h):cx + h]\n        flat = a[::4, ::4].ravel()\n        lo, hi = np.percentile(flat, [1, 99]) if flat.size else (0.0, 1.0)\n        a = np.clip((a - lo) / max(hi - lo, 1e-6), 0, 1)\n        out[k] = (cv2.resize(a, (HW, HW),\n                             interpolation=cv2.INTER_AREA) * 255).astype(np.uint8)\n    return out, None\n\n\ndef build_study(args):\n    study, chosen = args\n    vol = np.zeros((PLANES, N_SLICES, HW, HW), dtype=np.uint8)\n    for pi, plane in enumerate(PLANES_ORDER):\n        sid = chosen.get(plane)\n        if sid is None:\n            continue\n        v, _ = load_volume(\"test\", study, sid)\n        if v is not None:\n            vol[pi] = v\n    return study, vol\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- model: must match training token for token -----------------------------\n# If this drifts from the training definition, the state_dict keys or shapes\n# stop lining up and inference breaks quietly.\n\nPRETRAINED = False       # backbone weights live in the checkpoint\nUSE_TRIPLET_POS = False\n\n\ndef build_encoder(backbone):\n    kw = dict(pretrained=PRETRAINED, num_classes=0, in_chans=3)\n    if backbone.startswith(\"vit_\"):\n        # DINOv2 is patch14 and pretrained at 518px; timm resamples the\n        # positional embedding when img_size is given.\n        assert HW % 14 == 0, f\"patch14 ViT needs HW divisible by 14: {HW}\"\n        kw[\"img_size\"] = HW\n    return timm.create_model(backbone, **kw)\n\n\nclass Net(nn.Module):\n    \"\"\"2.5D encoder with per-finding attention over triplets.\n\n    The attention produces one weight per (plane, triplet, finding), so the\n    slices that matter for the meniscus can differ from the ones that matter\n    for the ACL. Collapsing to a single vector before the classifier, as an\n    earlier version did, throws that away.\n    \"\"\"\n\n    def __init__(self, backbone, n_out=len(TARGETS)):\n        super().__init__()\n        self.enc = build_encoder(backbone)\n        d = self.enc.num_features\n        self.n_out = n_out\n        self.att = nn.Sequential(nn.Linear(d, 128), nn.Tanh(), nn.Linear(128, n_out))\n        self.pos = nn.Parameter(torch.zeros(1, 1, TRIPLETS, d))\n        self.drop = nn.Dropout(0.2)\n        self.cls_w = nn.Parameter(torch.empty(n_out, PLANES * d))\n        self.cls_b = nn.Parameter(torch.zeros(n_out))\n        nn.init.normal_(self.cls_w, std=0.02)\n\n    def forward(self, x):                      # (B, 3, 24, H, W)\n        B = x.shape[0]\n        x = x.reshape(B * PLANES * TRIPLETS, 3, HW, HW)\n        f = self.enc(x).reshape(B, PLANES, TRIPLETS, -1)\n        if USE_TRIPLET_POS:\n            f = f + self.pos\n        w = torch.softmax(self.att(f), dim=2)          # normalise over triplets\n        p = torch.einsum(\"bptd,bptc->bpcd\", f, w)      # per-finding plane embedding\n        p = self.drop(p.permute(0, 2, 1, 3).reshape(B, self.n_out, -1))\n        return (p * self.cls_w.unsqueeze(0)).sum(-1) + self.cls_b\n\n\nnets = []\nfor f in range(5):\n    p = os.path.join(WEIGHTS, f\"model_fold{f}.pt\")\n    if not os.path.exists(p):\n        continue\n    m = Net(BACKBONE).to(dev).eval()\n    m.load_state_dict(torch.load(p, map_location=dev))\n    nets.append(m)\nassert nets, f\"no weights loaded from {WEIGHTS}\"\nprint(f\"loaded {len(nets)} folds of {BACKBONE}\")\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- run inference and write submission.csv --------------------------------\nt0 = time.time()\n\ntest = pd.read_csv(f\"{INPUT}/test.csv\")\nte_series = pd.read_csv(f\"{INPUT}/test_series.csv\")\nprint(f\"test: {len(test):,} studies / {len(te_series):,} series\")\n\nsel = pick_series(te_series)\ntasks = [(s, dict(zip(g.Anatomical_Plane, g.SeriesInstanceUID)))\n         for s, g in sel.groupby(\"StudyInstanceUID\")]\nprint(f\"selected {len(sel):,} series for {len(tasks):,} studies\")\n\norder, chunks = [], []\nworkers = max(1, (os.cpu_count() or 4))\nwith ProcessPoolExecutor(max_workers=workers) as ex:\n    buf_u, buf_x = [], []\n    for i, (study, vol) in enumerate(ex.map(build_study, tasks, chunksize=4), 1):\n        buf_u.append(study)\n        buf_x.append(vol)\n        last = (i == len(tasks))\n        if len(buf_u) < BATCH and not last:\n            continue\n\n        x = torch.from_numpy(np.stack(buf_x)).float().div_(255.0)\n        x = ((x - 0.449) / 0.226).to(dev)\n        with torch.no_grad(), torch.amp.autocast(\"cuda\", enabled=(dev == \"cuda\")):\n            # average the 5 folds in probability space, as at training time\n            p = sum(torch.sigmoid(m(x)) for m in nets) / len(nets)\n        chunks.append(p.float().cpu().numpy())\n        order += buf_u\n        buf_u, buf_x = [], []\n\n        if i % 200 == 0 or last:\n            el = time.time() - t0\n            print(f\"    {i:,}/{len(tasks):,}  {el/60:.1f} min  \"\n                  f\"eta {(len(tasks)-i)*el/max(i,1)/60:.0f} min\", flush=True)\n        if time.time() - t0 > TIME_BUDGET:\n            print(\"[budget] out of time; remaining studies get 0.5\")\n            break\n\nP = np.concatenate(chunks) if chunks else np.zeros((0, len(TARGETS)))\nn = len(order)\nassert len(P) == n, f\"{len(P)} != {n}\"\n\n# AUC only sees ranking, so convert to per-column ranks before writing.\n# With a single model this is a monotone transform and does not change the\n# score; it keeps the code identical to the multi-model blend.\npred = np.column_stack([rankdata(P[:, j]) / max(n, 1) for j in range(len(TARGETS))])\npreds = dict(zip(order, pred))\n\nsub = pd.DataFrame({\"StudyInstanceUID\": test.StudyInstanceUID})\narr = np.full((len(sub), len(TARGETS)), 0.5, np.float64)\nfor i, u in enumerate(sub.StudyInstanceUID):\n    if u in preds:\n        arr[i] = preds[u]\nfor j, t in enumerate(TARGETS):\n    sub[t] = arr[:, j]\n\nsample = pd.read_csv(f\"{INPUT}/sample_submission.csv\")\nsub = sub[list(sample.columns)]\nassert list(sub.columns) == list(sample.columns)\nassert len(sub) == len(test)\nassert sub[TARGETS].notna().all().all()\nsub.to_csv(f\"{OUT}/submission.csv\", index=False)\n\nn_pred = sum(1 for u in sub.StudyInstanceUID if u in preds)\nprint(f\"\\nwrote submission.csv {sub.shape}  predicted {n_pred:,}/{len(sub):,} \"\n      f\"(filled with 0.5: {len(sub)-n_pred:,})\")\nprint(f\"done in {(time.time()-t0)/60:.1f} min\")\n","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11"}},"nbformat":4,"nbformat_minor":5}