{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","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"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"3232d69d-071c-4c77-8051-565ed794bea8","cell_type":"markdown","source":"# RSNA Knee MRI — Preprocessing Pipeline (v1, built from the EDA results)\n\nThis notebook turns 4 407 studies / 24 371 series of raw DICOM into a compact, geometry-normalised\n`uint8` cache plus study-level metadata, ready for training. **No model is trained here.**\n\n## What the EDA told us, and what each finding forced\n\n| Finding (measured) | Consequence for this pipeline |\n|---|---|\n| **Only 58 / 4 407 studies are labelled (1.3 %)** | The images alone cannot be trained supervised. We cache **all** studies and also emit a cleaned report table, because the report→label module is now the critical path. |\n| Rule-based label recovery is weak (Medial/Lateral OA recall **0.00**, best F1 ≈ 0.6) | Regex weak-labelling is not enough — we keep the raw + cleaned report so a multilingual encoder can be trained on it later. |\n| Reports: 36 % English, rest Spanish / Turkish / German / Dutch / French / Greek / Cyrillic | Report cleaning is language-agnostic (placeholder stripping, whitespace, unicode normalisation only). No English-specific tokenisation. |\n| **Plane label agrees with the direction cosines 100 %**, only 0.9 % oblique | We can trust `Anatomical_Plane` from the CSV (available for train *and* test) and skip plane inference. |\n| **In-plane orientation is 100 % consistent** — Sag (row→+y, col→−z), Cor (row→+x, col→−z), Ax (row→+x, col→+y) | No per-study in-plane flip logic is needed. The only mirroring required is **laterality**. |\n| **`InstanceNumber` order matches geometric order only 63.8 %** of the time | Slices are sorted by projection onto the slice normal. Never by filename or InstanceNumber. |\n| 0 % duplicate slice positions, 0 % mid-series matrix changes, step std ≈ 0 | Volume assembly is simple; a de-dup guard is kept but will rarely fire. |\n| **`Laterality` present on only 50 % of slices**; L/R separate cleanly on the isocentre x-coordinate (L ≈ +31 mm, R ≈ −130 mm) | Three-source laterality resolution: DICOM tag → report regex → geometry, with the geometry threshold **fitted and validated on the tagged half**. |\n| FOV 130–220 mm (1st–99th pct) but matrix 192–1280 → **mm/pixel varies 0.137–1.172** | A plain `resize()` would put the same anatomy at up to 2× different scales. We **resample to a fixed mm/pixel and crop a fixed physical FOV**. |\n| In-plane spacing is isotropic in 100 % of slices | Single scale factor per slice; no anisotropic warping. |\n| Slice step median 3.5 mm, stack span median 101 mm, 30 slices | Through-plane: fixed physical extent, fixed slice count, **nearest-position sampling** (no interpolation — it blurs menisci). |\n| PhotometricInterpretation MONOCHROME2 100 %; RescaleIntercept always 0; 25 % signed pixels, vmin down to −622 | Apply slope, clamp negatives, keep a MONOCHROME1 guard for the unseen test set. |\n| p99 intensity ranges 198 → 5 717 across series | **Per-series** percentile normalisation is mandatory. No global constants. |\n| Only 6 (plane × contrast) cells exist; `Axial\\|F0\\|S0` present in just 19.5 % of studies | We cache **5 canonical cells** and drop Axial/non-fat-sat. Availability mask is stored per study. |\n| All 3 planes present in 100 % of studies | Plane-level fallbacks will effectively never fire, but are implemented anyway. |\n| `SeriesDescription` is `DUMMYSERIESDESC!` or blank for ~21 % | Series selection must **not** depend on the description. We use geometry + slice count. |\n| `PatientID` fill rate 1.0 | Collected for **grouped CV** — some patients may contribute several studies. |\n| Decode 2.2 ms/slice, all sampled series Explicit VR LE | Full-train preprocessing ≈ 30 min single-thread; trivially parallel. |\n| ⚠️ **`pylibjpeg` / `gdcm` are NOT installed** in this image | Every sampled series was uncompressed, but the data description promises JPEG Lossless and JPEG 2000 somewhere. See §0.1 — this must be fixed before submitting. |\n\n## What comes out\n\n```\n/kaggle/working/prep_out/\n├── cache/<StudyInstanceUID>.npz     # 5 cells, uint8, canonical geometry + availability mask\n├── study_meta.parquet               # laterality (+source), scanner, patient, cell availability, QC flags\n├── series_index.parquet             # which series was chosen for each cell, and why\n├── reports_clean.parquet            # cleaned report text + language guess + report-derived laterality\n└── preprocess_config.json           # every constant used, so the test-time path is reproducible\n```","metadata":{}},{"id":"3f0b598d-9795-45d9-83c4-7075c5eecf8b","cell_type":"markdown","source":"## 0. Configuration\n\nEverything the pipeline does is determined by the constants below, and they are all written to\n`preprocess_config.json` so the inference kernel can reproduce the exact same transform.\n\n**On the geometry targets.** `FOV_MM = 160` is the median measured field of view, and the measured content\nbounding box was 159 mm tall × 138 mm wide on average — so a 160 mm window captures the whole knee with\nmargin, while the 99th-percentile 220 mm studies simply get their empty periphery cropped away. `OUT_SIZE`\nthen fixes the scale: 192 px over 160 mm = **0.833 mm/pixel for every study in the dataset**, which is the\nwhole point.\n\n**On the slice extents.** Measured stack span is 101 mm (median) at 3.5 mm steps. The extents below stay\ninside that so we are sub-sampling, not padding, for the typical study.","metadata":{}},{"id":"28e3f255-e029-4014-a72c-b61395ac3e4e","cell_type":"code","source":"import os, sys, re, gc, json, time, math, glob, unicodedata, warnings, traceback\nfrom pathlib import Path\nfrom collections import Counter, defaultdict\nwarnings.filterwarnings(\"ignore\")\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport pydicom, cv2\n\ndef find_root():\n    for c in [\"/kaggle/input/competitions/rsna-knee-abnormality-detection\",\n              \"/kaggle/input/rsna-knee-abnormality-detection\"]:\n        if os.path.exists(os.path.join(c, \"train.csv\")): return c\n    for p in sorted(glob.glob(\"/kaggle/input/**/train.csv\", recursive=True)):\n        return os.path.dirname(p)\n    raise FileNotFoundError(\"train.csv not found\")\n\nDATA = find_root()\nOUT  = Path(\"/kaggle/working/prep_out\"); (OUT / \"cache\").mkdir(parents=True, exist_ok=True)\n\nLABELS = [\"ACL\",\"MCL\",\"Medial Meniscus\",\"Lateral Meniscus\",\"Medial OA\",\"Lateral OA\",\n          \"PF OA\",\"Effusion\",\"Synovitis\",\"Baker's\",\"Contusion\",\"Fracture\"]\n\nCFG = {\n    # ---- spatial normalisation ------------------------------------------------\n    \"OUT_SIZE\": 192,          # output matrix (px). 160mm / 192px = 0.833 mm/px for EVERY study\n    \"FOV_MM\": 160.0,          # physical field of view kept, in mm (measured median FOV = 160)\n    # ---- the 5 canonical (plane, fluid, fatsat) cells -------------------------\n    # name        plane       fluid fatsat  n_slices  extent_mm   fallback cell\n    \"CELLS\": [\n        (\"SAG_FS\", \"Sagittal\", 1, 1, 24, 92.0, \"SAG_T1\"),\n        (\"SAG_T1\", \"Sagittal\", 0, 0, 24, 92.0, \"SAG_FS\"),\n        (\"COR_FS\", \"Coronal\",  1, 1, 20, 76.0, \"COR_T1\"),\n        (\"COR_T1\", \"Coronal\",  0, 0, 20, 76.0, \"COR_FS\"),\n        (\"AX_FS\",  \"Axial\",    1, 1, 16, 96.0, None),\n    ],\n    # ---- intensity ------------------------------------------------------------\n    \"CLIP_LO_PCT\": 0.5,       # per-SERIES percentiles over foreground voxels\n    \"CLIP_HI_PCT\": 99.5,\n    \"FG_THRESH_FRAC\": 0.06,   # foreground = value > lo + frac*(hi-lo), used for centroid/percentiles\n    # ---- laterality -----------------------------------------------------------\n    \"CANONICAL_SIDE\": \"R\",    # mirror everything into right-knee convention\n    \"LAT_GEOM_THRESH\": None,  # fitted in section 3 from the tagged half; None until then\n    # ---- run control ----------------------------------------------------------\n    \"N_WORKERS\": 4,\n    \"PILOT_N\": 40,            # studies processed in the pilot before the full run\n    \"COMPRESS\": True,         # np.savez_compressed vs np.savez\n    \"SEED\": 42,\n}\nnp.random.seed(CFG[\"SEED\"])\n\nCELL_NAMES = [c[0] for c in CFG[\"CELLS\"]]\nCELL_BY_NAME = {c[0]: dict(zip([\"name\",\"plane\",\"fluid\",\"fatsat\",\"n_slices\",\"extent_mm\",\"fallback\"], c))\n                for c in CFG[\"CELLS\"]}\n\nprint(\"DATA:\", DATA)\nprint(\"mm per pixel  :\", round(CFG[\"FOV_MM\"] / CFG[\"OUT_SIZE\"], 4))\nprint(\"cells         :\", CELL_NAMES)\nprint(\"slices/study  :\", sum(c[4] for c in CFG[\"CELLS\"]))\nprint(\"bytes/study   :\", f\"{sum(c[4] for c in CFG['CELLS']) * CFG['OUT_SIZE']**2 / 1e6:.2f} MB (uncompressed)\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:43.263961Z","iopub.execute_input":"2026-08-06T22:16:43.264696Z","iopub.status.idle":"2026-08-06T22:16:44.922567Z","shell.execute_reply.started":"2026-08-06T22:16:43.264653Z","shell.execute_reply":"2026-08-06T22:16:44.921788Z"}},"outputs":[],"execution_count":null},{"id":"63a77e9e-3178-4992-8ac4-970c1584aeeb","cell_type":"markdown","source":"### 0.1 ⚠️ Decoder check — read this before you submit anything\n\nThe EDA found `pylibjpeg`, `pylibjpeg-libjpeg`, `pylibjpeg-openjpeg` and `gdcm` all **absent**, while every one\nof the 2 193 sampled series was uncompressed Explicit VR Little Endian. The competition description explicitly\nlists JPEG Lossless and JPEG 2000 among the transfer syntaxes present. So either the compressed series are a\nrare tail we did not sample, or they live in the hidden test set — and we cannot tell which.\n\nEither way the fix is the same and it is cheap: attach the codec wheels as an offline Kaggle dataset and\n`pip install --no-index` them in the inference kernel. The cell below fails loudly rather than silently\nproducing an all-`0.5` submission six hours into scoring.","metadata":{}},{"id":"c715abff-49eb-4df5-92cb-9aa55c2d8156","cell_type":"code","source":"def decoder_status():\n    st = {}\n    for m in [\"pylibjpeg\", \"libjpeg\", \"openjpeg\", \"gdcm\"]:\n        try:\n            __import__(m); st[m] = True\n        except Exception:\n            st[m] = False\n    return st\n\nDEC = decoder_status()\nprint(DEC)\nif not (DEC[\"gdcm\"] or (DEC[\"libjpeg\"] and DEC[\"openjpeg\"])):\n    print(\"\\n\" + \"!\" * 78)\n    print(\"NO JPEG/JPEG2000 DECODER AVAILABLE.\")\n    print(\"Uncompressed DICOMs will still work, so this notebook will run to completion,\")\n    print(\"but any compressed series in the hidden test set will raise on .pixel_array.\")\n    print(\"Before submitting: attach a dataset containing the wheels for\")\n    print(\"  pylibjpeg, pylibjpeg-libjpeg, pylibjpeg-openjpeg   (or python-gdcm)\")\n    print(\"and pip install --no-index --find-links=<that dataset> in the inference kernel.\")\n    print(\"!\" * 78)\n\n# the pipeline records every decode failure instead of crashing, see decode stats in section 6","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:44.924238Z","iopub.execute_input":"2026-08-06T22:16:44.924677Z","iopub.status.idle":"2026-08-06T22:16:44.933663Z","shell.execute_reply.started":"2026-08-06T22:16:44.924647Z","shell.execute_reply":"2026-08-06T22:16:44.932734Z"}},"outputs":[],"execution_count":null},{"id":"760475cc-7868-4441-a9fd-b7ec98422918","cell_type":"markdown","source":"## 1. Load the tables and pick one series per canonical cell\n\nThere are 24 371 series for 4 407 studies and 5 cells — so several studies have **more than one series per\ncell** (`Sagittal|F0|S0` alone appears 5 197 times for 4 407 studies). We must choose deterministically.\n\nThe selection score deliberately **ignores `SeriesDescription`**, because 21 % of them are `DUMMYSERIESDESC!`\nor blank. Instead we rank on slice count: a knee series has ~30 slices (measured median), so we penalise the\ndistance from 30 in log space. That naturally demotes 3-slice localisers and 320-slice dynamic/3D acquisitions\nwithout hard-coding thresholds. Ties break on the lexicographically smallest `SeriesInstanceUID`, so the\nchoice is reproducible across runs and across machines.\n\nWe count slices from the directory listing (one `.dcm` = one slice — confirmed by the EDA: `NumberOfFrames`\nis 1 or absent everywhere), which is far cheaper than opening files.","metadata":{}},{"id":"f975d6ed-963a-41f3-beb0-0007ca4dbf60","cell_type":"code","source":"train  = pd.read_csv(os.path.join(DATA, \"train.csv\"))\ntser   = pd.read_csv(os.path.join(DATA, \"train_series.csv\"))\ntest   = pd.read_csv(os.path.join(DATA, \"test.csv\"))\ntests  = pd.read_csv(os.path.join(DATA, \"test_series.csv\"))\nprint(train.shape, tser.shape, test.shape, tests.shape)\n\ndef count_slices(root, study, series):\n    d = os.path.join(root, study, series)\n    try:\n        return sum(1 for f in os.scandir(d) if f.name.endswith(\".dcm\"))\n    except FileNotFoundError:\n        return 0\n\nfrom concurrent.futures import ThreadPoolExecutor\n\ndef add_slice_counts(df, root):\n    with ThreadPoolExecutor(max_workers=16) as ex:\n        n = list(ex.map(lambda r: count_slices(root, r[0], r[1]),\n                        df[[\"StudyInstanceUID\", \"SeriesInstanceUID\"]].values))\n    out = df.copy(); out[\"n_slices\"] = n\n    return out\n\nt0 = time.time()\ntser = add_slice_counts(tser, os.path.join(DATA, \"train_series\"))\nprint(f\"slice counts in {time.time()-t0:.0f}s; empty series: {(tser['n_slices']==0).sum()}\")\nprint(tser[\"n_slices\"].describe().round(1))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:44.934655Z","iopub.execute_input":"2026-08-06T22:16:44.934956Z","iopub.status.idle":"2026-08-06T22:16:58.045777Z","shell.execute_reply.started":"2026-08-06T22:16:44.934931Z","shell.execute_reply":"2026-08-06T22:16:58.044968Z"}},"outputs":[],"execution_count":null},{"id":"a72b6837-7a92-49ad-8a04-114126eb3994","cell_type":"code","source":"def cell_key(row):\n    for name, plane, fl, fs, *_ in CFG[\"CELLS\"]:\n        if row[\"Anatomical_Plane\"] == plane and int(row[\"Fluid_Sensitive\"]) == fl and int(row[\"Fat_Suppression\"]) == fs:\n            return name\n    return None\n\ndef build_series_index(series_df):\n    df = series_df.copy()\n    df[\"cell\"] = df.apply(cell_key, axis=1)\n    df = df[df[\"cell\"].notna() & (df[\"n_slices\"] > 0)].copy()\n    # prefer series whose slice count is closest to the measured median of 30, in log space\n    df[\"score\"] = -np.abs(np.log(df[\"n_slices\"].clip(lower=1) / 30.0))\n    df = df.sort_values([\"StudyInstanceUID\", \"cell\", \"score\", \"SeriesInstanceUID\"],\n                        ascending=[True, True, False, True])\n    chosen = df.groupby([\"StudyInstanceUID\", \"cell\"], as_index=False).first()\n    return df, chosen\n\nall_cells, SEL = build_series_index(tser)\nprint(\"chosen rows:\", len(SEL))\npiv = SEL.pivot_table(index=\"StudyInstanceUID\", columns=\"cell\", values=\"n_slices\", aggfunc=\"first\")\npiv = piv.reindex(columns=CELL_NAMES)\nprint(\"\\ncell availability after selection:\")\nprint((piv.notna().mean() * 100).round(1).astype(str) + \" %\")\nprint(\"\\ncells per study:\")\nprint(piv.notna().sum(axis=1).value_counts().sort_index())\n\n# how many series were discarded as duplicates within a cell?\ndups = all_cells.groupby([\"StudyInstanceUID\",\"cell\"]).size()\nprint(\"\\nstudy-cells with >1 candidate series:\", int((dups > 1).sum()),\n      f\"({(dups > 1).mean()*100:.1f}% of study-cells)\")\nprint(\"discarded series:\", len(all_cells) - len(SEL))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:58.047866Z","iopub.execute_input":"2026-08-06T22:16:58.048129Z","iopub.status.idle":"2026-08-06T22:16:58.544221Z","shell.execute_reply.started":"2026-08-06T22:16:58.048106Z","shell.execute_reply":"2026-08-06T22:16:58.543474Z"}},"outputs":[],"execution_count":null},{"id":"1526a473-3e86-4b08-96a2-cb6b527e2ea1","cell_type":"markdown","source":"## 2. Reports — cleaning and language tagging\n\nWe do **not** attempt to derive labels here. The EDA showed regex rules top out around 0.6 precision and score\n0.00 recall on both OA compartment labels, so labelling is a separate modelling problem. What this section\nproduces is the clean input for that model:\n\n* de-identification placeholders (`[DATE]`, `[TIME]`, `[ID]`, `[NAME]`, …) removed — they are pure noise and\n  appear in a large fraction of reports;\n* NFKC unicode normalisation, so that Greek and Cyrillic text tokenises consistently;\n* whitespace collapsed, but **line structure preserved** as `\\n`, because many reports are sectioned\n  (`FINDINGS:` / `Hallazgos:` / `ΕΥΡΗΜΑΤΑ`) and that structure is signal;\n* a language guess, so a multilingual model can be evaluated per language later.\n\nCrucially it also extracts **laterality from the report text**, which is our second-best source when the\nDICOM `Laterality` tag is missing (it is missing half the time).","metadata":{}},{"id":"6d1e74d1-e652-4662-a464-487316c137c2","cell_type":"code","source":"PLACEHOLDER = re.compile(r\"\\[(DATE|TIME|ID|NAME|AGE|PHONE|ADDRESS|DATETIME|MRN|ACCESSION)\\]\", re.I)\nWS = re.compile(r\"[ \\t\\u00a0]+\")\nNL = re.compile(r\"\\n{3,}\")\n\ndef clean_report(t):\n    if not isinstance(t, str): return \"\"\n    t = unicodedata.normalize(\"NFKC\", t)\n    t = PLACEHOLDER.sub(\" \", t)\n    t = t.replace(\"\\r\", \"\\n\")\n    t = WS.sub(\" \", t)\n    t = NL.sub(\"\\n\\n\", t)\n    return t.strip()\n\ndef script_of(t):\n    if re.search(r\"[\\u0370-\\u03ff\\u1f00-\\u1fff]\", t): return \"greek\"\n    if re.search(r\"[\\u0400-\\u04ff]\", t): return \"cyrillic\"\n    if re.search(r\"[\\u0600-\\u06ff]\", t): return \"arabic\"\n    return \"latin\"\n\nSTOP = {\n \"en\": {\"the\",\"and\",\"with\",\"there\",\"normal\",\"tear\",\"joint\",\"signal\",\"findings\",\"knee\"},\n \"es\": {\"el\",\"la\",\"los\",\"con\",\"del\",\"sin\",\"rodilla\",\"articular\",\"menisco\",\"hallazgos\"},\n \"fr\": {\"le\",\"les\",\"des\",\"avec\",\"sans\",\"genou\",\"articulaire\",\"menisque\",\"aspect\"},\n \"de\": {\"der\",\"die\",\"das\",\"und\",\"mit\",\"ohne\",\"kniegelenk\",\"kein\",\"nachweis\",\"befund\"},\n \"it\": {\"il\",\"del\",\"della\",\"con\",\"senza\",\"ginocchio\",\"menisco\",\"articolare\",\"non\"},\n \"pt\": {\"do\",\"da\",\"com\",\"sem\",\"joelho\",\"menisco\",\"articular\",\"nao\",\"em\"},\n \"nl\": {\"de\",\"het\",\"een\",\"met\",\"zonder\",\"knie\",\"geen\",\"van\",\"bij\",\"links\"},\n \"tr\": {\"ve\",\"ile\",\"yok\",\"diz\",\"eklem\",\"izlendi\",\"olan\",\"mevcut\",\"sinyal\"},\n}\ndef guess_lang(t):\n    s = script_of(t)\n    if s != \"latin\": return s\n    toks = set(re.findall(r\"[a-zA-Zçğşöüıáéíóúñäöüß]+\", t.lower()))\n    best, n = \"unk\", 0\n    for k, v in STOP.items():\n        m = len(toks & v)\n        if m > n: best, n = k, m\n    return best if n >= 2 else \"unk\"\n\n# ---- laterality from report text, multilingual ---------------------------\nLAT_PAT = {\n \"R\": r\"\\b(right|rt\\.?|derech[ao]|direita|droit[e]?|rechts?|destr[ao]|sağ|sag\\b|правый|правая|справа|δεξι)\",\n \"L\": r\"\\b(left|lt\\.?|izquierd[ao]|esquerda|gauche|links?|sinistr[ao]|sol\\b|левый|левая|слева|αριστερ)\",\n}\ndef lat_from_report(t):\n    tl = t.lower()\n    r = len(re.findall(LAT_PAT[\"R\"], tl)); l = len(re.findall(LAT_PAT[\"L\"], tl))\n    if r > l: return \"R\"\n    if l > r: return \"L\"\n    return None\n\nREP = pd.DataFrame({\"StudyInstanceUID\": train[\"StudyInstanceUID\"]})\nREP[\"report_raw\"]   = train[\"Report\"].fillna(\"\")\nREP[\"report_clean\"] = REP[\"report_raw\"].map(clean_report)\nREP[\"lang\"]         = REP[\"report_clean\"].map(guess_lang)\nREP[\"n_chars\"]      = REP[\"report_clean\"].str.len()\nREP[\"lat_report\"]   = REP[\"report_clean\"].map(lat_from_report)\nif \"PatientSex\" in train.columns:\n    REP[\"PatientSex\"] = train[\"PatientSex\"]\n\nprint(REP[\"lang\"].value_counts())\nprint(\"\\nreport laterality found:\", REP[\"lat_report\"].notna().mean().round(3))\nprint(REP[\"lat_report\"].value_counts(dropna=False))\ndisplay(REP[[\"report_clean\",\"lang\",\"lat_report\"]].head(3))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:58.545177Z","iopub.execute_input":"2026-08-06T22:16:58.545704Z","iopub.status.idle":"2026-08-06T22:16:59.69872Z","shell.execute_reply.started":"2026-08-06T22:16:58.54568Z","shell.execute_reply":"2026-08-06T22:16:59.698014Z"}},"outputs":[],"execution_count":null},{"id":"779e49bb-4b5d-4f28-8ac7-69fa3efedb78","cell_type":"markdown","source":"## 3. Laterality resolution — the single most important transform\n\nFour of the twelve labels are compartment-specific (Medial/Lateral Meniscus, Medial/Lateral OA). In a coronal\nor axial image the medial compartment sits on the image *left* for one knee and the image *right* for the\nother; in a sagittal stack, the medial end is at slice 0 for one knee and slice N for the other. Without\nmirroring, the network has to learn a mirrored copy of every medial/lateral feature — an unaffordable waste\nagainst **58 labelled studies**.\n\nWhy mirroring is *sufficient* here, and why it is just an x-flip: the EDA found the direction cosines are\n100 % consistent, with exactly one sign combination per plane —\n\n* **Sagittal**: column axis → +y (posterior), row axis → −z (inferior) ⇒ normal = **−x** (patient right)\n* **Coronal**: column axis → +x (patient left), row axis → −z ⇒ normal = **+y** (posterior)\n* **Axial**: column axis → +x, row axis → +y ⇒ normal = **+z** (superior)\n\nSo after sorting slices by their projection on the normal, patient-space direction is fully determined, and\ncanonicalising to a right knee is a single mirror of the patient x-axis:\n\n* Coronal / Axial → flip the **column** axis (`np.flip(..., axis=2)`)\n* Sagittal → flip the **slice** axis (`np.flip(..., axis=0)`) — the in-plane axes are y,z and untouched\n\nResolution order per study: **DICOM `Laterality` tag** (majority vote over the study) → **report text** →\n**geometry**. For geometry we compute the image-centre x in patient coordinates and fit the decision threshold\non the ~50 % of studies that *do* carry the tag, then report held-out accuracy. If the geometry rule is not\nconvincingly accurate we would rather mark the study `unknown` and leave it unmirrored (with a flag the model\ncan consume) than mirror it wrongly.","metadata":{}},{"id":"3ad9d001-3a68-4673-a736-6ef9604adb73","cell_type":"code","source":"def study_lat_probe(root, study, series_ids):\n    # cheap: read one header from each of up to 3 series\n    tags, xs = [], []\n    for se in series_ids[:3]:\n        d = os.path.join(root, study, se)\n        try:\n            files = sorted(f.name for f in os.scandir(d) if f.name.endswith(\".dcm\"))\n        except FileNotFoundError:\n            continue\n        if not files: continue\n        for f in [files[0], files[len(files)//2]]:\n            try:\n                ds = pydicom.dcmread(os.path.join(d, f), stop_before_pixels=True, force=True)\n            except Exception:\n                continue\n            lt = str(getattr(ds, \"Laterality\", \"\") or \"\").strip().upper()\n            if lt in (\"L\", \"R\"): tags.append(lt)\n            iop = getattr(ds, \"ImageOrientationPatient\", None)\n            ipp = getattr(ds, \"ImagePositionPatient\", None)\n            ps  = getattr(ds, \"PixelSpacing\", None)\n            cols = getattr(ds, \"Columns\", None); rows = getattr(ds, \"Rows\", None)\n            if iop is None or ipp is None or ps is None: continue\n            iop = np.array(iop, float); ipp = np.array(ipp, float); ps = np.array(ps, float)\n            # centre of the image in patient space\n            centre = ipp + iop[:3] * (float(cols) / 2 * ps[1]) + iop[3:] * (float(rows) / 2 * ps[0])\n            xs.append(float(centre[0]))\n    tag = Counter(tags).most_common(1)[0][0] if tags else None\n    return tag, (float(np.mean(xs)) if xs else np.nan)\n\nstudies = SEL[\"StudyInstanceUID\"].unique()\nser_by_study = SEL.groupby(\"StudyInstanceUID\")[\"SeriesInstanceUID\"].apply(list).to_dict()\nROOT_TR = os.path.join(DATA, \"train_series\")\n\nt0 = time.time(); probe = {}\nwith ThreadPoolExecutor(max_workers=12) as ex:\n    futs = {ex.submit(study_lat_probe, ROOT_TR, s, ser_by_study[s]): s for s in studies}\n    for i, f in enumerate(futs):\n        pass\n    for f, s in futs.items():\n        probe[s] = f.result()\nprint(f\"probed {len(probe)} studies in {time.time()-t0:.0f}s\")\n\nLAT = pd.DataFrame([{\"StudyInstanceUID\": s, \"lat_tag\": v[0], \"centre_x\": v[1]} for s, v in probe.items()])\nLAT = LAT.merge(REP[[\"StudyInstanceUID\", \"lat_report\"]], on=\"StudyInstanceUID\", how=\"left\")\nprint(\"\\ntag available :\", LAT[\"lat_tag\"].notna().mean().round(3))\nprint(\"report available:\", LAT[\"lat_report\"].notna().mean().round(3))\nprint(\"centre_x available:\", LAT[\"centre_x\"].notna().mean().round(3))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:16:59.699726Z","iopub.execute_input":"2026-08-06T22:16:59.700068Z","iopub.status.idle":"2026-08-06T22:18:11.013463Z","shell.execute_reply.started":"2026-08-06T22:16:59.700032Z","shell.execute_reply":"2026-08-06T22:18:11.012551Z"}},"outputs":[],"execution_count":null},{"id":"62541373-1d85-41fb-9349-51f2e523ea75","cell_type":"code","source":"# ---- how good is the report against the tag? ----------------------------\nboth = LAT.dropna(subset=[\"lat_tag\", \"lat_report\"])\nprint(f\"report vs tag agreement: {(both['lat_tag'] == both['lat_report']).mean()*100:.1f}%  (n={len(both)})\")\nprint(pd.crosstab(both[\"lat_tag\"], both[\"lat_report\"]))\n\n# ---- fit the geometry threshold on the tagged studies -------------------\ng = LAT.dropna(subset=[\"lat_tag\", \"centre_x\"])\ncand = np.arange(-200, 100, 1.0)\nacc = [( (np.where(g[\"centre_x\"].values > t, \"L\", \"R\") == g[\"lat_tag\"].values).mean(), t) for t in cand]\nbest_acc, best_t = max(acc)\nCFG[\"LAT_GEOM_THRESH\"] = float(best_t)\nprint(f\"\\ngeometry rule: centre_x > {best_t:.1f} mm  =>  LEFT\")\nprint(f\"accuracy on tagged studies: {best_acc*100:.2f}%  (n={len(g)})\")\n\nplt.figure(figsize=(8, 3.4))\nfor k, grp in g.groupby(\"lat_tag\"):\n    plt.hist(grp[\"centre_x\"], bins=60, alpha=0.6, label=f\"tag={k}\")\nplt.axvline(best_t, color=\"k\", ls=\"--\", label=f\"threshold {best_t:.0f}\")\nplt.legend(); plt.xlabel(\"image centre x in patient space (mm)\")\nplt.title(\"Geometry-based laterality separation\"); plt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:18:11.015124Z","iopub.execute_input":"2026-08-06T22:18:11.015992Z","iopub.status.idle":"2026-08-06T22:18:11.533294Z","shell.execute_reply.started":"2026-08-06T22:18:11.015956Z","shell.execute_reply":"2026-08-06T22:18:11.532471Z"}},"outputs":[],"execution_count":null},{"id":"dde416d5-e586-46b4-a5e9-9326432ede2d","cell_type":"code","source":"GEOM_MIN_ACC = 0.95   # below this we refuse to trust geometry and mark the study unknown\n\ndef resolve_lat(row):\n    if isinstance(row[\"lat_tag\"], str):    return row[\"lat_tag\"], \"tag\"\n    if isinstance(row[\"lat_report\"], str): return row[\"lat_report\"], \"report\"\n    if best_acc >= GEOM_MIN_ACC and row[\"centre_x\"] == row[\"centre_x\"]:\n        return (\"L\" if row[\"centre_x\"] > CFG[\"LAT_GEOM_THRESH\"] else \"R\"), \"geometry\"\n    return \"U\", \"unknown\"\n\nres = LAT.apply(resolve_lat, axis=1, result_type=\"expand\")\nLAT[\"laterality\"], LAT[\"lat_source\"] = res[0], res[1]\nprint(LAT[\"lat_source\"].value_counts())\nprint()\nprint(LAT[\"laterality\"].value_counts())\nprint(f\"\\nstudies that will be mirrored (side != {CFG['CANONICAL_SIDE']}): \"\n      f\"{(LAT['laterality'] == ('L' if CFG['CANONICAL_SIDE']=='R' else 'R')).sum()}\")\nprint(f\"studies left unmirrored due to unknown side: {(LAT['laterality']=='U').sum()}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:18:11.534491Z","iopub.execute_input":"2026-08-06T22:18:11.534856Z","iopub.status.idle":"2026-08-06T22:18:11.692835Z","shell.execute_reply.started":"2026-08-06T22:18:11.534829Z","shell.execute_reply":"2026-08-06T22:18:11.692121Z"}},"outputs":[],"execution_count":null},{"id":"3122bc70-8fd1-46d7-846b-16a332ffdb4f","cell_type":"markdown","source":"## 4. The core transform\n\nOne function, used identically for train and test. Steps, in order:\n\n1. **Read headers only** for every slice in the series (`stop_before_pixels=True`, ~0.1 ms each) to get\n   `ImagePositionPatient`, `ImageOrientationPatient`, `PixelSpacing`, matrix, rescale and photometric fields.\n2. **Sort by projection onto the slice normal.** Not by filename, not by `InstanceNumber` — that would be wrong\n   for 36 % of series. De-duplicate identical positions (never observed, kept as a guard).\n3. **Select slice indices** at a fixed physical step so that `n_slices` slices span `extent_mm`, centred on the\n   stack centre, using **nearest position** rather than interpolation — interpolating along a 3.5 mm axis\n   smears the meniscus, which is the thing we are trying to see. If the stack is shorter than the requested\n   extent we sample the whole stack uniformly and record the actual step in the metadata.\n4. **Decode only the selected slices** (16–24 of ~30, so we skip ~30 % of the decode cost), apply\n   `RescaleSlope`, invert if `MONOCHROME1`, clamp negatives to 0.\n5. **Resample in-plane** to `FOV_MM / OUT_SIZE` mm per pixel with `INTER_AREA` when downscaling (correct\n   anti-aliasing; measured spacings mean we downscale in the large majority of cases) and `INTER_LINEAR` when\n   upscaling.\n6. **Crop / pad** to `OUT_SIZE`, centred on the foreground centroid of the middle slice — the same window for\n   every slice in the series, so the volume stays geometrically coherent. Measured content centroids sat at\n   (0.50, 0.51) of the frame, so this is a small correction, but it protects the off-centre tail.\n7. **Canonicalise laterality** by mirroring the patient x-axis (column flip for Coronal/Axial, slice flip for\n   Sagittal).\n8. **Normalise intensity per series**: percentiles 0.5/99.5 over foreground voxels only, clip, scale to\n   `uint8`. Per-series, because the measured p99 spans 198 → 5 717 across series and there is no\n   Hounsfield-like absolute scale in MR.","metadata":{}},{"id":"b8d9246f-8103-43a2-9b17-f97b80298370","cell_type":"code","source":"def _read_geom(path):\n    ds = pydicom.dcmread(path, stop_before_pixels=True, force=True)\n    iop = getattr(ds, \"ImageOrientationPatient\", None)\n    ipp = getattr(ds, \"ImagePositionPatient\", None)\n    ps  = getattr(ds, \"PixelSpacing\", None)\n    if iop is None or ipp is None or ps is None: return None\n    return dict(path=path,\n                iop=np.array(iop, float), ipp=np.array(ipp, float), ps=np.array(ps, float),\n                rows=int(ds.Rows), cols=int(ds.Columns),\n                slope=float(getattr(ds, \"RescaleSlope\", 1) or 1),\n                inter=float(getattr(ds, \"RescaleIntercept\", 0) or 0),\n                photo=str(getattr(ds, \"PhotometricInterpretation\", \"MONOCHROME2\")),\n                inst=float(getattr(ds, \"InstanceNumber\", 0) or 0))\n\ndef _pixels(g):\n    ds = pydicom.dcmread(g[\"path\"], force=True)\n    a = ds.pixel_array.astype(np.float32)\n    if g[\"slope\"] != 1 or g[\"inter\"] != 0:\n        a = a * g[\"slope\"] + g[\"inter\"]\n    if g[\"photo\"].upper() == \"MONOCHROME1\":\n        a = a.max() - a\n    np.clip(a, 0, None, out=a)\n    return a\n\ndef _resample_crop(a, ps_rc, out_size, mm_per_px, centre_rc=None):\n    # a: (H,W) float32; ps_rc: (row_spacing, col_spacing) in mm\n    sh = max(1, int(round(a.shape[0] * ps_rc[0] / mm_per_px)))\n    sw = max(1, int(round(a.shape[1] * ps_rc[1] / mm_per_px)))\n    interp = cv2.INTER_AREA if (sh < a.shape[0] or sw < a.shape[1]) else cv2.INTER_LINEAR\n    r = cv2.resize(a, (sw, sh), interpolation=interp)\n    if centre_rc is None:\n        centre_rc = (sh / 2.0, sw / 2.0)\n    r0 = int(round(centre_rc[0] - out_size / 2)); c0 = int(round(centre_rc[1] - out_size / 2))\n    out = np.zeros((out_size, out_size), np.float32)\n    sr0, sc0 = max(0, r0), max(0, c0)\n    sr1, sc1 = min(sh, r0 + out_size), min(sw, c0 + out_size)\n    if sr1 > sr0 and sc1 > sc0:\n        out[sr0 - r0:sr1 - r0, sc0 - c0:sc1 - c0] = r[sr0:sr1, sc0:sc1]\n    return out\n\ndef _fg_centroid(a, frac):\n    lo, hi = float(a.min()), float(np.percentile(a, 99.5))\n    m = a > lo + frac * (hi - lo)\n    if m.sum() < 64: return None\n    rs, cs = np.nonzero(m)\n    return (float(rs.mean()), float(cs.mean()))\n\ndef load_cell(root, study, series, cell, laterality, cfg=CFG):\n    spec = CELL_BY_NAME[cell]\n    d = os.path.join(root, study, series)\n    files = sorted(f.path for f in os.scandir(d) if f.name.endswith(\".dcm\"))\n    if not files: raise RuntimeError(\"no dcm files\")\n\n    geoms = [g for g in (_read_geom(p) for p in files) if g is not None]\n    if len(geoms) < 2: raise RuntimeError(\"insufficient geometry\")\n\n    iop0 = geoms[0][\"iop\"]\n    n = np.cross(iop0[:3], iop0[3:]); n /= (np.linalg.norm(n) + 1e-9)\n    t = np.array([float(np.dot(g[\"ipp\"], n)) for g in geoms])\n    order = np.argsort(t); geoms = [geoms[i] for i in order]; t = t[order]\n\n    keep = [0] + [i for i in range(1, len(t)) if abs(t[i] - t[i-1]) > 1e-3]   # de-dup guard\n    geoms = [geoms[i] for i in keep]; t = t[keep]\n    n_dup = len(order) - len(keep)\n\n    N, ext = spec[\"n_slices\"], spec[\"extent_mm\"]\n    span = float(t[-1] - t[0])\n    if span >= ext and len(t) > 1:\n        centre = 0.5 * (t[0] + t[-1])\n        want = centre + np.linspace(-ext / 2, ext / 2, N)\n        step_mm = ext / max(N - 1, 1)\n    else:\n        want = np.linspace(t[0], t[-1], N)\n        step_mm = span / max(N - 1, 1)\n    idx = [int(np.argmin(np.abs(t - w))) for w in want]\n\n    mm = cfg[\"FOV_MM\"] / cfg[\"OUT_SIZE\"]\n    mid = geoms[idx[len(idx) // 2]]\n    cen = _fg_centroid(_pixels(mid), cfg[\"FG_THRESH_FRAC\"])\n    if cen is not None:\n        cen = (cen[0] * mid[\"ps\"][0] / mm, cen[1] * mid[\"ps\"][1] / mm)\n\n    vol = np.empty((N, cfg[\"OUT_SIZE\"], cfg[\"OUT_SIZE\"]), np.float32)\n    for k, i in enumerate(idx):\n        g = geoms[i]\n        vol[k] = _resample_crop(_pixels(g), g[\"ps\"], cfg[\"OUT_SIZE\"], mm, cen)\n\n    # laterality canonicalisation: mirror the patient x-axis\n    mirrored = False\n    if laterality in (\"L\", \"R\") and laterality != cfg[\"CANONICAL_SIDE\"]:\n        vol = np.flip(vol, axis=0) if spec[\"plane\"] == \"Sagittal\" else np.flip(vol, axis=2)\n        vol = np.ascontiguousarray(vol); mirrored = True\n\n    # per-series intensity normalisation over foreground voxels\n    lo0, hi0 = float(vol.min()), float(np.percentile(vol, 99.5))\n    fg = vol[vol > lo0 + cfg[\"FG_THRESH_FRAC\"] * (hi0 - lo0)]\n    if fg.size < 1000: fg = vol.reshape(-1)\n    lo = float(np.percentile(fg, cfg[\"CLIP_LO_PCT\"])); hi = float(np.percentile(fg, cfg[\"CLIP_HI_PCT\"]))\n    if hi - lo < 1e-6: hi = lo + 1.0\n    vol = np.clip((vol - lo) / (hi - lo), 0, 1)\n    vol = (vol * 255.0 + 0.5).astype(np.uint8)\n\n    meta = dict(cell=cell, series=series, n_src_slices=len(t), span_mm=round(span, 2),\n                step_mm=round(float(step_mm), 3),\n                px_mm_src=round(float(geoms[0][\"ps\"][0]), 4),\n                src_rows=geoms[0][\"rows\"], src_cols=geoms[0][\"cols\"],\n                mirrored=mirrored, n_dup=n_dup,\n                clip_lo=round(lo, 2), clip_hi=round(hi, 2),\n                short_stack=bool(span < ext))\n    return vol, meta","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:18:11.693971Z","iopub.execute_input":"2026-08-06T22:18:11.69423Z","iopub.status.idle":"2026-08-06T22:18:11.72225Z","shell.execute_reply.started":"2026-08-06T22:18:11.694205Z","shell.execute_reply":"2026-08-06T22:18:11.721348Z"}},"outputs":[],"execution_count":null},{"id":"3b9f7896-c1b5-409a-a550-25649edd1dbf","cell_type":"code","source":"def process_study(args):\n    study, rows, laterality, root, cfg = args\n    out, meta, errs = {}, {}, {}\n    have = {}\n    for cell in CELL_NAMES:\n        se = rows.get(cell)\n        if se is None: continue\n        try:\n            v, m = load_cell(root, study, se, cell, laterality, cfg)\n            out[cell] = v; meta[cell] = m; have[cell] = 1\n        except Exception as e:\n            errs[cell] = f\"{type(e).__name__}: {e}\"[:160]\n    # fallback: fill a missing cell from its same-plane partner, flagged in the mask\n    for cell in CELL_NAMES:\n        if cell in out: continue\n        fb = CELL_BY_NAME[cell][\"fallback\"]\n        spec = CELL_BY_NAME[cell]\n        if fb and fb in out and out[fb].shape[0] == spec[\"n_slices\"]:\n            out[cell] = out[fb].copy(); have[cell] = 2      # 2 = filled from fallback\n            meta[cell] = dict(meta[fb], cell=cell, filled_from=fb)\n        else:\n            out[cell] = np.zeros((spec[\"n_slices\"], cfg[\"OUT_SIZE\"], cfg[\"OUT_SIZE\"]), np.uint8)\n            have[cell] = 0                                   # 0 = absent, zero-filled\n    return study, out, meta, have, errs\n\ndef save_study(study, out, have, cfg=CFG):\n    p = OUT / \"cache\" / f\"{study}.npz\"\n    payload = {c: out[c] for c in CELL_NAMES}\n    payload[\"avail\"] = np.array([have[c] for c in CELL_NAMES], np.uint8)\n    (np.savez_compressed if cfg[\"COMPRESS\"] else np.savez)(p, **payload)\n    return p.stat().st_size\n\ndef make_jobs(sel_df, lat_df, root, cfg=CFG):\n    lat = dict(zip(lat_df[\"StudyInstanceUID\"], lat_df[\"laterality\"]))\n    grp = sel_df.groupby(\"StudyInstanceUID\").apply(\n        lambda d: dict(zip(d[\"cell\"], d[\"SeriesInstanceUID\"]))).to_dict()\n    return [(s, rows, lat.get(s, \"U\"), root, cfg) for s, rows in grp.items()]\n\nJOBS = make_jobs(SEL, LAT, ROOT_TR)\nprint(\"jobs:\", len(JOBS))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:18:11.725185Z","iopub.execute_input":"2026-08-06T22:18:11.725557Z","iopub.status.idle":"2026-08-06T22:18:12.044253Z","shell.execute_reply.started":"2026-08-06T22:18:11.725523Z","shell.execute_reply":"2026-08-06T22:18:12.043412Z"}},"outputs":[],"execution_count":null},{"id":"37fb405d-a404-40ed-b8fc-4a53e7cc0cf5","cell_type":"markdown","source":"## 5. Pilot run — verify the transform before spending an hour on it\n\nWe process a handful of studies, time them, measure the on-disk size, and then **look at the result**. The\nqualitative checks that matter:\n\n1. Every cell of a study should be the same anatomy at the same scale.\n2. Mirroring should work: the mean canonicalised coronal image of originally-left knees and of\n   originally-right knees should look like *the same knee*. If mirroring is inverted they will be mirror\n   images of each other — this is the check that catches a sign error, and it is much more reliable than\n   eyeballing individual slices.","metadata":{}},{"id":"9aa8a93e-e84d-4aed-886a-e751896e256d","cell_type":"code","source":"pilot = JOBS[:CFG[\"PILOT_N\"]]\nt0 = time.time(); sizes = []; pilot_meta = []; pilot_err = []\nresults = {}\nfor j in pilot:\n    st, out, meta, have, errs = process_study(j)\n    sizes.append(save_study(st, out, have))\n    pilot_meta.append(dict(StudyInstanceUID=st, **{f\"avail_{c}\": have[c] for c in CELL_NAMES}))\n    for c, e in errs.items(): pilot_err.append((st, c, e))\n    if len(results) < 3: results[st] = (out, have, meta)\ndt = time.time() - t0\nprint(f\"{len(pilot)} studies in {dt:.1f}s  ->  {dt/len(pilot)*1000:.0f} ms/study (1 worker)\")\nprint(f\"mean npz size: {np.mean(sizes)/1e6:.2f} MB\")\nprint(f\"projected cache for {len(JOBS)} studies: {np.mean(sizes)*len(JOBS)/1e9:.1f} GB\")\nprint(f\"projected wall time at {CFG['N_WORKERS']} workers: \"\n      f\"{dt/len(pilot)*len(JOBS)/CFG['N_WORKERS']/60:.1f} min\")\nif pilot_err:\n    print(\"\\nerrors:\"); [print(\" \", e) for e in pilot_err[:10]]\nelse:\n    print(\"\\nno errors in pilot\")\nprint(pd.DataFrame(pilot_meta)[[f\"avail_{c}\" for c in CELL_NAMES]].apply(pd.Series.value_counts).fillna(0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-06T22:18:12.045362Z","iopub.execute_input":"2026-08-06T22:18:12.045684Z"}},"outputs":[],"execution_count":null},{"id":"dc26c2bf-2a2d-4f85-acd3-ab76ae9a2774","cell_type":"code","source":"# visual check 1: all five cells of one study, mid-slice\nst, (out, have, meta) = list(results.items())[0]\nfig, ax = plt.subplots(2, len(CELL_NAMES), figsize=(3.0*len(CELL_NAMES), 6.4))\nfor i, c in enumerate(CELL_NAMES):\n    v = out[c]; z = v.shape[0] // 2\n    ax[0, i].imshow(v[z], cmap=\"gray\", vmin=0, vmax=255); ax[0, i].axis(\"off\")\n    ax[0, i].set_title(f\"{c}\\navail={have[c]} z={z}\", fontsize=8)\n    ax[1, i].imshow(v[max(0, z - 5)], cmap=\"gray\", vmin=0, vmax=255); ax[1, i].axis(\"off\")\n    ax[1, i].set_title(f\"z={max(0,z-5)}\", fontsize=8)\nlat_here = LAT.set_index(\"StudyInstanceUID\").loc[st]\nfig.suptitle(f\"{st[:26]}...  laterality={lat_here['laterality']} ({lat_here['lat_source']})\", fontsize=10)\nplt.tight_layout(); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"9ffc3b70-bf71-46eb-a600-ac8b1ef49300","cell_type":"code","source":"# visual check 2: through-stack view of one sagittal cell\nv = out[\"SAG_FS\"]\nk = v.shape[0]; ncol = 8; nrow = int(np.ceil(k / ncol))\nfig, axes = plt.subplots(nrow, ncol, figsize=(1.9*ncol, 2.0*nrow))\nfor a, i in zip(np.atleast_1d(axes).ravel(), range(k)):\n    a.imshow(v[i], cmap=\"gray\", vmin=0, vmax=255); a.set_title(str(i), fontsize=7); a.axis(\"off\")\nfor a in np.atleast_1d(axes).ravel()[k:]: a.axis(\"off\")\nfig.suptitle(\"SAG_FS canonicalised: index 0 should be MEDIAL for every study\", fontsize=10)\nplt.tight_layout(); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"16925c24-e4da-4578-99ef-aed96e6ae710","cell_type":"code","source":"# visual check 3: the mirroring test -- mean coronal image, L-origin vs R-origin studies\nlat_map = dict(zip(LAT[\"StudyInstanceUID\"], LAT[\"laterality\"]))\nacc = {\"L\": [], \"R\": []}\nfor j in JOBS[:80]:\n    s = j[0]; side = lat_map.get(s, \"U\")\n    if side not in acc or len(acc[side]) >= 12: continue\n    try:\n        v, _ = load_cell(ROOT_TR, s, j[1][\"COR_FS\"], \"COR_FS\", side)\n        acc[side].append(v[v.shape[0] // 2].astype(np.float32))\n    except Exception:\n        pass\n    if len(acc[\"L\"]) >= 12 and len(acc[\"R\"]) >= 12: break\n\nfig, ax = plt.subplots(1, 3, figsize=(12, 4.2))\nmL = np.mean(acc[\"L\"], axis=0) if acc[\"L\"] else np.zeros((CFG[\"OUT_SIZE\"],)*2)\nmR = np.mean(acc[\"R\"], axis=0) if acc[\"R\"] else np.zeros((CFG[\"OUT_SIZE\"],)*2)\nax[0].imshow(mL, cmap=\"gray\"); ax[0].set_title(f\"mean COR_FS, originally LEFT (n={len(acc['L'])})\"); ax[0].axis(\"off\")\nax[1].imshow(mR, cmap=\"gray\"); ax[1].set_title(f\"mean COR_FS, originally RIGHT (n={len(acc['R'])})\"); ax[1].axis(\"off\")\nax[2].imshow(np.abs(mL - mR), cmap=\"magma\"); ax[2].set_title(\"absolute difference\"); ax[2].axis(\"off\")\nplt.tight_layout(); plt.show()\n\nd_same = float(np.abs(mL - mR).mean())\nd_flip = float(np.abs(mL - mR[:, ::-1]).mean())\nprint(f\"mean |L - R|            = {d_same:.2f}\")\nprint(f\"mean |L - flip(R)|      = {d_flip:.2f}\")\nprint(\"=> mirroring is CORRECT\" if d_same < d_flip else\n      \"=> WARNING: flipping R matches better; the mirroring convention is inverted\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"cd472c38-2c24-48f9-ab75-c3eb735c27cc","cell_type":"markdown","source":"## 6. Full run\n\nParallel over studies, resumable (existing `.npz` files are skipped), and it records every failure rather than\ndying on one bad series. Expect roughly the pilot's per-study time divided by the worker count.\n\nIf the projected cache size printed in the pilot exceeds ~18 GB, drop `OUT_SIZE` to 160 or reduce the\nper-cell slice counts before running this — `/kaggle/working` is capped at 20 GB.","metadata":{}},{"id":"9280df7a-8249-4b32-a32f-556311563ca5","cell_type":"code","source":"from concurrent.futures import ProcessPoolExecutor\n\ntodo = [j for j in JOBS if not (OUT / \"cache\" / f\"{j[0]}.npz\").exists()]\nprint(f\"{len(JOBS) - len(todo)} already cached, {len(todo)} to do\")\n\nrows_meta, rows_series, failures = [], [], []\nt0 = time.time()\nwith ProcessPoolExecutor(max_workers=CFG[\"N_WORKERS\"]) as ex:\n    for i, (st, out, meta, have, errs) in enumerate(ex.map(process_study, todo, chunksize=4)):\n        try:\n            sz = save_study(st, out, have)\n        except Exception as e:\n            failures.append((st, \"save\", str(e)[:150])); continue\n        rows_meta.append(dict(StudyInstanceUID=st, npz_bytes=sz,\n                              **{f\"avail_{c}\": have[c] for c in CELL_NAMES}))\n        for c, m in meta.items():\n            rows_series.append(dict(StudyInstanceUID=st, **m))\n        for c, e in errs.items():\n            failures.append((st, c, e))\n        if (i + 1) % 250 == 0:\n            el = time.time() - t0\n            print(f\"  {i+1}/{len(todo)}  {el/60:.1f} min  eta {el/(i+1)*(len(todo)-i-1)/60:.1f} min\")\n\nprint(f\"\\ndone in {(time.time()-t0)/60:.1f} min | failures: {len(failures)}\")\nif failures:\n    F = pd.DataFrame(failures, columns=[\"StudyInstanceUID\",\"cell\",\"error\"])\n    print(F[\"error\"].str.split(\":\").str[0].value_counts())\n    F.to_csv(OUT / \"failures.csv\", index=False)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"404b3a05-fc70-4c37-8de6-758d6935a4f7","cell_type":"code","source":"META = pd.DataFrame(rows_meta)\nSERM = pd.DataFrame(rows_series)\n\nSTUDY_META = (LAT[[\"StudyInstanceUID\",\"laterality\",\"lat_source\",\"lat_tag\",\"lat_report\",\"centre_x\"]]\n              .merge(META, on=\"StudyInstanceUID\", how=\"right\"))\nSTUDY_META = STUDY_META.merge(REP[[\"StudyInstanceUID\",\"lang\",\"n_chars\"]], on=\"StudyInstanceUID\", how=\"left\")\nif \"PatientSex\" in train.columns:\n    STUDY_META = STUDY_META.merge(train[[\"StudyInstanceUID\",\"PatientSex\"]], on=\"StudyInstanceUID\", how=\"left\")\n\n# gold labels, where they exist\nlab_mask = train[LABELS].notna().all(axis=1)\nSTUDY_META = STUDY_META.merge(train.loc[lab_mask, [\"StudyInstanceUID\"] + LABELS],\n                              on=\"StudyInstanceUID\", how=\"left\")\nSTUDY_META[\"has_gold\"] = STUDY_META[LABELS].notna().all(axis=1).astype(int)\n\nprint(\"cached studies:\", len(STUDY_META), \"| with gold labels:\", int(STUDY_META[\"has_gold\"].sum()))\nprint(\"\\ncell availability (0=absent, 1=real, 2=fallback-filled):\")\nprint(STUDY_META[[f\"avail_{c}\" for c in CELL_NAMES]].apply(pd.Series.value_counts).fillna(0).astype(int))\nprint(f\"\\ntotal cache: {STUDY_META['npz_bytes'].sum()/1e9:.2f} GB \"\n      f\"({STUDY_META['npz_bytes'].mean()/1e6:.2f} MB/study)\")\n\nSTUDY_META.to_parquet(OUT / \"study_meta.parquet\", index=False)\nSERM.to_parquet(OUT / \"series_index.parquet\", index=False)\nREP.to_parquet(OUT / \"reports_clean.parquet\", index=False)\nwith open(OUT / \"preprocess_config.json\", \"w\") as f:\n    json.dump({**CFG, \"geom_thresh_accuracy\": float(best_acc), \"data_root\": DATA,\n               \"cells\": CELL_NAMES, \"generated\": time.strftime(\"%Y-%m-%d %H:%M\")}, f, indent=2, default=str)\nprint(\"\\nwritten:\", [p.name for p in sorted(OUT.iterdir())])","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"f958555d-5096-4e2d-8641-81d386d94e56","cell_type":"markdown","source":"## 7. Post-hoc QC on the cache\n\nThree things worth checking after the fact, because they are cheap and each one has bitten real pipelines:\n\n* **Intensity distribution across studies** — after per-series normalisation the histograms should look\n  similar regardless of scanner. If Siemens and GE still separate, normalisation is not doing its job.\n* **Short stacks** — series whose physical span was smaller than the requested extent got the whole stack\n  resampled instead. Too many of those means the extent constants are too aggressive.\n* **Blank / near-blank volumes** — a cell that came out almost all zeros indicates a decode or crop failure\n  that did not raise.","metadata":{}},{"id":"6936be65-c054-47ee-ab6c-c1f60295b773","cell_type":"code","source":"samp = STUDY_META.sample(min(200, len(STUDY_META)), random_state=0)[\"StudyInstanceUID\"]\nmeans, blanks = [], []\nfor s in samp:\n    z = np.load(OUT / \"cache\" / f\"{s}.npz\")\n    for c in CELL_NAMES:\n        v = z[c]\n        means.append(dict(StudyInstanceUID=s, cell=c, mean=float(v.mean()),\n                          p99=float(np.percentile(v, 99)), frac_zero=float((v == 0).mean())))\n        if v.mean() < 3: blanks.append((s, c))\nQ = pd.DataFrame(means)\nprint(\"near-blank cells:\", len(blanks))\nprint(Q.groupby(\"cell\")[[\"mean\",\"p99\",\"frac_zero\"]].describe().round(2).T.head(20))\n\nfig, ax = plt.subplots(1, 3, figsize=(16, 3.8))\nfor c in CELL_NAMES:\n    ax[0].hist(Q.loc[Q.cell == c, \"mean\"], bins=30, alpha=0.5, label=c)\nax[0].legend(fontsize=7); ax[0].set_title(\"mean intensity per cell (uint8)\")\nax[1].hist(Q[\"frac_zero\"], bins=40, color=\"#a85c48\"); ax[1].set_title(\"fraction of exactly-zero voxels\")\nif len(SERM):\n    ax[2].hist(SERM[\"step_mm\"].clip(0, 10), bins=40, color=\"#487a5c\")\n    ax[2].set_title(\"effective slice step after resampling (mm)\")\nplt.tight_layout(); plt.show()\n\nif len(SERM):\n    print(\"\\nshort stacks (span < requested extent):\",\n          f\"{SERM['short_stack'].mean()*100:.1f}% of cells\")\n    print(SERM.groupby(\"cell\")[\"short_stack\"].mean().round(3))\n    print(\"\\nmirrored cells:\", f\"{SERM['mirrored'].mean()*100:.1f}%\")\n    print(\"source in-plane spacing after selection:\")\n    print(SERM[\"px_mm_src\"].describe().round(3))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"c7264539-3558-4b48-bfad-9ada6de4bf8f","cell_type":"markdown","source":"## 8. Loader + the test-time path\n\n`load_study` is what the training `Dataset` will call. The tensor layout is\n`(n_cells, n_slices, OUT_SIZE, OUT_SIZE)` — ragged in the slice dimension across cells, so it is returned as a\ndict rather than one array; a 2.5-D model will consume each cell as its own multi-channel input, and a 3-D\nmodel will take them as separate volumes.\n\n`preprocess_test_study` runs the identical transform straight from DICOM with no cache, which is what the\ninference kernel needs. Note it resolves laterality from tag → geometry only (no report is available at test\ntime), which is why §3 validated the geometry rule independently.","metadata":{}},{"id":"53557ffa-8668-4b03-b117-02b52fa4e64e","cell_type":"code","source":"def load_study(study, cache_dir=OUT / \"cache\"):\n    z = np.load(Path(cache_dir) / f\"{study}.npz\")\n    return {c: z[c] for c in CELL_NAMES}, dict(zip(CELL_NAMES, z[\"avail\"]))\n\ndef preprocess_test_study(study, series_df, root, cfg=CFG):\n    rows = series_df[series_df[\"StudyInstanceUID\"] == study].copy()\n    rows[\"n_slices\"] = [count_slices(root, study, s) for s in rows[\"SeriesInstanceUID\"]]\n    rows[\"cell\"] = rows.apply(cell_key, axis=1)\n    rows = rows[rows[\"cell\"].notna() & (rows[\"n_slices\"] > 0)]\n    rows[\"score\"] = -np.abs(np.log(rows[\"n_slices\"].clip(lower=1) / 30.0))\n    rows = rows.sort_values([\"cell\",\"score\",\"SeriesInstanceUID\"], ascending=[True, False, True])\n    chosen = dict(zip(rows.groupby(\"cell\").first().index,\n                      rows.groupby(\"cell\").first()[\"SeriesInstanceUID\"]))\n    tag, cx = study_lat_probe(root, study, list(chosen.values()))\n    if tag in (\"L\", \"R\"):\n        side, src = tag, \"tag\"\n    elif cx == cx and cfg[\"LAT_GEOM_THRESH\"] is not None:\n        side, src = (\"L\" if cx > cfg[\"LAT_GEOM_THRESH\"] else \"R\"), \"geometry\"\n    else:\n        side, src = \"U\", \"unknown\"\n    _, out, meta, have, errs = process_study((study, chosen, side, root, cfg))\n    return out, have, dict(laterality=side, lat_source=src, errors=errs)\n\n# smoke-test on the three example test studies\nROOT_TE = os.path.join(DATA, \"test_series\")\nfor s in test[\"StudyInstanceUID\"]:\n    t0 = time.time()\n    out, have, info = preprocess_test_study(s, tests, ROOT_TE)\n    print(f\"{s[:24]}...  {time.time()-t0:.2f}s  side={info['laterality']} ({info['lat_source']})  \"\n          f\"avail={have}  errors={info['errors']}\")\n\nout, have, info = preprocess_test_study(test['StudyInstanceUID'].iloc[0], tests, ROOT_TE)\nfig, ax = plt.subplots(1, len(CELL_NAMES), figsize=(3.0*len(CELL_NAMES), 3.3))\nfor a, c in zip(ax, CELL_NAMES):\n    v = out[c]; a.imshow(v[v.shape[0]//2], cmap=\"gray\", vmin=0, vmax=255); a.axis(\"off\")\n    a.set_title(f\"{c} ({have[c]})\", fontsize=8)\nfig.suptitle(\"Test-time path, identical transform\", fontsize=10)\nplt.tight_layout(); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"34f5d42d-d07b-451c-a0ce-0951de227c7f","cell_type":"markdown","source":"---\n\n## Notes for the next stage\n\n**What this cache is, in one line:** every study is now five volumes of a right knee at 0.833 mm/pixel over a\n160 mm field of view, slices ordered medial→lateral (sagittal), anterior→posterior (coronal),\ninferior→superior (axial), intensity-normalised per series to `uint8`.\n\n**The unresolved problem is labels, not images.** 58 gold studies against 12 targets — for MCL that is\n9 positives — will not train anything on its own, and the AUC is macro-averaged so the rare classes dominate\nthe score variance. The realistic path is:\n\n1. Train a **multilingual text classifier** on the reports, calibrated against the 58 gold studies, and use it\n   to produce soft labels for all 4 407. The EDA showed 64 % of reports are not in English, so an\n   XLM-R / multilingual-e5 backbone is the natural choice; the regex ruleset in the EDA notebook is a baseline\n   to beat, not a solution (recall 0.00 on both OA compartment labels).\n2. Train the image model on those soft labels, holding the 58 gold studies out entirely as the only honest\n   validation signal. Group any CV split by `PatientID` — it is present on 100 % of files and some patients\n   almost certainly contributed multiple studies.\n3. Treat the label noise explicitly: soft targets, and consider a two-stage scheme where the text model's\n   confidence weights the image loss.\n\n**Two things to verify before the first submission:**\n\n* the JPEG/JPEG2000 codec situation from §0.1 — the training data appears to be entirely uncompressed, so this\n  will only surface in the hidden test set, which is the worst place to discover it;\n* the `avail` mask semantics — a cell filled from its same-plane partner (value 2) is not the same as a real\n  observation, and the model should be given that flag rather than being silently fed a duplicate.\n\n**If you want more spatial resolution later**, the knobs are `OUT_SIZE` (192 → 256 costs 1.8× the cache) and\nthe per-cell `n_slices`. Everything else in the transform is scale-free.","metadata":{}}]}