{"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":"e3dbd9f3-fd7c-4f94-b743-a17064d3068c","cell_type":"markdown","source":"# RSNA Knee — Preprocessing (224 px)\n\nYour original preprocessing pipeline, unchanged except for the output resolution.\n\n**The only change: `OUT_SIZE` 192 → 224.**\n\n```\nv1:  160 mm / 192 px = 0.833 mm/pixel\nnow: 160 mm / 224 px = 0.714 mm/pixel\n```\n\n224 is chosen because it is divisible by both **16** (DINOv3 patch size) and **14** (DINOv2), so the same\ncache feeds either backbone with no resampling inside the model.\n\nEverything else is identical: the same five cells, the same 160 mm crop, the same slice ordering by\nprojection onto the slice normal, the same laterality resolution (tag → report → geometry), the same\nper-series intensity normalisation, the same fallback rules.\n\n### One thing to watch: cache size\n\n```\n4407 studies x 104 slices x 224^2 = 23.0 GB uncompressed   (v1 at 192 px: 16.9 GB)\n```\n\n`np.savez_compressed` typically brings that to 13–16 GB on disk, which fits `/kaggle/working`. §7 projects\nthe real number from a pilot **before** the full run and stops if it would overflow — a cache that\novershoots kills the kernel rather than slowing it down.\n\nIf the gate trips, lower the per-cell slice counts in `CFG[\"CELLS\"]` (the 4th field of each tuple).\n\n### A note on slice counts, since you train with 21\n\nThe cells hold 24 / 24 / 20 / 20 / 16 slices. Your training notebook takes `N_SLICES_PER_VIEW = 21` centred\nslices, so for `COR_FS`, `COR_T1` (20) and `AX_FS` (16) the index picker clamps at the stack edges and\nrepeats the outermost slice to reach 21. That is your existing v1 behaviour and it is preserved here\nunchanged — but if you want 21 genuinely distinct slices everywhere, raise those three counts to 24 (cache\ngrows to ~26.5 GB raw, so check the gate).\n","metadata":{}},{"id":"bffc724c-9967-45bf-84c0-d74162fccbf0","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\nfrom concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor\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\": 224,          # CHANGED from 192. 160mm/224px = 0.714 mm/px\n    \"FOV_MM\": 160.0,\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, 21, 76.0, \"COR_T1\"),\n        (\"COR_T1\", \"Coronal\",  0, 0, 21, 76.0, \"COR_FS\"),\n        (\"AX_FS\",  \"Axial\",    1, 1, 21, 96.0, None),\n    ],\n    # ---- intensity ------------------------------------------------------------\n    \"CLIP_LO_PCT\": 0.5,\n    \"CLIP_HI_PCT\": 99.5,\n    \"FG_THRESH_FRAC\": 0.06,\n    # ---- laterality -----------------------------------------------------------\n    \"CANONICAL_SIDE\": \"R\",\n    \"LAT_GEOM_THRESH\": None,\n    # ---- run control ----------------------------------------------------------\n    \"N_WORKERS\": 4,\n    \"PILOT_N\": 40,\n    \"COMPRESS\": True,\n    \"CACHE_GB_BUDGET\": 18.0,\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), \"  (v1 at 192px: 0.8333)\")\nprint(\"cells         :\", CELL_NAMES)\nprint(\"slices/study  :\", sum(c[4] for c in CFG[\"CELLS\"]))\nprint(\"raw bytes/study:\", f\"{sum(c[4] for c in CFG['CELLS']) * CFG['OUT_SIZE']**2 / 1e6:.2f} MB\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"612a5c65-cd8d-45de-a8b4-a351259c3df6","cell_type":"markdown","source":"## 1. Decoder check\n\nTraining data appears to be entirely uncompressed, so a missing JPEG/JPEG2000 codec can only surface in the\nhidden test set — the worst possible place to discover it.\n","metadata":{}},{"id":"0bf3363c-864e-433b-ba09-2b91ca3cf08c","cell_type":"code","source":"def decoder_status():\n    st = {}\n    for m in [\"pylibjpeg\", \"libjpeg\", \"openjpeg\", \"gdcm\"]:\n        try: __import__(m); st[m] = True\n        except Exception: st[m] = False\n    return st\n\nDEC = decoder_status(); print(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 still work so this notebook will finish, but compressed series\")\n    print(\"in the hidden test set will fail. Attach codec wheels before submitting:\")\n    print(\"  pip install --no-index --find-links=<dataset> pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg\")\n    print(\"!\"*78)\n","metadata":{},"outputs":[],"execution_count":null},{"id":"a04519ba-78e0-4213-89ed-29fda1bc7a93","cell_type":"markdown","source":"## 2. Tables and series selection\n\nSelection ignores `SeriesDescription` (~21% are `DUMMYSERIESDESC!` or blank) and ranks on slice count in log\nspace: a knee series has ~30 slices, so this demotes localisers and dynamic acquisitions without hard-coded\nthresholds. Ties break on the smallest UID so the choice is reproducible.\n","metadata":{}},{"id":"0c2d68fe-5dde-412c-a459-51ec958fe4c4","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    try: return sum(1 for f in os.scandir(os.path.join(root, study, series))\n                    if f.name.endswith(\".dcm\"))\n    except FileNotFoundError: return 0\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; return out\n\nROOT_TR = os.path.join(DATA, \"train_series\")\nt0 = time.time(); tser = add_slice_counts(tser, ROOT_TR)\nprint(f\"slice counts in {time.time()-t0:.0f}s | empty series: {(tser['n_slices']==0).sum()}\")\nprint(tser[\"n_slices\"].describe().round(1).to_string())\n","metadata":{},"outputs":[],"execution_count":null},{"id":"673500c3-3769-4b36-aa5f-6d63d00aca06","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 \\\n           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    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    return df, df.groupby([\"StudyInstanceUID\",\"cell\"], as_index=False).first()\n\nall_cells, SEL = build_series_index(tser)\nprint(\"chosen rows:\", len(SEL))\npiv = SEL.pivot_table(index=\"StudyInstanceUID\", columns=\"cell\",\n                      values=\"n_slices\", aggfunc=\"first\").reindex(columns=CELL_NAMES)\nprint(\"\\ncell availability after selection:\")\nfor c in CELL_NAMES:\n    print(f\"  {c:<10} {100*piv[c].notna().mean():5.1f} %\")\nprint(\"\\ncells per study:\")\nprint(piv.notna().sum(axis=1).value_counts().sort_index().to_string())\n","metadata":{},"outputs":[],"execution_count":null},{"id":"2325637d-5717-45b9-ab89-7c1f1772a8a5","cell_type":"markdown","source":"## 3. Reports and laterality\n\nFour of the twelve labels are medial/lateral pairs and a fifth (MCL) is named for its side, so mirroring\nevery knee into one convention is the transform that matters most. Resolution order: DICOM `Laterality` tag\n→ report text → geometry, with the geometry threshold **fitted and accuracy-reported** on the tagged half.\n","metadata":{}},{"id":"fe253f02-3d80-4886-b608-b54d6e493874","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]+\"); NL3 = re.compile(r\"\\n{3,}\")\n\ndef clean_report(t):\n    if not isinstance(t, str): return \"\"\n    t = unicodedata.normalize(\"NFKC\", t); t = PLACEHOLDER.sub(\" \", t); t = t.replace(\"\\r\", \"\\n\")\n    return NL3.sub(\"\\n\\n\", WS.sub(\" \", 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 = {\"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\"}}\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\nLAT_PAT = {\"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|левый|левая|слева|αριστερ)\"}\ndef lat_from_report(t):\n    tl = t.lower()\n    r_, l_ = len(re.findall(LAT_PAT[\"R\"], tl)), len(re.findall(LAT_PAT[\"L\"], tl))\n    return \"R\" if r_ > l_ else (\"L\" if l_ > r_ else 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: REP[\"PatientSex\"] = train[\"PatientSex\"]\nprint(REP[\"lang\"].value_counts().to_string())\nprint(\"\\nreport laterality found:\", round(REP[\"lat_report\"].notna().mean(), 3))\n","metadata":{},"outputs":[],"execution_count":null},{"id":"084ff750-c4a4-45cc-bb71-0f9577516326","cell_type":"code","source":"def study_lat_probe(root, study, series_ids):\n    tags, xs = [], []\n    for se in list(series_ids)[:3]:\n        d = os.path.join(root, study, se)\n        try: files = sorted(f.name for f in os.scandir(d) if f.name.endswith(\".dcm\"))\n        except FileNotFoundError: continue\n        if not files: continue\n        for f in [files[0], files[len(files)//2]]:\n            try: ds = pydicom.dcmread(os.path.join(d, f), stop_before_pixels=True, force=True)\n            except Exception: 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, rows = getattr(ds, \"Columns\", None), 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 = ipp + iop[:3]*(float(cols)/2*ps[1]) + iop[3:]*(float(rows)/2*ps[0])\n            xs.append(float(centre[0]))\n    return (Counter(tags).most_common(1)[0][0] if tags else None,\n            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()\nt0 = time.time()\nwith ThreadPoolExecutor(max_workers=12) as ex:\n    res = list(ex.map(lambda s: study_lat_probe(ROOT_TR, s, ser_by_study[s]), studies))\nprint(f\"probed {len(studies)} studies in {time.time()-t0:.0f}s\")\n\nLAT = pd.DataFrame({\"StudyInstanceUID\": studies,\n                    \"lat_tag\": [r[0] for r in res], \"centre_x\": [r[1] for r in res]})\nLAT = LAT.merge(REP[[\"StudyInstanceUID\",\"lat_report\"]], on=\"StudyInstanceUID\", how=\"left\")\nprint(\"tag available   :\", round(LAT[\"lat_tag\"].notna().mean(), 3))\nprint(\"report available:\", round(LAT[\"lat_report\"].notna().mean(), 3))\n\nboth = LAT.dropna(subset=[\"lat_tag\",\"lat_report\"])\nprint(f\"\\nreport vs tag agreement: {(both['lat_tag']==both['lat_report']).mean()*100:.1f}% (n={len(both)})\")\n\ng = LAT.dropna(subset=[\"lat_tag\",\"centre_x\"])\ncand = np.arange(-200, 100, 1.0)\nbest_acc, best_t = max(((np.where(g[\"centre_x\"].values > t, \"L\", \"R\") == g[\"lat_tag\"].values).mean(), t)\n                       for t in cand)\nCFG[\"LAT_GEOM_THRESH\"] = float(best_t)\nprint(f\"geometry rule: centre_x > {best_t:.1f} mm => LEFT | accuracy {best_acc*100:.2f}% (n={len(g)})\")\n\nplt.figure(figsize=(8, 3.2))\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()\n\nGEOM_MIN_ACC = 0.95\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\nr_ = LAT.apply(resolve_lat, axis=1, result_type=\"expand\")\nLAT[\"laterality\"], LAT[\"lat_source\"] = r_[0], r_[1]\nprint(\"\\n\" + LAT[\"lat_source\"].value_counts().to_string())\nprint(LAT[\"laterality\"].value_counts().to_string())\nprint(f\"studies to be mirrored: {(LAT['laterality'] == ('L' if CFG['CANONICAL_SIDE']=='R' else 'R')).sum()}\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"80b58e11-01fe-4995-a08b-6d32644f5bb2","cell_type":"markdown","source":"## 4. The transform\n\nUnchanged from v1. Slices are ordered by projection onto the slice normal — never by filename, since the\nfilename is a SOP Instance UID assigned to be unique rather than ordered, and your v1 EDA measured\n`InstanceNumber` agreeing with geometric order only 63.8% of the time.\n","metadata":{}},{"id":"5d8f01d4-0b60-4954-bdfb-32e411908fb4","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, iop=np.array(iop, float), ipp=np.array(ipp, float),\n                ps=np.array(ps, float), 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\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: a = a * g[\"slope\"] + g[\"inter\"]\n    if g[\"photo\"].upper() == \"MONOCHROME1\": 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    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: 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]\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]); 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); 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: 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    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    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), 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, short_stack=bool(span < ext))\n    return vol, meta\n","metadata":{},"outputs":[],"execution_count":null},{"id":"de1128d3-05f7-4766-9629-816158e38236","cell_type":"code","source":"def process_study(args):\n    study, rows, laterality, root, cfg = args\n    out, meta, errs, 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], meta[cell], have[cell] = v, m, 1\n        except Exception as e:\n            errs[cell] = f\"{type(e).__name__}: {e}\"[:160]\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        # `have[fb] == 1` (not `fb in out`): a cell zero-filled earlier in THIS loop is also\n        # \"in out\", and copying from it would flag an all-zero cell as available -- and would\n        # KeyError on meta[fb], which is only written in the try block above.\n        if fb and have.get(fb) == 1 and out[fb].shape[0] == spec[\"n_slices\"]:\n            out[cell] = out[fb].copy(); have[cell] = 2\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\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\nlat_map = dict(zip(LAT[\"StudyInstanceUID\"], LAT[\"laterality\"]))\ngrp = SEL.groupby(\"StudyInstanceUID\").apply(\n    lambda d: dict(zip(d[\"cell\"], d[\"SeriesInstanceUID\"]))).to_dict()\nJOBS = [(s, rows, lat_map.get(s, \"U\"), ROOT_TR, CFG) for s, rows in grp.items()]\nprint(\"jobs:\", len(JOBS))\n","metadata":{},"outputs":[],"execution_count":null},{"id":"07615dbb-c9f0-4c38-a0fb-f73b99a16309","cell_type":"markdown","source":"## 5. Pilot and size gate\n","metadata":{}},{"id":"42eb9084-bdbc-4811-9cd3-826bf3e45b18","cell_type":"code","source":"pilot = JOBS[:CFG[\"PILOT_N\"]]\nt0 = time.time(); sizes = []; shown = {}; pilot_err = []\nfor j in pilot:\n    st, out, meta, have, errs = process_study(j)\n    sizes.append(save_study(st, out, have))\n    for c, e in errs.items(): pilot_err.append((st, c, e))\n    if len(shown) < 2: shown[st] = (out, have)\ndt = time.time() - t0\nproj_gb = float(np.mean(sizes))*len(JOBS)/1e9\nprint(f\"{len(pilot)} studies in {dt:.1f}s -> {dt/len(pilot)*1000:.0f} ms/study\")\nprint(f\"mean npz {np.mean(sizes)/1e6:.2f} MB | projected cache {proj_gb:.1f} GB \"\n      f\"(budget {CFG['CACHE_GB_BUDGET']} 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[:8]]\nelse:\n    print(\"\\nno errors in pilot\")\n\nif proj_gb > CFG[\"CACHE_GB_BUDGET\"]:\n    raise RuntimeError(\n        f\"projected cache {proj_gb:.1f} GB exceeds the {CFG['CACHE_GB_BUDGET']} GB budget. \"\n        f\"Lower the per-cell slice counts in CFG['CELLS'] and re-run. Overshooting \"\n        f\"/kaggle/working kills the kernel rather than slowing it.\")\nprint(\"size gate passed\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"7e16dfe5-c223-47d6-9632-5061383af8ca","cell_type":"code","source":"st, (out, have) = list(shown.items())[0]\nfig, ax = plt.subplots(2, len(CELL_NAMES), figsize=(3.0*len(CELL_NAMES), 6.2))\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\")\nlh = LAT.set_index(\"StudyInstanceUID\").loc[st]\nfig.suptitle(f\"{st[:24]}...  laterality={lh['laterality']} ({lh['lat_source']})  |  \"\n             f\"{CFG['OUT_SIZE']}px, {CFG['FOV_MM']/CFG['OUT_SIZE']:.3f} mm/px\", fontsize=10)\nplt.tight_layout(); plt.show()\n","metadata":{},"outputs":[],"execution_count":null},{"id":"9f342e2e-c60e-4543-99e1-e62af3076011","cell_type":"code","source":"# Mirroring test: mean coronal of originally-left vs originally-right knees.\n# If the convention is inverted these come out as mirror images of each other.\nacc = {\"L\": [], \"R\": []}\nfor j in JOBS[:120]:\n    s = j[0]; side = lat_map.get(s, \"U\")\n    if side not in acc or len(acc[side]) >= 12: continue\n    if \"COR_FS\" not in j[1]: 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\nif acc[\"L\"] and acc[\"R\"]:\n    mL, mR = np.mean(acc[\"L\"], 0), np.mean(acc[\"R\"], 0)\n    fig, ax = plt.subplots(1, 3, figsize=(11, 4))\n    ax[0].imshow(mL, cmap=\"gray\"); ax[0].set_title(f\"mean COR_FS, originally LEFT (n={len(acc['L'])})\"); ax[0].axis(\"off\")\n    ax[1].imshow(mR, cmap=\"gray\"); ax[1].set_title(f\"mean COR_FS, originally RIGHT (n={len(acc['R'])})\"); ax[1].axis(\"off\")\n    ax[2].imshow(np.abs(mL-mR), cmap=\"magma\"); ax[2].set_title(\"|difference|\"); ax[2].axis(\"off\")\n    plt.tight_layout(); plt.show()\n    d_same, d_flip = float(np.abs(mL-mR).mean()), float(np.abs(mL-mR[:, ::-1]).mean())\n    print(f\"mean |L-R| = {d_same:.2f} | mean |L-flip(R)| = {d_flip:.2f}\")\n    print(\"=> mirroring is CORRECT\" if d_same < d_flip else\n          \"=> WARNING: flipping R matches better, the convention is inverted\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"c8c83b40-24f8-4d7d-9a1e-7d4bd9146f58","cell_type":"markdown","source":"## 6. Full run\n","metadata":{}},{"id":"798390c1-63b7-4e4d-9a4f-f6fadce3f3b2","cell_type":"code","source":"todo = [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_cell, 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: sz = save_study(st, out, have)\n        except Exception as e:\n            failures.append((st, \"save\", str(e)[:120])); 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(): rows_cell.append(dict(StudyInstanceUID=st, **m))\n        for c, e in errs.items(): 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\")\nprint(f\"\\ndone in {(time.time()-t0)/60:.1f} min | failures: {len(failures)}\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"af9764dc-4a21-40eb-9ec2-d1592c7687c3","cell_type":"code","source":"META = pd.DataFrame(rows_meta); CELLM = pd.DataFrame(rows_cell)\n\nSTUDY = LAT[[\"StudyInstanceUID\",\"laterality\",\"lat_source\",\"lat_tag\",\"lat_report\",\"centre_x\"]].merge(\n    META, on=\"StudyInstanceUID\", how=\"right\")\nSTUDY = STUDY.merge(REP[[\"StudyInstanceUID\",\"lang\",\"n_chars\"]], on=\"StudyInstanceUID\", how=\"left\")\nif \"PatientSex\" in train.columns:\n    STUDY = STUDY.merge(train[[\"StudyInstanceUID\",\"PatientSex\"]], on=\"StudyInstanceUID\", how=\"left\")\n\nlab_mask = train[LABELS].notna().all(axis=1)\nSTUDY = STUDY.merge(train.loc[lab_mask, [\"StudyInstanceUID\"]+LABELS],\n                    on=\"StudyInstanceUID\", how=\"left\")\nSTUDY[\"has_gold\"] = STUDY[LABELS].notna().all(axis=1).astype(int)\n\nprint(\"cached studies:\", len(STUDY), \"| with gold labels:\", int(STUDY[\"has_gold\"].sum()))\nprint(\"\\ncell availability (0=absent, 1=real, 2=fallback-filled):\")\nprint(STUDY[[f\"avail_{c}\" for c in CELL_NAMES]].apply(pd.Series.value_counts).fillna(0).astype(int))\nprint(f\"\\ntotal cache: {STUDY['npz_bytes'].sum()/1e9:.2f} GB \"\n      f\"({STUDY['npz_bytes'].mean()/1e6:.2f} MB/study)\")\n\nSTUDY.to_parquet(OUT/\"study_meta.parquet\", index=False)\nREP.to_parquet(OUT/\"reports_clean.parquet\", index=False)\nCELLM.to_parquet(OUT/\"series_index.parquet\", index=False)\nif failures:\n    pd.DataFrame(failures, columns=[\"StudyInstanceUID\",\"cell\",\"error\"]).to_csv(OUT/\"failures.csv\", index=False)\n\nwith open(OUT/\"preprocess_config.json\", \"w\") as f:\n    json.dump({**CFG, \"cells\": CELL_NAMES, \"mm_per_px\": CFG[\"FOV_MM\"]/CFG[\"OUT_SIZE\"],\n               \"geom_thresh_accuracy\": float(best_acc),\n               \"generated\": time.strftime(\"%Y-%m-%d %H:%M\")}, f, indent=2, default=str)\nprint(\"\\nwrote:\", [p.name for p in sorted(OUT.iterdir())])\n","metadata":{},"outputs":[],"execution_count":null},{"id":"db238651-8477-4003-a303-c485cec81568","cell_type":"markdown","source":"---\n\n**Check before moving on**\n\n1. §5's size gate passed and the total in §6 is comfortably under 20 GB.\n2. §5's mirroring test prints `mirroring is CORRECT`.\n3. Cell availability looks like your previous run — `COR_T1` will still show a large number of\n   fallback-filled (`avail = 2`) entries, which is expected.\n\nSave this notebook's output and attach it to the training notebook. `preprocess_config.json` now records\n`OUT_SIZE = 224`, and both downstream notebooks read the geometry from it rather than hard-coding it.\n","metadata":{}}]}