{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"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":[{"cell_type":"markdown","source":"# RSNA Knee -- Stage B-F: DICOM reconstruction through DINOv2 feature caching\n\n**This notebook produces the DINOv2 feature cache; it does not train anything.** Training lives in\na separate notebook, `v1_train.ipynb`, which takes this notebook's published output (the\n`feature_cache/` directory, published as a Kaggle Dataset) as its own input. The split exists\nbecause a full Stage F extraction pass takes multiple commit rounds across multiple sessions (see\nthe Stage F section, far below) -- running the training loop in the same notebook risked it\nexecuting against an incomplete cache and crashing partway through a long commit.\n\n## Stage B: reconstruct the real 3D volume from scrambled 2D DICOM files\n\nEach series under `train_series/<StudyInstanceUID>/<SeriesInstanceUID>/` is a folder of `.dcm` files whose\nfilenames carry no information about physical slice order. This stage turns that folder into a correctly\nordered 3D array plus the geometry needed to map any voxel back to a real (Left, Posterior, Superior)\ncoordinate, using:\n\n- `ImageOrientationPatient` (shared across the series) -> row/column direction vectors `r`, `c`\n- `n = r x c` -> the stacking axis (perpendicular to the imaging plane)\n- `ImagePositionPatient` (per slice) projected onto `n` -> the true spatial order to sort slices by\n\nStage B itself reconstructs at runtime, on every access, with no caching to disk -- no `.npz` files\nare written here. Given a `(StudyInstanceUID, SeriesInstanceUID)`, `reconstruct_series` reads the\nraw `.dcm` files and returns the ordered volume directly, every time it's called. (Stage F, far\nbelow, does cache something -- DINOv2's *output features*, not raw pixels -- see that section for\nwhy the two are treated differently.)","metadata":{}},{"cell_type":"code","source":"!pip install -q pydicom pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg python-gdcm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:04.457124Z","iopub.execute_input":"2026-08-25T17:07:04.457355Z","iopub.status.idle":"2026-08-25T17:07:10.857459Z","shell.execute_reply.started":"2026-08-25T17:07:04.457281Z","shell.execute_reply":"2026-08-25T17:07:10.856632Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nfrom pathlib import Path\n\n# OpenCV spawns its own internal thread pool per process. Combined with multiple\n# DataLoader worker processes (see the training cell), that causes CPU oversubscription\n# and can slow things down rather than speed them up -- let DataLoader's process-level\n# parallelism be the only parallelism, and disable cv2's own threading.\ncv2.setNumThreads(0)\n\nROOT_DIR = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection/train_series\")\nMANIFEST_PATH = Path(\"/kaggle/input/datasets/naveedlihazi/rsna-datasets/stage_a_manifest.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:10.860969Z","iopub.execute_input":"2026-08-25T17:07:10.861703Z","iopub.status.idle":"2026-08-25T17:07:12.311814Z","shell.execute_reply.started":"2026-08-25T17:07:10.861661Z","shell.execute_reply":"2026-08-25T17:07:12.311115Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## `reconstruct_series` — the Stage B function\n\nReturns a dict with the ordered `volume` plus everything needed to map a voxel index back to a real-world\nposition: `r`, `c` (in-plane direction vectors), `pixel_spacing` (mm/pixel), the sorted per-slice `ipp`\nanchors, and `rescale` (`RescaleSlope`/`RescaleIntercept`).\n\n`resize_to` is optional (`None` keeps native resolution). If set, `pixel_spacing` is rescaled by the same\nfactor as the resize — resizing changes how many millimeters each pixel represents, so the stored spacing\nhas to stay consistent with the array's actual dimensions for any later geometry (e.g. cross-plane\nalignment) to be correct.","metadata":{}},{"cell_type":"code","source":"def reconstruct_series(series_dir: Path, resize_to: int | None = None) -> dict:\n    dcm_paths = sorted(series_dir.glob(\"*.dcm\"))\n    if not dcm_paths:\n        raise ValueError(f\"no .dcm files found in {series_dir}\")\n\n    datasets = [pydicom.dcmread(p) for p in dcm_paths]\n\n    # orientation is shared across the whole series -> read once\n    iop = np.array(datasets[0].ImageOrientationPatient, dtype=np.float64)\n    r, c = iop[:3], iop[3:]\n    n = np.cross(r, c)\n    pixel_spacing = np.array(datasets[0].PixelSpacing, dtype=np.float64)  # [row_spacing, col_spacing]\n    slope = float(getattr(datasets[0], \"RescaleSlope\", 1.0))\n    intercept = float(getattr(datasets[0], \"RescaleIntercept\", 0.0))\n\n    ipp_list, pixel_list = [], []\n    for ds in datasets:\n        ipp_list.append(np.array(ds.ImagePositionPatient, dtype=np.float64))\n        arr = ds.pixel_array.astype(np.int16)\n        if resize_to is not None and arr.shape != (resize_to, resize_to):\n            row_scale = resize_to / arr.shape[0]\n            col_scale = resize_to / arr.shape[1]\n            arr = cv2.resize(arr, (resize_to, resize_to), interpolation=cv2.INTER_AREA)\n            pixel_spacing = pixel_spacing / np.array([row_scale, col_scale])\n        pixel_list.append(arr)\n\n    ipp_arr = np.stack(ipp_list)          # (depth, 3)\n    proj = ipp_arr @ n                    # scalar position along the stacking axis, per slice\n    order = np.argsort(proj)\n\n    volume = np.stack(pixel_list)[order]  # (depth, height, width), correctly ordered\n    ipp_sorted = ipp_arr[order]\n\n    return {\n        \"volume\": volume,\n        \"ipp\": ipp_sorted.astype(np.float32),\n        \"r\": r.astype(np.float32),\n        \"c\": c.astype(np.float32),\n        \"pixel_spacing\": pixel_spacing.astype(np.float32),\n        \"rescale\": np.array([slope, intercept], dtype=np.float32),\n    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:12.313434Z","iopub.execute_input":"2026-08-25T17:07:12.314987Z","iopub.status.idle":"2026-08-25T17:07:12.333987Z","shell.execute_reply.started":"2026-08-25T17:07:12.314847Z","shell.execute_reply":"2026-08-25T17:07:12.331998Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Stage C: locate the knee and crop away irrelevant background\n\nEvery reconstructed slice includes coil/background area outside the joint, and the acquired field\nof view (`Rows × PixelSpacing`) varies a lot across this corpus — a fixed-pixel crop would mean a\nwildly different amount of real anatomy per series. Two approaches were evaluated for finding the\ncrop region:\n\n1. **A learned localizer** — segment the knee with a general-purpose foundation model (MedSAM),\n   prompted with a box built from this series' own header data, then derive a crop from the mask.\n   Fully explored and validated end-to-end in `notebooks/working_notebook/stage_c_medsam_experiment.ipynb`.\n2. **A fixed physical-extent crop, centered on the frame** — no model, no extra inference cost:\n   convert a constant real-world window (`CROP_MM`) into pixels via this series' own\n   `PixelSpacing`, and crop around the array's geometric center.\n\n**Decision: (2), the fixed-mm center crop — the same approach as\n`notebooks/reference_notebooks/rsna-knee-baseline-v1.ipynb`.** This was not a default choice; it\nfollowed a real evaluation of (1) across six series spanning every plane/contrast slot, which\nfound: IoU confidence consistently moderate-to-low (0.37–0.56, never confident, in *every* slot\ntested, not one bad sample); derived crop sizes converging to 95–100% of the seed box's own size\nrather than a meaningfully tighter, independently-found boundary; one series where the derived\ncrop came out *larger than the source image itself* despite having the highest confidence of the\nsix; and no ground truth in this competition's data to check whether any of this was actually more\naccurate than the plain centered assumption, only whether it looked plausible. Full writeup and\nthe four-point evidence trail: `docs/progress_so_far.md`, \"Stage C: knee localization\" section.\nThis conclusion is about this corpus and this compute budget as evaluated — not a permanent\nban on revisiting a learned localizer if a first training run later shows an error pattern that\ncorrelates with off-center or unusual-FOV acquisitions.\n\n`CROP_MM = 130.0` matches the reference notebook's own calibration: the acquired field of view\nacross the corpus has median ~160mm and ranges 70–320mm, and 130mm sits below the field of view of\n~99.6% of series while still comfortably containing the joint. Row and column spacing are handled\n**separately** here (`s_row` sizes the height crop, `s_col` sizes the width crop) rather than\nassuming isotropic spacing, since `PixelSpacing`'s row/column order is easy to get backwards and\nthis project has already hit that gotcha once (in the MedSAM box-prompt work) — this is a small,\ndeliberate refinement over blindly porting the reference notebook's exact line, not a departure\nfrom its logic.","metadata":{}},{"cell_type":"code","source":"CROP_MM = 130.0\n\n\ndef crop_to_knee(volume: np.ndarray, pixel_spacing: np.ndarray, crop_mm: float = CROP_MM) -> np.ndarray:\n    depth, H, W = volume.shape\n    s_row, s_col = pixel_spacing\n\n    want_h = int(round(crop_mm / s_row))\n    want_w = int(round(crop_mm / s_col))\n\n    # If the requested window isn't smaller than this series' actual field of view, cropping\n    # would either do nothing or ask for more than exists — leave the volume untouched rather\n    # than silently misbehave.\n    if not (16 < want_h < H and 16 < want_w < W):\n        return volume\n\n    cy, cx = H // 2, W // 2\n    half_h, half_w = want_h // 2, want_w // 2\n\n    y1, y2 = max(0, cy - half_h), min(H, cy + half_h)\n    x1, x2 = max(0, cx - half_w), min(W, cx + half_w)\n    return volume[:, y1:y2, x1:x2]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:12.33537Z","iopub.execute_input":"2026-08-25T17:07:12.33573Z","iopub.status.idle":"2026-08-25T17:07:12.350057Z","shell.execute_reply.started":"2026-08-25T17:07:12.335706Z","shell.execute_reply":"2026-08-25T17:07:12.349452Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Sanity check on one real series\n\nGrabs the first `filled` row from the Stage A manifest, reconstructs it (Stage B), then crops it\nto the knee region (Stage C). Prints shapes before/after the crop plus the consecutive-slice gap\ncheck described above (should be small and roughly consistent for a clean series).","metadata":{}},{"cell_type":"code","source":"manifest = pd.read_csv(MANIFEST_PATH)\nfilled = manifest[manifest[\"status\"] == \"filled\"]\n\nsample = filled.iloc[0]\nseries_dir = ROOT_DIR / sample[\"StudyInstanceUID\"] / sample[\"SeriesInstanceUID\"]\n\nresult = reconstruct_series(series_dir)\n\nprint(\"Study:\", sample[\"StudyInstanceUID\"])\nprint(\"Series:\", sample[\"SeriesInstanceUID\"], \"| Slot:\", sample[\"Slot\"])\nprint(\"volume shape (pre-crop):\", result[\"volume\"].shape, result[\"volume\"].dtype)\nprint(\"r:\", result[\"r\"], \"c:\", result[\"c\"])\nprint(\"pixel_spacing:\", result[\"pixel_spacing\"])\nprint(\"rescale (slope, intercept):\", result[\"rescale\"])\n\nn = np.cross(result[\"r\"], result[\"c\"])\nproj = result[\"ipp\"] @ n\ngaps = np.diff(proj)\nprint(\"consecutive slice gaps (mm):\", gaps)\n\ncropped = crop_to_knee(result[\"volume\"], result[\"pixel_spacing\"])\nprint(\"\\nvolume shape (post-crop):\", cropped.shape, cropped.dtype)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:12.350928Z","iopub.execute_input":"2026-08-25T17:07:12.352051Z","iopub.status.idle":"2026-08-25T17:07:12.932624Z","shell.execute_reply.started":"2026-08-25T17:07:12.352027Z","shell.execute_reply":"2026-08-25T17:07:12.93171Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Stage E: synthetic cross-plane volumes (the paper's \"rotation\" trick)\n\nFor a given study, and for every ordered pair of *present* fluid planes (e.g. Sagittal → Coronal),\nrelabel the source volume's axes into the target plane's frame by comparing their `r`/`c`/`n`\ndirection vectors with dot products, then `transpose` + `flip` — no pixel values are computed or\ninvented, only relabeled.\n\nTwo things worth being explicit about:\n\n- **Rotation is only ever attempted when the target plane's own native volume already exists for\n  that study.** If a plane has no native series, there's no branch for it to feed anyway (see the\n  Study X walkthrough in `docs/progress_so_far.md`), so there's never a case where a target frame\n  is needed but unavailable.\n- **Spacing has to be relabeled along with the pixel data.** A volume's two in-plane axes have a\n  real `PixelSpacing`, but its depth axis doesn't — it has a slice *gap* instead (computed from\n  `ipp @ n`, the same projection Stage B already uses for ordering). Rotating axes means whichever\n  new axis absorbs the old depth axis needs to carry that slice-gap value forward as its\n  \"spacing,\" not a `PixelSpacing` number that was never measured for that direction.","metadata":{}},{"cell_type":"code","source":"def _axis_correspondence(src_axes: np.ndarray, tgt_axes: np.ndarray):\n    \"\"\"src_axes, tgt_axes: (3,3) arrays, each row a unit direction vector, ordered\n    [n, r, c] to match a volume's own (depth, height, width) axis order.\n\n    Returns (transpose_order, flip_axes) such that\n        np.flip(np.transpose(volume, transpose_order), flip_axes)\n    sits in the target frame.\n    \"\"\"\n    dot = src_axes @ tgt_axes.T                      # dot[i, j]: source axis i vs target axis j\n    src_for_tgt = np.argmax(np.abs(dot), axis=0)      # for target axis j, best-matching source axis\n    if len(set(src_for_tgt.tolist())) != 3:\n        raise ValueError(\"degenerate axis correspondence -- geometry too close to non-orthogonal\")\n\n    transpose_order = tuple(int(i) for i in src_for_tgt)\n    flip_axes = tuple(j for j in range(3) if dot[src_for_tgt[j], j] < 0)\n    return transpose_order, flip_axes\n\n\ndef rotate_to_frame(volume: np.ndarray, r_src, c_src, n_src, slice_gap_src, pixel_spacing_src,\n                     r_tgt, c_tgt, n_tgt):\n    \"\"\"volume: (depth, H, W), already cropped + resized. Returns the same pixel data\n    relabeled into the target plane's frame, plus the correspondingly relabeled\n    (depth-spacing, H-spacing, W-spacing) tuple. No pixel values are invented.\"\"\"\n    src_axes = np.stack([n_src, r_src, c_src]).astype(np.float64)\n    tgt_axes = np.stack([n_tgt, r_tgt, c_tgt]).astype(np.float64)\n    transpose_order, flip_axes = _axis_correspondence(src_axes, tgt_axes)\n\n    rotated = np.transpose(volume, transpose_order)\n    if flip_axes:\n        rotated = np.flip(rotated, axis=flip_axes)\n    rotated = np.ascontiguousarray(rotated)\n\n    spacing_src_full = np.array([slice_gap_src, pixel_spacing_src[0], pixel_spacing_src[1]])\n    rotated_spacing = spacing_src_full[list(transpose_order)]\n    return rotated, rotated_spacing","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:12.933827Z","iopub.execute_input":"2026-08-25T17:07:12.934481Z","iopub.status.idle":"2026-08-25T17:07:12.941796Z","shell.execute_reply.started":"2026-08-25T17:07:12.934456Z","shell.execute_reply":"2026-08-25T17:07:12.940913Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Normalization, resizing, and the frozen DINOv2 encoder\n\n`reconstruct_series` deliberately does **not** apply `RescaleSlope`/`RescaleIntercept` -- that has\nto happen explicitly before any intensity-sensitive use. The order settled on is: reconstruct at\nnative resolution -> apply rescale + percentile-normalize to `[0,1]` -> `crop_to_knee` (uses the\n*native* `PixelSpacing`, which is what the 130mm physical window is defined against) -> **then**\nresize the already-cropped result to `RESIZE_TO=224`. Resizing before cropping was considered and\nrejected -- it would make the crop math operate on already-downscaled pixels rather than the real\nacquired resolution the physical window was calibrated against.\n\n`RESIZE_TO=224` and the DINOv2 variant (`dinov2_vitb14_reg` -- Base, with the register-token fix,\n~86M params, 768-dim output) were both decided with the reasoning already recorded in this\nproject's history: 224 matches DINOv2's own pretraining resolution (avoiding positional-embedding\nmismatch on a *frozen* encoder); Base balances frozen-feature quality against the throughput cost\nof encoding every slice. Both are defaults to revisit with evidence, not final commitments -- and\nnow, once Stage F's cache is built against them, revisiting either means rebuilding the entire\ncache from scratch.\n\n`dino` is loaded here for Stage F's one-time feature-extraction pass below -- it never appears in\nthe training notebook (`v1_train.ipynb`) at all; that notebook consumes Stage F's cached output\ndirectly and has no DINOv2 dependency of its own.","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nFLUID_PLANES = [\"Sagittal\", \"Coronal\", \"Axial\"]\nRESIZE_TO = 224\n\nDEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n\n\ndef apply_rescale_and_normalize(volume: np.ndarray, rescale: np.ndarray) -> np.ndarray:\n    slope, intercept = rescale\n    vol = volume.astype(np.float32) * slope + intercept\n    lo, hi = np.percentile(vol, [1, 99])\n    return np.clip((vol - lo) / max(hi - lo, 1e-6), 0, 1)\n\n\ndef resize_volume(volume: np.ndarray, size: int = RESIZE_TO,\n                   interpolation: int = cv2.INTER_AREA) -> np.ndarray:\n    out = np.empty((volume.shape[0], size, size), dtype=np.float32)\n    for i in range(volume.shape[0]):\n        out[i] = cv2.resize(volume[i], (size, size), interpolation=interpolation)\n    return out\n\n\nIMAGENET_MEAN = torch.tensor([0.485, 0.456, 0.406]).view(1, 3, 1, 1)\nIMAGENET_STD = torch.tensor([0.229, 0.224, 0.225]).view(1, 3, 1, 1)\n\ndino = torch.hub.load(\"facebookresearch/dinov2\", \"dinov2_vitb14_reg\")\ndino.eval()\nfor p in dino.parameters():\n    p.requires_grad = False\ndino = dino.to(DEVICE)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:12.943978Z","iopub.execute_input":"2026-08-25T17:07:12.944356Z","iopub.status.idle":"2026-08-25T17:07:31.758463Z","shell.execute_reply.started":"2026-08-25T17:07:12.944279Z","shell.execute_reply":"2026-08-25T17:07:31.757781Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Stage F: caching DINOv2's output features (one-time, not per-epoch)\n\nDINOv2 is frozen and never updated by the optimizer -- for a fixed pixel input, it returns the\nsame 768-dim feature vector every time, on every epoch. But the pipeline as originally written\nrecomputed it from scratch on every single epoch (measured: ~21.7s/study of GPU time), which is\nwhy a real 12-epoch run over the full corpus extrapolated to ~319 hours -- far outside any\nrealistic Kaggle budget. Nothing about the *model* changes here; what changes is *when* DINOv2\nruns: once, ever, per `(SeriesInstanceUID, role)`, instead of once per study per epoch.\n\n**Cache key**: `(SeriesInstanceUID, role)`. A `SeriesInstanceUID` belongs to exactly one study,\nand within a study the rotation target reference is picked the same deterministic way every time\n(`tgt_data[\"native\"][0]`, same as the original per-epoch Dataset), so this key has no cross-study\nambiguity. Three roles:\n\n- `native` -- a fluid series' own plane, cropped/resized, no rotation\n- `structural` -- a structural series' own plane, cropped/resized, no rotation (structural\n  series are never rotated -- only native fluid series ever become rotation companions)\n- `rotated_into_<Plane>` -- a native series, rotated into some *other* present plane's frame,\n  as that plane's cross-plane-attention companion\n\n**Storage**: one `.npz` per study (not per series-role -- a study's ~13 entries belong together,\nand `KneeStudyDataset` already loads one study at a time), each array `float16`, shape\n`(n_slices, 768)`, keyed by `f\"{SeriesInstanceUID}__{role}\"` inside the archive.\n\nThis is exactly the same Stage B -> normalize -> Stage C -> resize -> Stage E pipeline the\noriginal per-epoch `KneeStudyDataset` ran live, every time -- `extract_study_features` below just\nruns it once and saves the DINOv2 output instead of throwing it away.","metadata":{}},{"cell_type":"code","source":"FEATURE_CACHE_DIR = Path(\"/kaggle/working/feature_cache\")\nFEATURE_CACHE_DIR.mkdir(parents=True, exist_ok=True)\n\nENCODE_CHUNK = 64  # max slices per single DINOv2 forward call -- bounds peak GPU memory\n                    # on the largest studies while still batching far fewer, larger calls\n                    # than one-call-per-series\n\n\ndef _load_prepped_series(study_uid: str, series_uid: str, root_dir: Path,\n                          crop_mm: float = CROP_MM, resize_to: int = RESIZE_TO) -> dict:\n    \"\"\"Stage B + normalize + Stage C + resize -- the shared pixel prep every role\n    (native/structural/rotated) starts from before diverging.\"\"\"\n    series_dir = root_dir / study_uid / series_uid\n    result = reconstruct_series(series_dir)\n    vol = apply_rescale_and_normalize(result[\"volume\"], result[\"rescale\"])\n    vol = crop_to_knee(vol, result[\"pixel_spacing\"], crop_mm)\n    vol = resize_volume(vol, resize_to)\n\n    n = np.cross(result[\"r\"], result[\"c\"])\n    proj = result[\"ipp\"] @ n\n    slice_gap = float(np.mean(np.diff(proj))) if len(proj) > 1 else 1.0\n\n    return {\"volume\": vol, \"r\": result[\"r\"], \"c\": result[\"c\"], \"n\": n,\n            \"pixel_spacing\": result[\"pixel_spacing\"], \"slice_gap\": slice_gap}\n\n\n@torch.no_grad()\ndef _encode_volumes(volumes: list[np.ndarray], dino: nn.Module, device) -> list[np.ndarray]:\n    \"\"\"One batched DINOv2 call for many volumes' worth of slices at once (chunked to\n    bound peak GPU memory), split back into one float16 array per volume. This is the\n    only place DINOv2 ever runs now -- once per (series, role), not once per epoch.\"\"\"\n    if not volumes:\n        return []\n    sizes = [v.shape[0] for v in volumes]\n    all_slices = np.concatenate(volumes, axis=0)\n\n    feats_chunks = []\n    for start in range(0, len(all_slices), ENCODE_CHUNK):\n        chunk = all_slices[start:start + ENCODE_CHUNK]\n        x = torch.from_numpy(chunk).unsqueeze(1).repeat(1, 3, 1, 1).float().to(device)\n        x = (x - IMAGENET_MEAN.to(device)) / IMAGENET_STD.to(device)\n        feats_chunks.append(dino(x))\n    feats = torch.cat(feats_chunks, dim=0).half().cpu().numpy()  # cast to fp16 only for storage\n\n    out, idx = [], 0\n    for s in sizes:\n        out.append(feats[idx:idx + s])\n        idx += s\n    return out\n\n\ndef extract_study_features(study_uid: str, manifest: pd.DataFrame, root_dir: Path,\n                            dino: nn.Module, device) -> dict:\n    \"\"\"Runs Stage B/C/E once per (series, role) needed for this study and encodes each\n    with DINOv2, returning {\"<SeriesInstanceUID>__<role>\": float16 array}. Exactly the\n    set of pixel volumes the original per-epoch Dataset used to reconstruct from scratch\n    on every single access.\"\"\"\n    rows = manifest[manifest[\"StudyInstanceUID\"] == study_uid]\n\n    fluid, structural = {}, {}\n    for plane in FLUID_PLANES:\n        fluid_rows = rows[(rows[\"Slot\"] == f\"{plane}-fluid\") & (rows[\"status\"] == \"filled\")]\n        if len(fluid_rows) == 0:\n            continue\n        fluid[plane] = [(r.SeriesInstanceUID,\n                          _load_prepped_series(study_uid, r.SeriesInstanceUID, root_dir))\n                         for r in fluid_rows.itertuples()]\n\n        struct_rows = rows[(rows[\"Slot\"] == f\"{plane}-structural\") & (rows[\"status\"] == \"filled\")]\n        structural[plane] = [(r.SeriesInstanceUID,\n                               _load_prepped_series(study_uid, r.SeriesInstanceUID, root_dir))\n                              for r in struct_rows.itertuples()]\n\n    keys, volumes = [], []\n\n    for plane, series in fluid.items():\n        for series_uid, s in series:\n            keys.append(f\"{series_uid}__native\")\n            volumes.append(s[\"volume\"])\n    for plane, series in structural.items():\n        for series_uid, s in series:\n            keys.append(f\"{series_uid}__structural\")\n            volumes.append(s[\"volume\"])\n\n    # every present plane's every native series becomes a rotation source into every\n    # OTHER present plane -- same rule, same reference-frame choice, as the original\n    # per-epoch Dataset used, just computed once instead of every epoch.\n    for tgt_plane, tgt_series in fluid.items():\n        tgt_ref = tgt_series[0][1]\n        for src_plane, src_series in fluid.items():\n            if src_plane == tgt_plane:\n                continue\n            for series_uid, s in src_series:\n                rotated_vol, _ = rotate_to_frame(\n                    s[\"volume\"], s[\"r\"], s[\"c\"], s[\"n\"], s[\"slice_gap\"], s[\"pixel_spacing\"],\n                    tgt_ref[\"r\"], tgt_ref[\"c\"], tgt_ref[\"n\"],\n                )\n                if rotated_vol.shape[1:] != (RESIZE_TO, RESIZE_TO):\n                    rotated_vol = resize_volume(rotated_vol, RESIZE_TO, interpolation=cv2.INTER_LINEAR)\n                keys.append(f\"{series_uid}__rotated_into_{tgt_plane}\")\n                volumes.append(rotated_vol)\n\n    feats = _encode_volumes(volumes, dino, device)\n    return dict(zip(keys, feats))\n\n\ndef cache_study(study_uid: str, manifest: pd.DataFrame, root_dir: Path, dino: nn.Module,\n                 device, cache_dir: Path = FEATURE_CACHE_DIR) -> str:\n    \"\"\"Wraps extract_study_features with the resumability check: skip a study whose\n    .npz already exists, catch and log a per-study failure rather than raising, so one\n    bad series can't take down a run that's already hours in.\"\"\"\n    out_path = cache_dir / f\"{study_uid}.npz\"\n    if out_path.exists():\n        return \"skipped\"\n    try:\n        feats = extract_study_features(study_uid, manifest, root_dir, dino, device)\n    except Exception as e:\n        print(f\"  FAILED {study_uid}: {e}\")\n        return \"failed\"\n    np.savez(out_path, **feats)\n    return \"done\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:31.759278Z","iopub.execute_input":"2026-08-25T17:07:31.759766Z","iopub.status.idle":"2026-08-25T17:07:31.776902Z","shell.execute_reply.started":"2026-08-25T17:07:31.759742Z","shell.execute_reply":"2026-08-25T17:07:31.776276Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Resuming from a previous session\n\n`/kaggle/working` starts **empty every session** -- it is not the persistent feature cache itself,\njust this session's scratch output. What actually persists across sessions is whatever gets\n**published as a Kaggle Dataset** after a commit finishes. That means `cache_study`'s\nskip-if-exists check, on its own, only ever sees studies cached *in the current session* -- it has\nnothing to skip on a fresh run, even if 4,000 studies were already cached in a previous one. This\ncell is what actually makes resumoption work: it copies a previously-published cache (if attached\nas input) back into `/kaggle/working/feature_cache` *before* the smoke test or driver run, so the\nskip-if-exists check has something real to check against.","metadata":{}},{"cell_type":"code","source":"import shutil\n\n# Update this once you've published a previous run's /kaggle/working/feature_cache as a\n# Kaggle Dataset and attached it as input to this session -- point it at wherever that\n# dataset's feature_cache/ folder lands under /kaggle/input/.\nPREVIOUS_CACHE_DIR = Path(\"/kaggle/input/datasets/naveedlihazi/dino-v2-extracted-features-cache/feature_cache\")\n\nif PREVIOUS_CACHE_DIR.exists():\n    resumed = 0\n    for f in PREVIOUS_CACHE_DIR.glob(\"*.npz\"):\n        shutil.copy2(f, FEATURE_CACHE_DIR / f.name)\n        resumed += 1\n    print(f\"resumed {resumed} already-cached studies from {PREVIOUS_CACHE_DIR}\")\nelse:\n    print(f\"no previous cache found at {PREVIOUS_CACHE_DIR} -- starting fresh\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:07:31.777896Z","iopub.execute_input":"2026-08-25T17:07:31.778169Z","iopub.status.idle":"2026-08-25T17:08:30.347963Z","shell.execute_reply.started":"2026-08-25T17:07:31.778138Z","shell.execute_reply":"2026-08-25T17:08:30.347283Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## The extraction driver -- resumable across multiple Kaggle sessions\n\nOnce the smoke test above looks correct, this is the real, resumable extraction pass. A full run\nover ~4,407 studies is still ~20+ hours of GPU time -- too long for one Kaggle session (~9-12h) or\neven one week's quota (~30h). `cache_study` makes this resumable rather than needing to fit in one\nrun: it skips any study whose `.npz` already exists in `FEATURE_CACHE_DIR`. The intended workflow,\nreusing the same pattern already proven necessary for Stage B's own disk-quota struggles: run a\nbounded chunk of studies -> commit -> publish `/kaggle/working/feature_cache` as a (new version of\na) Kaggle Dataset -> attach it as input in the next session -> the skip-if-exists check picks up\nwhere the last session left off -> repeat until the full corpus is cached.\n\n`TIME_BUDGET_HOURS` bounds a single run so it stops cleanly with time to spare before a commit's\ntime limit, rather than risking a mid-extraction kill that (per Kaggle's own persistence rules)\ncould lose a run that never finished committing -- **set this to something safely below your\nsession's actual time ceiling before running for real**, the `None` default has no self-imposed\nstop and will run until Kaggle kills it. A single bad series (corrupt DICOM, unexpected geometry)\nshouldn't take down a run that's already hours in, so failures are caught and logged per-study\nrather than raised.","metadata":{}},{"cell_type":"code","source":"import time\nfrom tqdm.auto import tqdm\n\nTIME_BUDGET_HOURS = 18.0  # e.g. 8.0 -- stop cleanly before a Kaggle commit's time limit.\n                           # None = no limit (fine for a quick manual/interactive run).\n\nall_study_uids = manifest[\"StudyInstanceUID\"].unique()\ncounts = {\"done\": 0, \"skipped\": 0, \"failed\": 0}\nstart = time.perf_counter()\n\nextract_bar = tqdm(all_study_uids, desc=\"Stage F: caching DINOv2 features\", unit=\"study\")\nfor study_uid in extract_bar:\n    if TIME_BUDGET_HOURS is not None and (time.perf_counter() - start) / 3600 > TIME_BUDGET_HOURS:\n        print(f\"\\nstopping early: {TIME_BUDGET_HOURS}h budget reached\")\n        break\n    result = cache_study(study_uid, manifest, ROOT_DIR, dino, DEVICE)\n    counts[result] += 1\n    extract_bar.set_postfix(**counts)\n\nprint(f\"\\ndone: {counts['done']}  skipped (already cached): {counts['skipped']}  failed: {counts['failed']}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-25T17:08:30.349001Z","iopub.execute_input":"2026-08-25T17:08:30.349359Z","iopub.status.idle":"2026-08-25T17:08:57.618934Z","shell.execute_reply.started":"2026-08-25T17:08:30.349326Z","shell.execute_reply":"2026-08-25T17:08:57.616554Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Next: run the driver above for real, then switch to `v1_train.ipynb`\n\nOnce the smoke test passed, run the driver cell above for real (after setting `TIME_BUDGET_HOURS`\nto a safe value below your session's actual commit time limit). After each commit finishes,\npublish `/kaggle/working/feature_cache` as a (new version of a) Kaggle Dataset, attach it as input\nto the next session, and rerun the driver -- the skip-if-cached check picks up automatically where\nthe last round left off. Once every study has a cached `.npz`, open `v1_train.ipynb`, point its\n`FEATURE_CACHE_DIR` at the published dataset, and run the training loop there.","metadata":{}}]}