{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":99552,"databundleVersionId":13851420,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# ============================================================\n# STEP 1: RSNA IA — Balanced CTA+MRA Subset Builder\n# Outputs: balanced_cta_mra_subset_50.csv\n# ============================================================\nimport os\nfrom pathlib import Path\nimport math\nimport json\nimport numpy as np\nimport pandas as pd\n\n# ---------------------- CONFIG ----------------------\nCONFIG = {\n    \"train_csv\": \"/kaggle/input/rsna-intracranial-aneurysm-detection/train.csv\",\n    \"localizers_csv\": \"/kaggle/input/rsna-intracranial-aneurysm-detection/train_localizers.csv\",\n    \"segmentations_dir\": \"/kaggle/input/rsna-intracranial-aneurysm-detection/segmentations\",\n    \"modalities_to_keep\": [\"CTA\", \"MRA\"],\n    \"fractions\": [0.20, 0.80],\n    \"seed\": 18,\n    \"multiframe_uids_csv\": None,\n    \"min_per_site_20\": 3,\n    \"min_per_site_50\": 6,\n}\n\n# Robust pathing\ndef _safe_path(p):\n    p = Path(p)\n    return p if p.exists() else Path(\"/kaggle/working\") / p.name\n\ntrain_csv = _safe_path(CONFIG[\"train_csv\"])\nlocalizers_csv = _safe_path(CONFIG[\"localizers_csv\"])\nsegs_dir = Path(CONFIG[\"segmentations_dir\"])\n\n# Load data\ntrain = pd.read_csv(train_csv)\nloc = pd.read_csv(localizers_csv)\n\nLABEL_COLS = [\n    'Left Infraclinoid Internal Carotid Artery', 'Right Infraclinoid Internal Carotid Artery',\n    'Left Supraclinoid Internal Carotid Artery', 'Right Supraclinoid Internal Carotid Artery',\n    'Left Middle Cerebral Artery', 'Right Middle Cerebral Artery',\n    'Anterior Communicating Artery', 'Left Anterior Cerebral Artery',\n    'Right Anterior Cerebral Artery', 'Left Posterior Communicating Artery',\n    'Right Posterior Communicating Artery', 'Basilar Tip', 'Other Posterior Circulation',\n]\nGLOBAL_COL = 'Aneurysm Present'\n\n# Keep only CTA/MRA\ntrain = train[train['Modality'].isin(CONFIG[\"modalities_to_keep\"])].copy().reset_index(drop=True)\n\n# Flags\ntrain['has_localizer'] = train['SeriesInstanceUID'].isin(loc['SeriesInstanceUID'])\nif segs_dir.exists():\n    seg_uids = set(p.stem for p in segs_dir.glob(\"*.nii*\"))\nelse:\n    seg_uids = set()\ntrain['has_segmentation'] = train['SeriesInstanceUID'].isin(seg_uids)\n\n# Helper: sample with priority\ndef _sample_class_with_priority(df_class: pd.DataFrame, n_target: int, seed: int):\n    if n_target <= 0:\n        return df_class.iloc[0:0].copy()\n\n    counts_by_mod = df_class['Modality'].value_counts()\n    total = counts_by_mod.sum()\n    desired = {m: int(round(n_target * (counts_by_mod[m] / total))) for m in counts_by_mod.index}\n\n    # Fix rounding drift\n    drift = n_target - sum(desired.values())\n    if drift != 0:\n        order = counts_by_mod.sort_values(ascending=False).index.tolist()\n        i = 0\n        while drift != 0 and i < len(order):\n            m = order[i]\n            desired[m] += 1 if drift > 0 else -1\n            drift += -1 if drift > 0 else 1\n            i = (i + 1) % len(order)\n\n    df_class = df_class.sort_values(['has_localizer', 'has_segmentation'], ascending=False)\n    out_parts = []\n    for m, want in desired.items():\n        block = df_class[df_class['Modality'] == m]\n        take = min(want, len(block))\n        out_parts.append(block.head(take))\n\n    sampled = pd.concat(out_parts).drop_duplicates(subset=['SeriesInstanceUID'])\n    if len(sampled) < n_target:\n        remain = df_class[~df_class['SeriesInstanceUID'].isin(sampled['SeriesInstanceUID'])]\n        topup = remain.head(n_target - len(sampled))\n        sampled = pd.concat([sampled, topup])\n\n    return sampled.sample(frac=1.0, random_state=seed).reset_index(drop=True)\n\ndef _enforce_min_site_coverage(df_pos: pd.DataFrame, min_per_site: int, seed: int):\n    if min_per_site <= 0 or df_pos.empty:\n        return df_pos\n\n    cov = {site: int(df_pos[site].sum()) for site in LABEL_COLS}\n    pool = df_pos.sort_values(['has_localizer', 'has_segmentation'], ascending=False).copy()\n    extra = []\n    seen = set(df_pos['SeriesInstanceUID'])\n\n    for site in LABEL_COLS:\n        need = max(0, min_per_site - cov[site])\n        if need == 0:\n            continue\n        cand = pool[(pool[site] == 1) & (~pool['SeriesInstanceUID'].isin(seen))]\n        if not cand.empty:\n            take = cand.head(need)\n            extra.append(take)\n            seen.update(take['SeriesInstanceUID'].tolist())\n            cov[site] += len(take)\n\n    if extra:\n        boosted = pd.concat([df_pos] + extra, ignore_index=True)\n        boosted = boosted.sort_values(['has_localizer', 'has_segmentation'], ascending=False)\n        boosted = boosted.head(len(df_pos)).copy()\n        return boosted.sample(frac=1.0, random_state=seed).reset_index(drop=True)\n    return df_pos\n\ndef build_balanced_subset(df_all: pd.DataFrame, frac: float, seed: int, min_per_site: int = 0):\n    df_all = df_all.copy()\n    pos = df_all[df_all[GLOBAL_COL] == 1].copy()\n    neg = df_all[df_all[GLOBAL_COL] == 0].copy()\n    n_total = int(math.floor(len(df_all) * frac))\n    n_per_class = min(len(pos), len(neg), n_total // 2)\n\n    pos_s = _sample_class_with_priority(pos, n_per_class, seed=seed)\n    neg_s = _sample_class_with_priority(neg, n_per_class, seed=seed)\n\n    if min_per_site > 0:\n        pos_s = _enforce_min_site_coverage(pos_s, min_per_site=min_per_site, seed=seed)\n\n    subset = pd.concat([pos_s, neg_s], ignore_index=True)\n    return subset.sample(frac=1.0, random_state=seed).reset_index(drop=True)\n\n# Build 50% subset\nout50 = build_balanced_subset(\n    train,\n    frac=CONFIG[\"fractions\"][1],\n    seed=CONFIG[\"seed\"],\n    min_per_site=CONFIG[\"min_per_site_50\"]\n)\nout50.to_csv(\"balanced_cta_mra_subset_50.csv\", index=False)\n\nprint(\"Subset created: balanced_cta_mra_subset_50.csv\")\nprint(f\"Total rows: {len(out50)}\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-03T05:32:45.497312Z","iopub.execute_input":"2025-12-03T05:32:45.497662Z","iopub.status.idle":"2025-12-03T05:32:45.59958Z","shell.execute_reply.started":"2025-12-03T05:32:45.497635Z","shell.execute_reply":"2025-12-03T05:32:45.598693Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# STEP 2: RSNA IA — 2.5D Preprocessing for \"Teacher\" Training\n# (Updated: central Z-slab + thinning + light in-plane auto-crop)\n# - Handles corrupt/empty DICOMs safely\n# - Handles single-frame AND multi-frame series\n# - Moderate Z-cropping (keeps central head/neck band)\n# - CTA: 3 HU windows as RGB (multi-window)\n# - MRA: z-score → clip → min-max, replicated to 3 channels\n# - Saves per-slice SOPInstanceUIDs for localizer-guided training\n# Produces:\n#   /kaggle/working/preproc_teacher/{SeriesUID}.npz\n#   /kaggle/working/preproc_teacher_index.csv\n# ============================================================\n\n# 1) Ensure DICOM + JPEG2000 decoding works on Kaggle\n!pip install \"numpy<2.0\" \"pandas<2.2.0\" \"scipy<1.13.0\" \\\n             \"pydicom\" \"pylibjpeg-libjpeg\" \"pylibjpeg-openjpeg\" \\\n             -U --force-reinstall --no-deps -q\n\nimport os, sys, json, warnings\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nfrom tqdm.auto import tqdm\n\nwarnings.filterwarnings(\"ignore\")\n\n# Try importing decoding libs (safe no-op if missing)\ntry:\n    import pylibjpeg  # noqa: F401\n    import pylibjpeg_openjpeg  # noqa: F401\nexcept Exception:\n    pass\n\n# -------------------- CONFIG --------------------\nSERIES_ROOT = Path(\"/kaggle/input/rsna-intracranial-aneurysm-detection/series\")\n\n# Point this to your subset file (e.g. your balanced 50% subset)\nSUBSET_CSV = Path(\"balanced_cta_mra_subset_50.csv\")  # change to _20 if needed\n\nOUT_DIR = Path(\"/kaggle/working/preproc_teacher\")\nOUT_DIR.mkdir(parents=True, exist_ok=True)\n\nINDEX_CSV_OUT = Path(\"/kaggle/working/preproc_teacher/preproc_teacher_index.csv\")\n\n# Geometry / size\nTARGET_HW = (224, 224)  # (H, W)\n\n# Z-cropping / thinning\nKEEP_Z_FRAC = 0.7          # if None, use adaptive fraction\nADAPTIVE_KEEP_Z = False\nMIN_KEEP_FRAC = 0.50        # never keep less than 50% of slices\nMAX_KEEP_FRAC = 0.80        # never keep more than 80% of slices\nTHINNING_STRIDE = 1         # keep every 2nd slice after mid-crop\nMAX_SLICES = 160            # cap total slices per series; set None to disable\n\n# In-plane auto-crop (light)\nAUTO_CROP = True\nAUTO_CROP_PAD = 8           # padding (pixels) around brain/body mask\n\n# CTA HU windows (3-channel RGB)\nCTA_WINDOWS_RGB = [\n    (40.0, 400.0),   # soft / angiographic-ish\n    (100.0, 700.0),  # vessel / contrast\n    (600.0, 2800.0)  # bone\n]\n\n# MRA normalization\nMRA_CLIP_STD = 5.0\n\n# -------------------- Helpers --------------------\n\ndef _safe_pixel_spacing(ds):\n    \"\"\"Return (dy, dx) spacing with fallbacks.\"\"\"\n    ps = getattr(ds, \"PixelSpacing\", None)\n    if ps is None:\n        return 1.0, 1.0\n    if isinstance(ps, (str, bytes)):\n        vals = [float(x) for x in str(ps).replace(\"\\\\\", \" \").split()]\n    else:\n        vals = [float(x) for x in ps]\n    if len(vals) == 1:\n        vals = [vals[0], vals[0]]\n    return float(vals[0]), float(vals[1])\n\ndef _safe_z_from_ds(ds):\n    \"\"\"\n    Try to extract a sortable z-position:\n    - Prefer ImagePositionPatient[2]\n    - Fallback to InstanceNumber\n    - As last resort, 0.0\n    \"\"\"\n    ipp = getattr(ds, \"ImagePositionPatient\", None)\n    if ipp is not None:\n        try:\n            if isinstance(ipp, (str, bytes)):\n                vals = [float(x) for x in str(ipp).replace(\"\\\\\", \" \").split()]\n            else:\n                vals = [float(v) for v in ipp]\n            if len(vals) >= 3:\n                return float(vals[2])\n        except Exception:\n            pass\n    # Fallback: InstanceNumber\n    inst = getattr(ds, \"InstanceNumber\", None)\n    try:\n        return float(inst)\n    except Exception:\n        return 0.0\n\ndef _load_singleframe_series(files):\n    \"\"\"\n    Load a single-frame DICOM series into a list of slices & SOP IDs.\n    Returns:\n        slices: list of np.ndarray [H,W] float32 (in HU or MR units)\n        sops:   list of str SOPInstanceUID\n        dz, dy, dx: spacing (approximate)\n    \"\"\"\n    slices = []\n    sops = []\n    zs = []\n    dy = dx = dz = 1.0\n\n    for f in files:\n        try:\n            ds = pydicom.dcmread(str(f))\n        except Exception:\n            continue\n\n        # Skip if no pixels\n        if not hasattr(ds, \"pixel_array\"):\n            continue\n\n        try:\n            px = ds.pixel_array\n        except Exception:\n            continue\n\n        if px is None or px.size == 0:\n            continue\n        # handle RGB or other weird dims by squeezing\n        if px.ndim > 2:\n            # If [H,W,3], convert to grayscale; if [1,H,W] etc, squeeze\n            if px.ndim == 3 and px.shape[-1] == 3:\n                r, g, b = px[..., 0], px[..., 1], px[..., 2]\n                px = (0.299 * r + 0.587 * g + 0.114 * b).astype(np.float32)\n            else:\n                px = np.squeeze(px)\n        if px.ndim != 2 or px.shape[0] == 0 or px.shape[1] == 0:\n            continue\n\n        try:\n            slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n            inter = float(getattr(ds, \"RescaleIntercept\", 0.0))\n        except Exception:\n            slope, inter = 1.0, 0.0\n\n        img = px.astype(np.float32) * slope + inter\n        z = _safe_z_from_ds(ds)\n        sop = str(getattr(ds, \"SOPInstanceUID\", f.name))\n\n        slices.append(img)\n        sops.append(sop)\n        zs.append(z)\n\n        # Spacings (rough)\n        dy, dx = _safe_pixel_spacing(ds)\n        try:\n            dz = float(getattr(ds, \"SliceThickness\", 1.0))\n        except Exception:\n            dz = 1.0\n\n    if len(slices) == 0:\n        raise ValueError(\"No valid single-frame slices loaded\")\n\n    zs = np.array(zs, dtype=float)\n    order = np.argsort(zs)\n    slices_sorted = [slices[i] for i in order]\n    sops_sorted = [sops[i] for i in order]\n    return slices_sorted, sops_sorted, (dz, dy, dx)\n\ndef _load_multiframe_series(file_path):\n    \"\"\"\n    Load a multi-frame DICOM file.\n    Returns:\n        slices: list of [H,W] float32 images (one per frame)\n        sops:   list of SOPInstanceUID (repeated per frame)\n        dz, dy, dx: spacing (approximate)\n    \"\"\"\n    ds = pydicom.dcmread(str(file_path))\n\n    if not hasattr(ds, \"pixel_array\"):\n        raise ValueError(\"Multi-frame: no pixel_array\")\n\n    try:\n        arr = ds.pixel_array\n    except Exception as e:\n        raise ValueError(f\"Multi-frame pixel_array error: {e}\")\n\n    if arr is None or arr.size == 0:\n        raise ValueError(\"Multi-frame: empty pixel_array\")\n\n    # shapes: [F,H,W] or [F,1,H,W] or [F,H,W,3]\n    if arr.ndim == 4:\n        # if [F,1,H,W] -> squeeze channel\n        if arr.shape[1] in (1,):\n            arr = arr[:, 0]\n        elif arr.shape[-1] == 3:\n            # convert to grayscale per frame\n            arr_list = []\n            for f in range(arr.shape[0]):\n                frame = arr[f]\n                if frame.ndim == 3 and frame.shape[-1] == 3:\n                    r, g, b = frame[..., 0], frame[..., 1], frame[..., 2]\n                    arr_list.append((0.299 * r + 0.587 * g + 0.114 * b).astype(np.float32))\n                else:\n                    arr_list.append(np.squeeze(frame).astype(np.float32))\n            arr = np.stack(arr_list, axis=0)\n        else:\n            arr = np.squeeze(arr)\n    elif arr.ndim == 3:\n        # [F,H,W] is fine\n        pass\n    else:\n        # Unexpected\n        arr = np.squeeze(arr)\n        if arr.ndim != 3:\n            raise ValueError(f\"Unexpected multi-frame shape: {arr.shape}\")\n\n    F, H, W = arr.shape\n\n    try:\n        slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n        inter = float(getattr(ds, \"RescaleIntercept\", 0.0))\n    except Exception:\n        slope, inter = 1.0, 0.0\n\n    vol = arr.astype(np.float32) * slope + inter\n\n    # For per-frame z ordering, try PerFrameFunctionalGroupsSequence\n    zs = None\n    pffg = getattr(ds, \"PerFrameFunctionalGroupsSequence\", None)\n    if pffg:\n        zvals = []\n        for f in pffg:\n            pos_seq = f.get(\"PlanePositionSequence\", None) or f.get(\"PlanePositionPatientSequence\", None)\n            if pos_seq:\n                ipp = pos_seq[0].get(\"ImagePositionPatient\", None)\n                if ipp is not None:\n                    try:\n                        zvals.append(float(ipp[2]))\n                    except Exception:\n                        zvals.append(0.0)\n                else:\n                    zvals.append(0.0)\n            else:\n                zvals.append(0.0)\n        if len(zvals) == F:\n            zs = np.array(zvals, dtype=float)\n\n    if zs is None:\n        # fallback: order as-is\n        zs = np.arange(F, dtype=float)\n\n    order = np.argsort(zs)\n    vol = vol[order]\n\n    try:\n        dy, dx = _safe_pixel_spacing(ds)\n    except Exception:\n        dy, dx = 1.0, 1.0\n    try:\n        dz = float(getattr(ds, \"SliceThickness\", 1.0))\n    except Exception:\n        dz = 1.0\n\n    sop = str(getattr(ds, \"SOPInstanceUID\", file_path.name))\n    sops = [sop for _ in range(vol.shape[0])]\n\n    slices = [vol[i] for i in range(vol.shape[0])]\n    return slices, sops, (dz, dy, dx)\n\ndef load_dicom_series(series_uid):\n    \"\"\"\n    Load a DICOM series (single- or multi-frame) into:\n        vol_hu: [Z,H,W] float32 (HU for CT, raw for MR)\n        sops:   [Z] list of SOPInstanceUID strings\n        spacing: (dz,dy,dx)\n    \"\"\"\n    sdir = SERIES_ROOT / series_uid\n    if not sdir.exists():\n        raise FileNotFoundError(f\"Series directory not found: {sdir}\")\n\n    files = sorted(list(sdir.glob(\"**/*.dcm\")))\n    if not files:\n        raise FileNotFoundError(f\"No DICOM files in {sdir}\")\n\n    # Peek first file header\n    ds0 = pydicom.dcmread(str(files[0]), stop_before_pixels=True)\n    is_multi = hasattr(ds0, \"NumberOfFrames\") or hasattr(ds0, \"PerFrameFunctionalGroupsSequence\")\n\n    if is_multi:\n        slices, sops, spacing = _load_multiframe_series(files[0])\n    else:\n        slices, sops, spacing = _load_singleframe_series(files)\n\n    if len(slices) == 0:\n        raise ValueError(f\"No valid slices loaded for series {series_uid}\")\n\n    vol = np.stack(slices, axis=0).astype(np.float32)  # [Z,H,W]\n    return vol, sops, spacing\n\n# --------- Normalization & geometry helpers ---------\n\ndef _hu_to01(x, lo, hi):\n    x = np.clip(x, lo, hi)\n    return (x - lo) / (hi - lo + 1e-6)\n\ndef normalize_cta_to_rgb(vol_hu):\n    \"\"\"\n    CTA: apply 3 HU windows → stack as RGB channels.\n    Input:\n        vol_hu: [Z,H,W] float32 (HU)\n    Output:\n        vol_rgb: [Z,3,H,W] float32 in [0,1]\n    \"\"\"\n    Z, H, W = vol_hu.shape\n    out = np.empty((Z, 3, H, W), dtype=np.float32)\n    for z in range(Z):\n        slice_hu = vol_hu[z]\n        chans = []\n        for (lo, hi) in CTA_WINDOWS_RGB:\n            x = np.clip(slice_hu, lo, hi)\n            x = (x - lo) / (hi - lo + 1e-6)\n            chans.append(x.astype(np.float32))\n        out[z] = np.stack(chans, axis=0)\n    return out\n\ndef normalize_mra_to_rgb(vol_raw):\n    \"\"\"\n    MRA: z-score → clip → min-max; replicate to 3 channels.\n    Input:\n        vol_raw: [Z,H,W] float32\n    Output:\n        vol_rgb: [Z,3,H,W] float32 in [0,1]\n    \"\"\"\n    v = vol_raw.astype(np.float32)\n    m = float(v.mean())\n    s = float(v.std() + 1e-6)\n    v = (v - m) / s\n    v = np.clip(v, -MRA_CLIP_STD, MRA_CLIP_STD)\n    v = (v - v.min()) / (v.max() - v.min() + 1e-6)\n    Z, H, W = v.shape\n    vol_rgb = np.repeat(v[:, None, :, :], 3, axis=1)  # [Z,3,H,W]\n    return vol_rgb\n\ndef compute_bbox_from_vol01(vol01, pad=AUTO_CROP_PAD):\n    \"\"\"\n    Compute a 2D bbox (y0,y1,x0,x1) from a normalized [Z,H,W] volume.\n    We OR across Z, then find non-zero area in XY.\n    \"\"\"\n    Z, H, W = vol01.shape\n    # Binary mask of \"tissue\" vs background\n    mask = (vol01 > 0.05).astype(np.uint8)  # simple, robust threshold\n    M = mask.any(axis=0).astype(np.uint8)   # [H,W]\n    ys, xs = np.where(M > 0)\n    if ys.size == 0:\n        return None\n    y0 = max(0, ys.min() - pad)\n    y1 = min(H, ys.max() + pad + 1)\n    x0 = max(0, xs.min() - pad)\n    x1 = min(W, xs.max() + pad + 1)\n    return (y0, y1, x0, x1)\n\ndef crop_with_bbox(vol, bbox):\n    if bbox is None:\n        return vol\n    y0, y1, x0, x1 = bbox\n    return vol[:, y0:y1, x0:x1]\n\ndef adaptive_keep_frac(Z):\n    \"\"\"\n    Adaptive fraction of slices to keep based on depth,\n    clamped into [MIN_KEEP_FRAC, MAX_KEEP_FRAC].\n    \"\"\"\n    if Z <= 0:\n        return 1.0\n    frac = 200.0 / float(max(Z, 1))  # more aggressive for thick volumes\n    frac = max(MIN_KEEP_FRAC, min(MAX_KEEP_FRAC, frac))\n    return float(frac)\n\ndef keep_middle_z(vol, sops, frac=None):\n    \"\"\"\n    Keep a central Z slab of the volume (and matching sops).\n    Returns vol_mid, (z_start,z_end), sops_mid\n    \"\"\"\n    Z = vol.shape[0]\n    if Z == 0:\n        raise ValueError(\"Empty volume in keep_middle_z\")\n\n    if frac is None:\n        frac = adaptive_keep_frac(Z) if ADAPTIVE_KEEP_Z else MAX_KEEP_FRAC\n    k = max(1, int(round(Z * frac)))\n    k = min(Z, k)\n    start = (Z - k) // 2\n    end = start + k\n    vol_mid = vol[start:end]\n    sops_mid = sops[start:end]\n    return vol_mid, (start, end), sops_mid\n\ndef thin_and_cap(vol, sops):\n    \"\"\"\n    Apply thinning stride and cap the number of slices; keep sops aligned.\n    \"\"\"\n    Z = vol.shape[0]\n    idx = np.arange(Z, dtype=int)\n\n    if THINNING_STRIDE and THINNING_STRIDE > 1:\n        idx = idx[::THINNING_STRIDE]\n\n    if MAX_SLICES is not None and len(idx) > MAX_SLICES:\n        idx = np.linspace(0, len(idx) - 1, MAX_SLICES).round().astype(int)\n\n    vol_out = vol[idx]\n    sops_out = [sops[i] for i in idx]\n    return vol_out, sops_out\n\ndef resize_volume(vol_3ch, target_hw):\n    \"\"\"\n    Resize [Z,3,H,W] to [Z,3,target_H,target_W] using cv2.INTER_AREA.\n    \"\"\"\n    Z, C, H, W = vol_3ch.shape\n    th, tw = target_hw\n    out = np.empty((Z, C, th, tw), dtype=np.float32)\n    for z in range(Z):\n        for c in range(C):\n            out[z, c] = cv2.resize(vol_3ch[z, c], (tw, th), interpolation=cv2.INTER_AREA)\n    return out\n\n# -------------------- Main Preprocessing --------------------\n\nif not SUBSET_CSV.exists():\n    raise FileNotFoundError(\n        f\"Subset CSV not found at {SUBSET_CSV}. \"\n        \"Make sure you have generated balanced_cta_mra_subset_XX.csv first.\"\n    )\n\nsubset = pd.read_csv(SUBSET_CSV)\nsubset[\"Modality\"] = subset[\"Modality\"].astype(str).str.upper().str.strip()\nsubset = subset[subset[\"Modality\"].isin([\"CTA\", \"MRA\"])].reset_index(drop=True)\n\nprint(f\"Preprocessing {len(subset)} CTA/MRA series from {SUBSET_CSV.name} ...\")\n\nindex_rows = []\n\nfor i, row in tqdm(subset.iterrows(), total=len(subset)):\n    uid = str(row[\"SeriesInstanceUID\"])\n    mod = row[\"Modality\"]\n\n    try:\n        # Load raw volume and SOPs\n        vol_raw, sops, spacing = load_dicom_series(uid)   # [Z,H,W], list[str], (dz,dy,dx)\n        Z0, H0, W0 = vol_raw.shape\n\n        # ---- In-plane normalization for bbox (not final input) ----\n        if mod == \"CTA\":\n            # use the softest CTA window for a robust tissue mask\n            lo_soft, hi_soft = CTA_WINDOWS_RGB[0]\n            base01 = _hu_to01(vol_raw, lo_soft, hi_soft)\n        else:\n            # quick MRA norm: z-score + clip + min-max\n            v = vol_raw.astype(np.float32)\n            m = float(v.mean())\n            s = float(v.std() + 1e-6)\n            v = (v - m) / s\n            v = np.clip(v, -MRA_CLIP_STD, MRA_CLIP_STD)\n            base01 = (v - v.min()) / (v.max() - v.min() + 1e-6)\n\n        # ---- In-plane auto-crop (XY) ----\n        if AUTO_CROP:\n            bbox = compute_bbox_from_vol01(base01, pad=AUTO_CROP_PAD)\n            if bbox is not None:\n                vol_raw = crop_with_bbox(vol_raw, bbox)\n            else:\n                bbox = (-1, -1, -1, -1)\n        else:\n            bbox = (-1, -1, -1, -1)\n\n        # ---- Z cropping: keep central band ----\n        vol_mid, z_bounds, sops_mid = keep_middle_z(vol_raw, sops, frac=KEEP_Z_FRAC)\n\n        # ---- Thinning + cap ----\n        vol_mid, sops_mid = thin_and_cap(vol_mid, sops_mid)\n\n        # Final Z/H/W after cropping + thinning\n        Z, H, W = vol_mid.shape\n\n        if Z == 0 or H == 0 or W == 0:\n            raise ValueError(\"Volume became empty after cropping / thinning\")\n\n        # ---- Final normalization to 3-channel [0,1] ----\n        if mod == \"CTA\":\n            vol_3ch = normalize_cta_to_rgb(vol_mid)       # [Z,3,H,W] in [0,1]\n        else:\n            vol_3ch = normalize_mra_to_rgb(vol_mid)       # [Z,3,H,W] in [0,1]\n\n        # ---- Resize ----\n        vol_resized = resize_volume(vol_3ch, TARGET_HW)   # [Z,3,Ht,Wt]\n\n        # Store as uint8 to save space\n        vol_uint8 = (np.clip(vol_resized, 0.0, 1.0) * 255.0).astype(np.uint8)\n\n        # Safety: ensure sops length = Z; if mismatch, fix by clipping or repeating\n        if len(sops_mid) != vol_uint8.shape[0]:\n            if len(sops_mid) == 1:\n                sops_mid = [sops_mid[0] for _ in range(vol_uint8.shape[0])]\n            else:\n                tmp = list(sops_mid)[:vol_uint8.shape[0]]\n                if len(tmp) < vol_uint8.shape[0]:\n                    while len(tmp) < vol_uint8.shape[0]:\n                        tmp.append(tmp[-1])\n                sops_mid = tmp\n\n        npz_path = OUT_DIR / f\"{uid}.npz\"\n        np.savez_compressed(\n            npz_path,\n            volume=vol_uint8,                          # [Z,3,H,W] uint8\n            sops=np.array(sops_mid, dtype=object),     # [Z] object/str\n            spacing=np.array(spacing, dtype=np.float32),  # (dz,dy,dx)\n            modality=np.array(mod),\n            series_uid=np.array(uid),\n            bbox=np.array(bbox, dtype=np.int32),       # (y0,y1,x0,x1) or (-1,...)\n            z_bounds=np.array(z_bounds, dtype=np.int32)  # (z_start,z_end) before thinning\n        )\n\n        index_rows.append({\n            \"SeriesInstanceUID\": uid,\n            \"Modality\": mod,\n            \"npz_path\": str(npz_path),\n            \"z_depth\": int(vol_uint8.shape[0]),\n            \"h\": int(vol_uint8.shape[2]),\n            \"w\": int(vol_uint8.shape[3]),\n            \"spacing_dz\": float(spacing[0]),\n            \"spacing_dy\": float(spacing[1]),\n            \"spacing_dx\": float(spacing[2]),\n            \"bbox_y0\": int(bbox[0]),\n            \"bbox_y1\": int(bbox[1]),\n            \"bbox_x0\": int(bbox[2]),\n            \"bbox_x1\": int(bbox[3]),\n            \"z_start\": int(z_bounds[0]),\n            \"z_end\": int(z_bounds[1]),\n            \"z_orig\": int(Z0),\n        })\n\n    except Exception as e:\n        print(f\"[WARN] {uid}: {type(e).__name__}: {e}\", file=sys.stderr)\n        continue\n\nindex_df = pd.DataFrame(index_rows)\nindex_df.to_csv(INDEX_CSV_OUT, index=False)\n\nprint(\"\\n=== DONE ===\")\nprint(f\"Series processed: {len(index_df)} / {len(subset)}\")\nprint(f\"NPZ dir: {OUT_DIR}\")\nprint(f\"Index CSV: {INDEX_CSV_OUT}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-03T06:44:04.892402Z","iopub.execute_input":"2025-12-03T06:44:04.892978Z","execution_failed":"2025-12-03T07:14:14.73Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}