{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceType":"competition","sourceId":99552,"databundleVersionId":13851420},{"sourceType":"datasetVersion","sourceId":15040734,"datasetId":9628522,"databundleVersionId":15919893},{"sourceType":"datasetVersion","sourceId":3610416,"datasetId":2126553,"databundleVersionId":3663963},{"sourceType":"datasetVersion","sourceId":15126861,"datasetId":9686495,"databundleVersionId":16015164},{"sourceType":"datasetVersion","sourceId":15126867,"datasetId":9686499,"databundleVersionId":16015170},{"sourceType":"datasetVersion","sourceId":15020557,"datasetId":9615023,"databundleVersionId":15897969},{"sourceType":"datasetVersion","sourceId":14998015,"datasetId":9600370,"databundleVersionId":15872863},{"sourceType":"modelInstanceVersion","sourceId":612683,"databundleVersionId":14140664,"modelInstanceId":460275,"modelId":476073}],"dockerImageVersionId":31286,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Physics-Grounded Synthetic Supervision for Controlled CT Restoration\n\n---\n\n## Research Map — This notebook argues for a **method**, supported by completed evidence blocks\n\n### Core problem\nIn medical imaging, perfectly paired training data are rare. We usually **do not** have the same patient scanned once in a pristine state and once again under a realistic degradation process with identical anatomy.\n\nSo the scientific question here is not simply:\n\n> “Can I build a stronger restoration network?”\n\nIt is:\n\n> **When true paired data are unavailable, can physics-grounded degradation synthesis provide supervision that is reliable enough to train a clinically relevant and conservative CT restoration model?**\n\n---\n\n## What is the main contribution here?\nThis notebook uses **V-Ultimate** as the implementation vehicle, but the broader contribution is the **methodology**:\n\n1. **Generate paired supervision from physics-grounded degradations** rather than arbitrary image corruption.\n2. **Train a restoration policy under controlled, do-no-harm constraints** so edits remain conservative rather than aggressive.\n3. **Judge the synthetic supervision through transfer, clinical-proxy behavior, and anatomy-aware evaluation**, not by visual sharpness alone.\n\nSo the scientific claim is **not** that this exact architecture is the final or universally best restoration system.\n\nThe claim is narrower and more useful:\n\n> **Physics-grounded synthetic supervision can be a credible way to train CT restoration when real paired data are missing, provided that the learned behavior is clinically useful, externally transferable, and anatomically structured.**\n\n---\n\n## What this notebook actually shows in this run\nThe completed results support a **three-part methodology story**:\n\n- On a **strict held-out RSNA OOD CT/CTA set (N = 50)**, V-Ultimate achieved **mean target gain = 0.0324**, **target win rate = 66%**, **iatrogenic rate = 32%**, and the **highest PSNR = 42.74 dB** among the compared methods.\n- On **external paired Mayo low-dose CT (10 unique matched series)**, V-Ultimate improved over Quarter Dose by **+2.35 dB PSNR**, **+0.00203 MAE reduction**, **+0.04099 SSIM**, and **100% PSNR / MAE / SSIM win rate**, with a bootstrap **95% CI of +2.03 to +2.67 dB** for PSNR.\n- In the same Mayo transfer test, V-Ultimate remained **more conservative than Gaussian smoothing**, with **mean absolute change vs Quarter = 0.00316** and **max absolute change = 0.02463**, compared with Gaussian’s **0.00690** mean and **0.15326** max.\n- In **TotalSegmentator analysis on the same 10 Mayo series**, the original threshold of **0.05** revealed essentially no changed voxels, showing that the model’s edits are very small. After lowering the analysis threshold to **0.01**, the edits became measurable and structured: mean Dice vs Full improved from **0.6137 (Quarter)** to **0.7395 (V-Ultimate)** over **238 masks**, with changes concentrating in specific anatomical structures rather than uniformly across the volume.\n\nThese results do **not** say that V-Ultimate wins every metric against every baseline.  \nInstead, they show that the synthetic-supervision pipeline produces a model with a strong **balance of usefulness, restraint, external transfer, and anatomical structure**, rather than a narrow one-metric win.\n\n---\n\n## Evidence chain used in this notebook\n\n| Evidence block | What this run shows | Why it matters for the methodology |\n|---|---|---|\n| **Clinical Rescue Matrix (strict RSNA OOD, N = 50)** | V-Ultimate achieves the **highest overall mean target gain (0.0324)**, the **highest target win rate (66%)**, and the **highest PSNR (42.74 dB)**, while keeping iatrogenic rate tied for the lowest among the stronger methods | Shows that the synthetic supervision is **internally useful**, not just visually plausible |\n| **Mayo external paired transfer (10 unique matched series)** | V-Ultimate achieves **+2.35 dB PSNR**, **+0.04099 SSIM**, and **100% PSNR win rate** over Quarter Dose, with consistently positive bootstrap confidence intervals | This is the strongest **sim-to-real transfer** evidence for the synthetic supervision claim |\n| **TotalSegmentator anatomy analysis (10 Mayo series)** | At a coarse threshold the edits are too small to register; at a finer threshold they become visible and are **anatomically structured**, while Dice improves from **0.6137 to 0.7395** | Shows that the restoration is **conservative but not random**: edits are small, structured, and anatomically meaningful |\n\n---\n\n## What this notebook does **not** claim\nTo stay scientifically honest, this notebook does **not** prove that synthetic low-quality CT is distributionally identical to every real scanner protocol.\n\nIt also does **not** claim that V-Ultimate is the best possible filter on every metric or every task.\n\nInstead, it argues something narrower and stronger:\n\n> If a model trained on physics-grounded synthetic pairs repeatedly shows clinically relevant behavior on strict OOD data, transfers to real external low-dose CT, and makes small but anatomically structured edits, then the synthetic supervision is **useful and empirically supported as a training signal**, even if it is not a perfect replica of reality.\n\n---\n\n## Reading guide\n- If you care about **how the method is trained**, start with the training and configuration sections.\n- If you care about **why the synthetic supervision should be believed**, focus on the completed evidence blocks: **CRM, Mayo, and TotalSegmentator**.\n- If you care about **stability**, read the Monte Carlo section as a **supplementary robustness probe** rather than a headline result in this run.\n- If you care about **Science Fair framing**, read the final “Judge Summary” and abstract-style closing cells.","metadata":{}},{"cell_type":"markdown","source":"# ═══════════════════════════════════════════════════════════════\n# Part I: Method Setup — Synthetic supervision + controlled restoration\n# ═══════════════════════════════════════════════════════════════\n\nThis section keeps the **original training code and model implementation intact**.  \nThe purpose of the rewritten narrative is to make the scientific framing clearer:\n\n- **Primary contribution:** a way to generate and validate training supervision when real pairs are absent\n- **Implementation vehicle:** V-Ultimate (the restoration model used to test that supervision)\n- **Evaluation principle:** safety and transfer matter more than “sharpest looking image”\n\nIn other words, the model is important, but it serves the larger question:\n\n> **Can physics-grounded synthetic pairs support a restoration system that generalizes beyond the simulator?**\n","metadata":{}},{"cell_type":"markdown","source":"# §2 Why the architecture is present, but not the whole claim\n\n## What to tell a judge\nA fair question is:\n> “Is your contribution the network, or the data-generation methodology?”\n\nThe most accurate answer is:\n\n> **The network is the experimental vehicle; the broader contribution is the methodology for creating supervision and validating it responsibly.**\n\nThis section therefore explains the model only as much as needed to understand the experiments:\n\n- the model must be able to ingest physically degraded inputs,\n- it must be able to restore conservatively,\n- and it must expose behavior that can be stress-tested downstream.\n\n## Why keep the model fixed here?\nBecause this notebook is **not** trying to win by endless architecture tweaking.  \nKeeping the training model unchanged helps preserve the main story:\n\n1. build synthetic supervision from physics,\n2. train a controlled restorer on it,\n3. ask whether the learned behavior survives external scrutiny.\n\nSo if judges are short on time, this cell can be treated as **implementation background**, not the headline contribution.\n","metadata":{}},{"cell_type":"code","source":"# =====================================================================\n# V-Ultimate FULL TRAINING CELL \n# =====================================================================\n\nimport os, gc, math, time, random\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport pydicom\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader, Sampler\nfrom contextlib import nullcontext\nimport torch.fft\n\n# -----------------------------\n# Environment\n# -----------------------------\ntry:\n    cv2.setNumThreads(0)\nexcept Exception:\n    pass\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(\"Device:\", device)\n\nUSE_AMP = (device.type == \"cuda\")\nAMP_CTX = lambda: torch.amp.autocast(\"cuda\") if USE_AMP else nullcontext()\nscaler = torch.amp.GradScaler(\"cuda\") if USE_AMP else None\n\ntorch.backends.cudnn.benchmark = True\nif device.type == \"cuda\":\n    try:\n        torch.backends.cuda.matmul.allow_tf32 = True\n        torch.backends.cudnn.allow_tf32 = True\n        torch.set_float32_matmul_precision(\"high\")\n    except Exception:\n        pass\n\n# ============================================================\n# Configuration\n# ============================================================\n\n# --- Data paths ---\nRSNA_DATA_ROOT = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/series\"\nTRAIN_LOCALIZERS_CSV = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train_localizers.csv\"\nMETA_CSV = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train.csv\"\n\nOUT_TRAIN_UIDS = \"/kaggle/working/train_uids_ultimate.csv\"\nOUT_VAL_UIDS   = \"/kaggle/working/val_uids_ultimate.csv\"\n\nSAVE_BEST = \"/kaggle/working/deblur_ultimate_film_auth_best.pt\"\nSAVE_LAST = \"/kaggle/working/deblur_ultimate_film_auth_last.pt\"\n\n# --- Training scale ---\nSEED = 2026\nN_TRAIN_UIDS = 300\nN_VAL_UIDS   = 20\nEPOCHS = 24\nBATCH_SIZE = 32\nBATCHES_PER_EPOCH = 400\nNUM_WORKERS = 4\n\n# --- Optimiser ---\nLR = 2e-4\nWEIGHT_DECAY = 1e-4\nGRAD_CLIP = 1.0\n\n# --- Image parameters ---\nTARGET_D, TARGET_H, TARGET_W = 64, 448, 448\nPATCH_SIZE = 112\nPATCHES_PER_SLICE = 1\n\nHU_MIN, HU_MAX = -1024.0, 3072.0\nHU_RANGE = HU_MAX - HU_MIN\n\n# --- Degradation parameters ---\nDIFFUSION_ALPHA = 0.20\nBLUR_LEVELS = [0, 1, 3, 5, 8]\nBLUR_LEVEL_MAX = float(max(BLUR_LEVELS))\nBLUR_T_MAX = BLUR_LEVEL_MAX\n\nP_IDENTITY = 0.20\nENABLE_MOTION = True\nP_MOTION = 0.15\n\nPEAK_RANGE_QUARTER, SIGMA_E_QUARTER = (3000.0, 6000.0), (0.01, 0.02)\nPEAK_RANGE_EXTREME, SIGMA_E_EXTREME = (1000.0, 3000.0), (0.02, 0.04)\n\nREGIME_PROBS = {\"clean\": 0.30, \"typical\": 0.45, \"hard\": 0.15, \"motion\": 0.10}\n\n# --- Anatomy-aware patch sampling ---\nANATOMY_REJECT_TRIES = 15\nANATOMY_MEAN_TH = 0.05\nANATOMY_STD_TH  = 0.02\nP_RANDOM_PATCH  = 0.10\n\n# --- Loss weights ---\nW_CHARBONNIER = 1.0\nW_SSIM = 0.20\nW_SOBEL = 0.10\nW_LAP = 0.05\nW_FFT = 0.05\nFFT_FCUTOFF = 0.20\nFFT_ONLY_IF_T_LE = 10.0\n\nW_CHANGE_ID = 10.0\nW_LOW_T_EDIT = 2.0\nLOW_T_THR = 0.20\nW_AUTH_TV = 1e-3\n\nRES_MIN, RES_MAX = 0.02, 0.15\n\n# ============================================================\n# Random seed\n# ============================================================\n\ndef seed_all(seed):\n    random.seed(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed_all(seed)\n\nseed_all(SEED)\n\ndef _normalize_probs(d):\n    s = float(sum(d.values()))\n    return {k: float(v) / s for k, v in d.items()}\n\nREGIME_PROBS = _normalize_probs(REGIME_PROBS)\n\n# ============================================================\n# UID list construction\n# ============================================================\n\ndef build_ct_only_uid_lists(meta_csv, rsna_series_root, localizers_csv, n_train, n_val, seed):\n    rsna_uids = set(\n        [u for u in os.listdir(rsna_series_root)\n         if os.path.isdir(os.path.join(rsna_series_root, u)) and not u.startswith(\".\")]\n    )\n\n    localizer_uids = set()\n    if localizers_csv and os.path.exists(localizers_csv):\n        df_loc = pd.read_csv(localizers_csv)\n        localizer_uids = set(df_loc[df_loc.columns[0]].astype(str).tolist())\n\n    meta = pd.read_csv(meta_csv)\n    meta_ct = meta[meta[\"Modality\"].isin({\"CT\", \"CTA\"})]\n\n    ct_candidates = [\n        u for u in set(meta_ct[\"SeriesInstanceUID\"].astype(str))\n        if u in rsna_uids and u not in localizer_uids\n    ]\n\n    rng = random.Random(seed)\n    rng.shuffle(ct_candidates)\n\n    train_uids = ct_candidates[:min(n_train, len(ct_candidates))]\n    rem_ct = [u for u in ct_candidates if u not in set(train_uids)]\n    rng.shuffle(rem_ct)\n    val_uids = rem_ct[:min(n_val, len(rem_ct))]\n    return train_uids, val_uids\n\n# ============================================================\n# DICOM loading\n# ============================================================\n\ndef get_sorted_dicom_files(series_path):\n    files = [f for f in os.listdir(series_path) if not f.startswith(\".\")]\n    pairs, ok = [], True\n    for f in files:\n        try:\n            ds = pydicom.dcmread(os.path.join(series_path, f), stop_before_pixels=True, force=True)\n            if getattr(ds, \"InstanceNumber\", None) is None:\n                ok = False\n                break\n            pairs.append((int(ds.InstanceNumber), os.path.join(series_path, f)))\n        except Exception:\n            ok = False\n            break\n    if ok and len(pairs) == len(files):\n        return [p[1] for p in sorted(pairs, key=lambda x: x[0])]\n    return [os.path.join(series_path, f) for f in sorted(files)]\n\ndef load_series_volume(uid, series_root, target_shape=(64, 448, 448)):\n    series_path = os.path.join(series_root, uid)\n    if not os.path.isdir(series_path):\n        return None\n\n    dcm_files = get_sorted_dicom_files(series_path)\n    tD, tH, tW = target_shape\n    if len(dcm_files) < 10:\n        return None\n\n    if len(dcm_files) != tD:\n        idxs = np.linspace(0, len(dcm_files) - 1, tD).astype(int)\n        dcm_files = [dcm_files[i] for i in idxs]\n\n    slices = []\n    for fp in dcm_files:\n        try:\n            ds = pydicom.dcmread(fp, force=True)\n            hu = ds.pixel_array.astype(np.float32) * float(getattr(ds, \"RescaleSlope\", 1.0)) \\\n                 + float(getattr(ds, \"RescaleIntercept\", 0.0))\n            x = (np.clip(hu, HU_MIN, HU_MAX) - HU_MIN) / HU_RANGE\n            x = cv2.resize(x, (tW, tH), interpolation=cv2.INTER_LINEAR)\n            slices.append(x.astype(np.float32))\n        except Exception:\n            continue\n\n    if len(slices) < int(0.8 * tD):\n        return None\n\n    while len(slices) < tD:\n        slices.append(slices[-1].copy())\n\n    return np.stack(slices[:tD], axis=0).astype(np.float32)\n\nclass VolumeLRU:\n    def __init__(self, max_items=12):\n        self.max_items = int(max_items)\n        self.cache = {}\n        self.order = []\n\n    def get(self, key):\n        if key not in self.cache:\n            return None\n        self.order.remove(key)\n        self.order.append(key)\n        return self.cache[key]\n\n    def put(self, key, value):\n        if key in self.cache:\n            self.order.remove(key)\n        self.cache[key] = value\n        self.order.append(key)\n        if len(self.order) > self.max_items:\n            old = self.order.pop(0)\n            self.cache.pop(old, None)\n\n# ============================================================\n# Degradation engine\n# ============================================================\n\ndef gaussian_psf_surrogate(img01, blur_level, alpha=0.20):\n    if blur_level <= 0:\n        return img01\n    sigma = math.sqrt(max(1e-8, 2.0 * alpha * float(blur_level)))\n    return np.clip(\n        cv2.GaussianBlur(img01, (0, 0), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REPLICATE),\n        0.0, 1.0\n    )\n\ndef motion_artifact_surrogate(img01, length=None, angle=None):\n    length = length or random.choice([3, 5, 7, 9, 11])\n    if length <= 1:\n        return img01\n    angle = angle or random.uniform(0, 180)\n\n    k = np.zeros((length, length), dtype=np.float32)\n    c = length // 2\n    cos_a, sin_a = np.cos(np.radians(angle)), np.sin(np.radians(angle))\n    for i in range(length):\n        x = int(c + (i - c) * cos_a)\n        y = int(c + (i - c) * sin_a)\n        if 0 <= x < length and 0 <= y < length:\n            k[y, x] = 1.0\n    if k.sum() > 0:\n        k /= k.sum()\n\n    return np.clip(cv2.filter2D(img01, -1, k, borderType=cv2.BORDER_REPLICATE), 0.0, 1.0)\n\ndef mixed_poisson_gaussian(img01, mode=\"quarter\"):\n    if mode == \"clean\":\n        return img01\n    peak_rng, sigma_rng = (\n        (PEAK_RANGE_EXTREME, SIGMA_E_EXTREME) if mode == \"extreme\"\n        else (PEAK_RANGE_QUARTER, SIGMA_E_QUARTER)\n    )\n    peak = random.uniform(*peak_rng)\n    sigma_e = random.uniform(*sigma_rng)\n    noisy_p = np.random.poisson(np.clip(img01 * peak, 0, None)).astype(np.float32) / peak\n    noisy_g = np.random.randn(*img01.shape).astype(np.float32) * sigma_e\n    return np.clip(noisy_p + noisy_g, 0.0, 1.0)\n\ndef _choice_weighted(prob_dict):\n    r = random.random()\n    acc = 0.0\n    for k, p in prob_dict.items():\n        acc += p\n        if r <= acc:\n            return k\n    return list(prob_dict.keys())[-1]\n\ndef _sample_regime_params():\n    reg = _choice_weighted(REGIME_PROBS)\n    if reg == \"clean\":\n        return 0, \"clean\", False\n    if reg == \"typical\":\n        return random.choice([1, 3, 5]), \"quarter\", False\n    if reg == \"hard\":\n        return random.choice([3, 5, 8]), \"extreme\", False\n    if reg == \"motion\":\n        return random.choice([1, 3, 5]), random.choice([\"quarter\", \"extreme\"]), True\n    return 3, \"quarter\", False\n\n# ============================================================\n# Dataset\n# ============================================================\n\nclass CTDeblur25D(Dataset):\n    def __init__(self, uids, series_root, target_shape=(64, 448, 448), patch_size=112, patches_per_slice=1):\n        self.uids = list(uids)\n        self.series_root = series_root\n        self.target_shape = target_shape\n        self.patch_size = int(patch_size)\n        self.patches_per_slice = int(patches_per_slice)\n        self.cache = VolumeLRU(max_items=8)\n        self.items = [(ui, z) for ui in range(len(self.uids)) for z in range(1, target_shape[0] - 1)]\n\n    def __len__(self):\n        return len(self.items)\n\n    def _sample_patch_xy(self, cent, ps):\n        H, W = cent.shape\n        if random.random() < P_RANDOM_PATCH:\n            return np.random.randint(0, H - ps + 1), np.random.randint(0, W - ps + 1)\n\n        for _ in range(ANATOMY_REJECT_TRIES):\n            y = np.random.randint(0, H - ps + 1)\n            x = np.random.randint(0, W - ps + 1)\n            patch = cent[y:y+ps, x:x+ps]\n            if patch.mean() > ANATOMY_MEAN_TH and patch.std() > ANATOMY_STD_TH:\n                return y, x\n\n        return (H - ps) // 2, (W - ps) // 2\n\n    def __getitem__(self, idx):\n        ui, z = self.items[idx]\n        uid = self.uids[ui]\n\n        vol = self.cache.get(uid)\n        if vol is None:\n            vol = load_series_volume(uid, self.series_root, self.target_shape)\n            if vol is None:\n                return self.__getitem__(random.randint(0, len(self.items) - 1))\n            self.cache.put(uid, vol)\n\n        ps = self.patch_size\n        y, x = self._sample_patch_xy(vol[z], ps)\n\n        clean = vol[z][y:y+ps, x:x+ps].copy()\n        pp = vol[z-1][y:y+ps, x:x+ps].copy()\n        cc = vol[z][y:y+ps, x:x+ps].copy()\n        nn_ = vol[z+1][y:y+ps, x:x+ps].copy()\n\n        is_identity = 0.0\n        if random.random() < P_IDENTITY:\n            blur_level, dose_mode, do_motion, is_identity = 0, \"clean\", False, 1.0\n        else:\n            blur_level, dose_mode, do_motion = _sample_regime_params()\n            do_motion = do_motion and ENABLE_MOTION and (random.random() < P_MOTION)\n\n        bp = gaussian_psf_surrogate(pp, blur_level, alpha=DIFFUSION_ALPHA)\n        bc = gaussian_psf_surrogate(cc, blur_level, alpha=DIFFUSION_ALPHA)\n        bn = gaussian_psf_surrogate(nn_, blur_level, alpha=DIFFUSION_ALPHA)\n\n        if do_motion:\n            L = random.choice([3, 5, 7, 9, 11])\n            A = random.uniform(0, 180)\n            bp = motion_artifact_surrogate(bp, L, A)\n            bc = motion_artifact_surrogate(bc, L, A)\n            bn = motion_artifact_surrogate(bn, L, A)\n\n        if dose_mode != \"clean\":\n            bp = mixed_poisson_gaussian(bp, dose_mode)\n            bc = mixed_poisson_gaussian(bc, dose_mode)\n            bn = mixed_poisson_gaussian(bn, dose_mode)\n\n        cp = clean.copy()\n        if random.random() > 0.5:\n            cp, bp, bc, bn = cp[::-1].copy(), bp[::-1].copy(), bc[::-1].copy(), bn[::-1].copy()\n        if random.random() > 0.5:\n            cp, bp, bc, bn = cp[:, ::-1].copy(), bp[:, ::-1].copy(), bc[:, ::-1].copy(), bn[:, ::-1].copy()\n        k = random.randint(0, 3)\n        if k > 0:\n            cp, bp, bc, bn = np.rot90(cp, k).copy(), np.rot90(bp, k).copy(), np.rot90(bc, k).copy(), np.rot90(bn, k).copy()\n\n        t_norm = float(blur_level) / BLUR_LEVEL_MAX if BLUR_LEVEL_MAX > 0 else 0.0\n\n        inp = np.stack([bp, bc, bn, np.full_like(bc, t_norm, dtype=np.float32)], axis=0).astype(np.float32)\n        tgt = cp[np.newaxis, ...].astype(np.float32)\n\n        meta = np.array([\n            t_norm,\n            1.0 if do_motion else 0.0,\n            1.0 if dose_mode == \"clean\" else 0.0,\n            is_identity\n        ], dtype=np.float32)\n\n        return (\n            torch.from_numpy(inp).float(),\n            torch.from_numpy(tgt).float(),\n            torch.from_numpy(meta).float()\n        )\n\nclass UIDBatchSampler(Sampler):\n    def __init__(self, dataset, batch_size, seed=42, batches_per_epoch=None):\n        self.dataset = dataset\n        self.batch_size = int(batch_size)\n        self.rng = random.Random(seed)\n\n        self.by_ui = {}\n        for idx, (ui, z) in enumerate(dataset.items):\n            self.by_ui.setdefault(ui, []).append(idx)\n\n        self.ui_keys = list(self.by_ui.keys())\n        self.batches_per_epoch = int(batches_per_epoch) if batches_per_epoch else len(dataset) // self.batch_size\n\n    def __len__(self):\n        return self.batches_per_epoch\n\n    def __iter__(self):\n        for _ in range(self.batches_per_epoch):\n            ui = self.rng.choice(self.ui_keys)\n            pool = self.by_ui[ui]\n            if len(pool) >= self.batch_size:\n                yield self.rng.sample(pool, self.batch_size)\n            else:\n                yield [self.rng.choice(pool) for _ in range(self.batch_size)]\n\n# ============================================================\n# New model: FiLM + authority map\n# ============================================================\n\nclass SEBlock(nn.Module):\n    def __init__(self, c, r=4):\n        super().__init__()\n        self.fc = nn.Sequential(\n            nn.AdaptiveAvgPool2d(1),\n            nn.Conv2d(c, max(1, c // r), 1, bias=False),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(max(1, c // r), c, 1, bias=False),\n            nn.Sigmoid(),\n        )\n\n    def forward(self, x):\n        return x * self.fc(x)\n\nclass ResBlockPhysics(nn.Module):\n    def __init__(self, ic, oc):\n        super().__init__()\n        self.conv = nn.Sequential(\n            nn.Conv2d(ic, oc, 3, padding=1, bias=True),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(oc, oc, 3, padding=1, bias=True),\n        )\n        self.shortcut = nn.Conv2d(ic, oc, 1, bias=True) if ic != oc else nn.Identity()\n\n    def forward(self, x):\n        return F.relu(self.conv(x) + self.shortcut(x), inplace=True)\n\nclass UpsamplePhysicsUltimate(nn.Module):\n    def __init__(self, ic, oc):\n        super().__init__()\n        self.up = nn.Sequential(\n            nn.Upsample(scale_factor=2.0, mode=\"nearest\"),\n            nn.Conv2d(ic, oc, 3, padding=1, bias=True),\n            nn.ReLU(inplace=True),\n        )\n\n    def forward(self, x):\n        return self.up(x)\n\nclass FiLM2d(nn.Module):\n    def __init__(self, channels, meta_dim=4, hidden=64):\n        super().__init__()\n        self.net = nn.Sequential(\n            nn.Linear(meta_dim, hidden),\n            nn.ReLU(inplace=True),\n            nn.Linear(hidden, channels * 2),\n        )\n        self.channels = channels\n\n    def forward(self, x, meta):\n        gb = self.net(meta)\n        gamma, beta = torch.chunk(gb, 2, dim=1)\n        gamma = gamma.view(-1, self.channels, 1, 1)\n        beta = beta.view(-1, self.channels, 1, 1)\n        return x * (1.0 + gamma) + beta\n\nclass DeblurUNet25D_Ultimate(nn.Module):\n    def __init__(\n        self,\n        in_ch=4,\n        out_ch=1,\n        base=32,\n        res_min=0.02,\n        res_max=0.15,\n        meta_dim=4,\n        film_hidden=64,\n        authority_bias_init=2.0,\n    ):\n        super().__init__()\n        self.res_min = float(res_min)\n        self.res_max = float(res_max)\n        self.meta_dim = int(meta_dim)\n\n        c = [base, base * 2, base * 4, base * 8]\n\n        self.se = SEBlock(in_ch)\n\n        self.enc1 = ResBlockPhysics(in_ch, c[0])\n        self.enc2 = ResBlockPhysics(c[0], c[1])\n        self.enc3 = ResBlockPhysics(c[1], c[2])\n        self.enc4 = ResBlockPhysics(c[2], c[3])\n        self.pool = nn.MaxPool2d(2)\n\n        self.up3 = UpsamplePhysicsUltimate(c[3], c[2])\n        self.dec3 = ResBlockPhysics(c[2] * 2, c[2])\n\n        self.up2 = UpsamplePhysicsUltimate(c[2], c[1])\n        self.dec2 = ResBlockPhysics(c[1] * 2, c[1])\n\n        self.up1 = UpsamplePhysicsUltimate(c[1], c[0])\n        self.dec1 = ResBlockPhysics(c[0] * 2, c[0])\n\n        self.film_e1 = FiLM2d(c[0], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_e2 = FiLM2d(c[1], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_e3 = FiLM2d(c[2], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_e4 = FiLM2d(c[3], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_d3 = FiLM2d(c[2], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_d2 = FiLM2d(c[1], meta_dim=meta_dim, hidden=film_hidden)\n        self.film_d1 = FiLM2d(c[0], meta_dim=meta_dim, hidden=film_hidden)\n\n        self.out_conv = nn.Conv2d(c[0], out_ch, 1, bias=True)\n\n        self.auth_head = nn.Sequential(\n            nn.Conv2d(c[0], c[0] // 2, 3, padding=1, bias=True),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(c[0] // 2, 1, 1, bias=True),\n        )\n        nn.init.constant_(self.auth_head[-1].bias, float(authority_bias_init))\n\n    def forward(self, x, meta=None, return_aux=False):\n        bc = x[:, 1:2]\n        tch = x[:, 3:4]\n\n        if meta is None:\n            B = x.shape[0]\n            t_scalar = torch.mean(tch, dim=(2, 3)).view(B, 1)\n            zeros = torch.zeros(B, 3, device=x.device, dtype=x.dtype)\n            meta = torch.cat([t_scalar, zeros], dim=1)\n\n        e1 = self.enc1(self.se(x))\n        e1 = self.film_e1(e1, meta)\n\n        e2 = self.enc2(self.pool(e1))\n        e2 = self.film_e2(e2, meta)\n\n        e3 = self.enc3(self.pool(e2))\n        e3 = self.film_e3(e3, meta)\n\n        e4 = self.enc4(self.pool(e3))\n        e4 = self.film_e4(e4, meta)\n\n        d3 = self.dec3(torch.cat([self.up3(e4), e3], dim=1))\n        d3 = self.film_d3(d3, meta)\n\n        d2 = self.dec2(torch.cat([self.up2(d3), e2], dim=1))\n        d2 = self.film_d2(d2, meta)\n\n        d1 = self.dec1(torch.cat([self.up1(d2), e1], dim=1))\n        d1 = self.film_d1(d1, meta)\n\n        residual = torch.tanh(self.out_conv(d1))\n\n        base_rmax = self.res_min + (self.res_max - self.res_min) * tch\n        authority = torch.sigmoid(self.auth_head(d1))\n        rmax_map = base_rmax * authority\n\n        pred_soft = (bc + residual * rmax_map).clamp(0.0, 1.0)\n        pred = torch.where(tch <= 1e-8, bc, pred_soft)\n\n        if return_aux:\n            return pred, {\n                \"authority\": authority,\n                \"rmax_map\": rmax_map,\n                \"residual\": residual,\n            }\n        return pred\n\n# ============================================================\n# Loss functions\n# ============================================================\n\ndef charbonnier_loss(pred, target, eps=1e-3):\n    return torch.mean(torch.sqrt((pred - target) ** 2 + eps ** 2))\n\ndef fft_spectrum_loss(pred, target, fft_mask):\n    pred_fft = torch.fft.rfft2(pred.float(), dim=(-2, -1), norm=\"ortho\")\n    tgt_fft = torch.fft.rfft2(target.float(), dim=(-2, -1), norm=\"ortho\")\n    m = fft_mask.view(1, 1, *fft_mask.shape)\n    return charbonnier_loss(torch.abs(pred_fft) * m, torch.abs(tgt_fft) * m)\n\ndef ssim_loss(pred, target, window_size=11):\n    C1, C2 = 0.01 ** 2, 0.03 ** 2\n    pad = window_size // 2\n\n    mu_x = F.avg_pool2d(pred, window_size, stride=1, padding=pad)\n    mu_y = F.avg_pool2d(target, window_size, stride=1, padding=pad)\n\n    sigma_x2 = F.avg_pool2d(pred ** 2, window_size, stride=1, padding=pad) - mu_x ** 2\n    sigma_y2 = F.avg_pool2d(target ** 2, window_size, stride=1, padding=pad) - mu_y ** 2\n    sigma_xy = F.avg_pool2d(pred * target, window_size, stride=1, padding=pad) - mu_x * mu_y\n\n    ssim_map = ((2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)) / (\n        (mu_x ** 2 + mu_y ** 2 + C1) * (sigma_x2 + sigma_y2 + C2) + 1e-8\n    )\n    return 1.0 - ssim_map.mean()\n\ndef _make_fft_mask(H, W, fcut=0.20, device=\"cpu\"):\n    fy = torch.fft.fftfreq(H, d=1.0, device=device).view(H, 1).abs()\n    fx = torch.fft.rfftfreq(W, d=1.0, device=device).view(1, W // 2 + 1).abs()\n    return (torch.sqrt(fx * fx + fy * fy) >= fcut).float()\n\nclass UltimatePhysicsLoss(nn.Module):\n    def __init__(self, patch_size=112):\n        super().__init__()\n        self.sobel_x = torch.tensor([[[-1., 0., 1.], [-2., 0., 2.], [-1., 0., 1.]]]).view(1, 1, 3, 3).to(device)\n        self.sobel_y = torch.tensor([[[-1., -2., -1.], [0., 0., 0.], [1., 2., 1.]]]).view(1, 1, 3, 3).to(device)\n        self.lap = torch.tensor([[[0., 1., 0.], [1., -4., 1.], [0., 1., 0.]]]).view(1, 1, 3, 3).to(device)\n        self.register_buffer(\"fft_mask\", _make_fft_mask(patch_size, patch_size, fcut=FFT_FCUTOFF, device=device))\n\n    def forward(self, pred, target, allow_fft=False):\n        total = W_CHARBONNIER * charbonnier_loss(pred, target)\n        total += W_SSIM * ssim_loss(pred, target)\n\n        p_pad = F.pad(pred, (1, 1, 1, 1), mode=\"replicate\")\n        t_pad = F.pad(target, (1, 1, 1, 1), mode=\"replicate\")\n\n        total += W_SOBEL * (\n            charbonnier_loss(F.conv2d(p_pad, self.sobel_x), F.conv2d(t_pad, self.sobel_x))\n            + charbonnier_loss(F.conv2d(p_pad, self.sobel_y), F.conv2d(t_pad, self.sobel_y))\n        )\n        total += W_LAP * charbonnier_loss(F.conv2d(p_pad, self.lap), F.conv2d(t_pad, self.lap))\n\n        if allow_fft:\n            total += W_FFT * fft_spectrum_loss(pred, target, self.fft_mask)\n        return total\n\ndef alg_humility_penalty(pred, inp, is_id):\n    center = inp[:, 1:2]\n    per_sample = torch.mean(torch.abs(pred - center), dim=(1, 2, 3))\n    return torch.mean(per_sample * is_id * W_CHANGE_ID)\n\ndef low_t_edit_penalty(pred, inp, meta):\n    center = inp[:, 1:2]\n    t_norm = meta[:, 0]\n    low_mask = (t_norm <= LOW_T_THR).float()\n    per_sample = torch.mean(torch.abs(pred - center), dim=(1, 2, 3))\n    return torch.mean(per_sample * low_mask * W_LOW_T_EDIT)\n\ndef authority_tv_penalty(authority):\n    dy = torch.abs(authority[:, :, 1:, :] - authority[:, :, :-1, :]).mean()\n    dx = torch.abs(authority[:, :, :, 1:] - authority[:, :, :, :-1]).mean()\n    return (dx + dy) * W_AUTH_TV\n\n# ============================================================\n# Validation\n# ============================================================\n\n@torch.no_grad()\ndef eval_model_psnr(model, val_uids, series_root, target_shape=(64, 448, 448), max_uids=8):\n    model.eval()\n    scores = []\n    cache = VolumeLRU(max_items=2)\n\n    for uid in list(val_uids)[:max_uids]:\n        vol = cache.get(uid)\n        if vol is None:\n            vol = load_series_volume(uid, series_root, target_shape)\n            if vol is None:\n                continue\n            cache.put(uid, vol)\n\n        D = vol.shape[0]\n        for blur_level in [3, 8]:\n            for z in range(1, D - 1, 8):\n                cl = vol[z].astype(np.float32)\n                prev = vol[z - 1].astype(np.float32)\n                cent = vol[z].astype(np.float32)\n                next_ = vol[z + 1].astype(np.float32)\n\n                bp = mixed_poisson_gaussian(gaussian_psf_surrogate(prev, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bc = mixed_poisson_gaussian(gaussian_psf_surrogate(cent, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bn = mixed_poisson_gaussian(gaussian_psf_surrogate(next_, blur_level, DIFFUSION_ALPHA), \"quarter\")\n\n                t_norm = float(blur_level) / BLUR_LEVEL_MAX\n                inp_np = np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32)\n                inp_t = torch.from_numpy(inp_np).unsqueeze(0).to(device)\n\n                meta_np = np.array([[t_norm, 0.0, 0.0, 0.0]], dtype=np.float32)\n                meta_t = torch.from_numpy(meta_np).to(device)\n\n                with AMP_CTX():\n                    pred = model(inp_t, meta=meta_t)[0, 0].float().cpu().numpy()\n\n                mse = float(np.mean((pred - cl) ** 2))\n                scores.append(99.0 if mse <= 0 else 10.0 * math.log10(1.0 / mse))\n\n    return float(np.mean(scores)) if scores else None\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-14T22:13:53.307046Z","iopub.execute_input":"2026-03-14T22:13:53.307565Z","iopub.status.idle":"2026-03-14T22:13:58.336032Z","shell.execute_reply.started":"2026-03-14T22:13:53.307536Z","shell.execute_reply":"2026-03-14T22:13:58.335201Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Scientific-Method Training Design for the Synthetic-Supervision Pipeline\n\n## Why this training design matters scientifically\nIf a project claims that synthetic supervision is trustworthy, then the evaluation protocol itself must also be trustworthy.\n\nA common failure mode in ML projects is:\n1. train a model,\n2. peek at the final test set,\n3. change the method based on those failures,\n4. report the same test set as “independent evidence.”\n\nThat weakens the credibility of any claimed generalization.\n\nThis notebook avoids that mistake by using a **firewalled training design**:\n\n- a main training split,\n- a separate refinement split,\n- and an external evaluation stage kept conceptually distinct.\n\n## Core idea in one sentence\n> **Use internal data to improve the training process, but do not let the external evidence become part of the tuning loop.**\n\n## Why this matters for the methodology claim\nBecause the project is not only proposing a model.  \nIt is proposing a **way to create supervision from simulated degradations** and then argue that the supervision is meaningful.\n\nThat argument only works if the external evidence remains genuinely external.\n\n## What this section contributes to the final story\nThis training design supports the statement that later results—especially on Mayo—are being used as **validation of transfer**, not as hidden tuning feedback.\n","metadata":{}},{"cell_type":"code","source":"# =====================================================================\n# V-Ultimate scientific-method training cell\n# Two-stage training:\n#   Stage A: base training on TRAIN split\n#   Stage B: failure-driven refinement on REFINE split only\n#\n# Goal:\n#   improve the model using observed failure modes\n#   without leaking information from the external test set\n# =====================================================================\n\nimport os, gc, math, time, random, json\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport pydicom\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader, Sampler\nfrom contextlib import nullcontext\nimport torch.fft\n\n# -----------------------------\n# environment\n# -----------------------------\ntry:\n    cv2.setNumThreads(0)\nexcept Exception:\n    pass\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(\"Device:\", device)\n\nUSE_AMP = (device.type == \"cuda\")\nAMP_CTX = lambda: torch.amp.autocast(\"cuda\") if USE_AMP else nullcontext()\nscaler = torch.amp.GradScaler(\"cuda\") if USE_AMP else None\n\ntorch.backends.cudnn.benchmark = True\nif device.type == \"cuda\":\n    try:\n        torch.backends.cuda.matmul.allow_tf32 = True\n        torch.backends.cudnn.allow_tf32 = True\n        torch.set_float32_matmul_precision(\"high\")\n    except Exception:\n        pass\n\n# -----------------------------\n# require architecture\n# -----------------------------\nif \"DeblurUNet25D_Ultimate\" not in globals():\n    raise RuntimeError(\"请先运行新版 architecture cell（FiLM + authority map）。\")\n\n# ============================================================\n# config\n# ============================================================\n\n# --- paths ---\nRSNA_DATA_ROOT = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/series\"\nTRAIN_LOCALIZERS_CSV = \"/kaggle/input/rsna-intracranial-aneurysm-detection/train_localizers.csv\"\nMETA_CSV = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train.csv\"\n\nOUTDIR = \"/kaggle/working/vultimate_scientific\"\nos.makedirs(OUTDIR, exist_ok=True)\n\nOUT_TRAIN_UIDS  = os.path.join(OUTDIR, \"train_uids.csv\")\nOUT_REFINE_UIDS = os.path.join(OUTDIR, \"refine_uids.csv\")\nOUT_VAL_UIDS    = os.path.join(OUTDIR, \"val_uids.csv\")\nOUT_HARD_UIDS   = os.path.join(OUTDIR, \"hard_uids_stageB.csv\")\n\nSAVE_STAGEA_BEST = os.path.join(OUTDIR, \"deblur_stageA_best.pt\")\nSAVE_STAGEA_LAST = os.path.join(OUTDIR, \"deblur_stageA_last.pt\")\nSAVE_STAGEB_BEST = os.path.join(OUTDIR, \"deblur_stageB_best.pt\")\nSAVE_STAGEB_LAST = os.path.join(OUTDIR, \"deblur_stageB_last.pt\")\n\n# --- dataset size ---\nSEED = 2026\nN_TRAIN_UIDS  = 300\nN_REFINE_UIDS = 60\nN_VAL_UIDS    = 20\n\n# --- stage A training ---\nEPOCHS_STAGE_A = 18\nBATCHES_PER_EPOCH_A = 300\n\n# --- stage B refinement ---\nEPOCHS_STAGE_B = 8\nBATCHES_PER_EPOCH_B = 220\nN_HARD_UIDS = 24           # 从 refine 集里挑多少个 hardest cases\nHARD_UID_REPEAT = 3        # Stage B 中对 hard UID 做重复强化\n\n# --- optimizer ---\nBATCH_SIZE = 32\nNUM_WORKERS = 4\nLR_STAGE_A = 2e-4\nLR_STAGE_B = 8e-5\nWEIGHT_DECAY = 1e-4\nGRAD_CLIP = 1.0\n\n# --- image ---\nTARGET_D, TARGET_H, TARGET_W = 64, 448, 448\nPATCH_SIZE = 112\nPATCHES_PER_SLICE = 1\n\nHU_MIN, HU_MAX = -1024.0, 3072.0\nHU_RANGE = HU_MAX - HU_MIN\n\n# --- degradation ---\nDIFFUSION_ALPHA = 0.20\nBLUR_LEVELS = [0, 1, 3, 5, 8]\nBLUR_LEVEL_MAX = float(max(BLUR_LEVELS))\nBLUR_T_MAX = BLUR_LEVEL_MAX\n\nP_IDENTITY = 0.20\nENABLE_MOTION = True\nP_MOTION = 0.15\n\nPEAK_RANGE_QUARTER, SIGMA_E_QUARTER = (3000.0, 6000.0), (0.01, 0.02)\nPEAK_RANGE_EXTREME, SIGMA_E_EXTREME = (1000.0, 3000.0), (0.02, 0.04)\n\nREGIME_PROBS_STAGE_A = {\"clean\": 0.30, \"typical\": 0.45, \"hard\": 0.15, \"motion\": 0.10}\nREGIME_PROBS_STAGE_B = {\"clean\": 0.20, \"typical\": 0.35, \"hard\": 0.30, \"motion\": 0.15}\n\n# --- anatomy-aware crop ---\nANATOMY_REJECT_TRIES = 15\nANATOMY_MEAN_TH = 0.05\nANATOMY_STD_TH  = 0.02\nP_RANDOM_PATCH  = 0.10\n\n# --- loss weights ---\nW_CHARBONNIER = 1.0\nW_SSIM = 0.20\nW_SOBEL = 0.10\nW_LAP = 0.05\nW_FFT = 0.05\nFFT_FCUTOFF = 0.20\nFFT_ONLY_IF_T_LE = 10.0\n\nW_CHANGE_ID = 10.0\nW_LOW_T_EDIT = 2.0\nLOW_T_THR = 0.20\nW_AUTH_TV = 1e-3\n\n# stage B 强一点，专门压过修复\nW_LOW_T_EDIT_STAGE_B = 3.0\nW_AUTH_TV_STAGE_B = 2e-3\n\nRES_MIN, RES_MAX = 0.02, 0.15\n\n# --- hard-case mining ---\nMINE_EVAL_MAX_UIDS = None        # None = all refine uids\nMINE_STRIDE_Z = 8\nMINE_BLUR_LEVELS = [3, 8]\nHARD_SCORE_W_PSNR = 1.0\nHARD_SCORE_W_OVEREDIT = 2.0\n\n# ============================================================\n# seed\n# ============================================================\n\ndef seed_all(seed):\n    random.seed(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed_all(seed)\n\nseed_all(SEED)\n\ndef _normalize_probs(d):\n    s = float(sum(d.values()))\n    return {k: float(v) / s for k, v in d.items()}\n\nREGIME_PROBS_STAGE_A = _normalize_probs(REGIME_PROBS_STAGE_A)\nREGIME_PROBS_STAGE_B = _normalize_probs(REGIME_PROBS_STAGE_B)\n\n# ============================================================\n# split logic\n# scientific-method version:\n#   train      -> stage A base learning\n#   refine     -> failure analysis + stage B targeted improvement\n#   val        -> held-out internal validation only\n# external Mayo is untouched\n# ============================================================\n\ndef build_ct_uid_lists_3way(meta_csv, rsna_series_root, localizers_csv, n_train, n_refine, n_val, seed):\n    rsna_uids = set(\n        [u for u in os.listdir(rsna_series_root)\n         if os.path.isdir(os.path.join(rsna_series_root, u)) and not u.startswith(\".\")]\n    )\n\n    localizer_uids = set()\n    if localizers_csv and os.path.exists(localizers_csv):\n        df_loc = pd.read_csv(localizers_csv)\n        localizer_uids = set(df_loc[df_loc.columns[0]].astype(str).tolist())\n\n    meta = pd.read_csv(meta_csv)\n    meta_ct = meta[meta[\"Modality\"].isin({\"CT\", \"CTA\"})]\n\n    ct_candidates = [\n        u for u in set(meta_ct[\"SeriesInstanceUID\"].astype(str))\n        if u in rsna_uids and u not in localizer_uids\n    ]\n\n    rng = random.Random(seed)\n    rng.shuffle(ct_candidates)\n\n    total_need = n_train + n_refine + n_val\n    picked = ct_candidates[:min(total_need, len(ct_candidates))]\n\n    train_uids = picked[:min(n_train, len(picked))]\n    rem = picked[len(train_uids):]\n\n    refine_uids = rem[:min(n_refine, len(rem))]\n    rem2 = rem[len(refine_uids):]\n\n    val_uids = rem2[:min(n_val, len(rem2))]\n\n    return train_uids, refine_uids, val_uids\n\n# ============================================================\n# DICOM loading\n# ============================================================\n\ndef get_sorted_dicom_files(series_path):\n    files = [f for f in os.listdir(series_path) if not f.startswith(\".\")]\n    pairs, ok = [], True\n    for f in files:\n        try:\n            ds = pydicom.dcmread(os.path.join(series_path, f), stop_before_pixels=True, force=True)\n            if getattr(ds, \"InstanceNumber\", None) is None:\n                ok = False\n                break\n            pairs.append((int(ds.InstanceNumber), os.path.join(series_path, f)))\n        except Exception:\n            ok = False\n            break\n\n    if ok and len(pairs) == len(files):\n        return [p[1] for p in sorted(pairs, key=lambda x: x[0])]\n\n    return [os.path.join(series_path, f) for f in sorted(files)]\n\ndef load_series_volume(uid, series_root, target_shape=(64, 448, 448)):\n    series_path = os.path.join(series_root, uid)\n    if not os.path.isdir(series_path):\n        return None\n\n    dcm_files = get_sorted_dicom_files(series_path)\n    tD, tH, tW = target_shape\n    if len(dcm_files) < 10:\n        return None\n\n    if len(dcm_files) != tD:\n        idxs = np.linspace(0, len(dcm_files) - 1, tD).astype(int)\n        dcm_files = [dcm_files[i] for i in idxs]\n\n    slices = []\n    for fp in dcm_files:\n        try:\n            ds = pydicom.dcmread(fp, force=True)\n            hu = ds.pixel_array.astype(np.float32) * float(getattr(ds, \"RescaleSlope\", 1.0)) \\\n                 + float(getattr(ds, \"RescaleIntercept\", 0.0))\n            x = (np.clip(hu, HU_MIN, HU_MAX) - HU_MIN) / HU_RANGE\n            x = cv2.resize(x, (tW, tH), interpolation=cv2.INTER_LINEAR)\n            slices.append(x.astype(np.float32))\n        except Exception:\n            continue\n\n    if len(slices) < int(0.8 * tD):\n        return None\n\n    while len(slices) < tD:\n        slices.append(slices[-1].copy())\n\n    return np.stack(slices[:tD], axis=0).astype(np.float32)\n\nclass VolumeLRU:\n    def __init__(self, max_items=12):\n        self.max_items = int(max_items)\n        self.cache = {}\n        self.order = []\n\n    def get(self, key):\n        if key not in self.cache:\n            return None\n        self.order.remove(key)\n        self.order.append(key)\n        return self.cache[key]\n\n    def put(self, key, value):\n        if key in self.cache:\n            self.order.remove(key)\n        self.cache[key] = value\n        self.order.append(key)\n        if len(self.order) > self.max_items:\n            old = self.order.pop(0)\n            self.cache.pop(old, None)\n\n# ============================================================\n# physics degradation\n# ============================================================\n\ndef gaussian_psf_surrogate(img01, blur_level, alpha=0.20):\n    if blur_level <= 0:\n        return img01\n    sigma = math.sqrt(max(1e-8, 2.0 * alpha * float(blur_level)))\n    return np.clip(\n        cv2.GaussianBlur(img01, (0, 0), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REPLICATE),\n        0.0, 1.0\n    )\n\ndef motion_artifact_surrogate(img01, length=None, angle=None):\n    length = length or random.choice([3, 5, 7, 9, 11])\n    if length <= 1:\n        return img01\n    angle = angle or random.uniform(0, 180)\n\n    k = np.zeros((length, length), dtype=np.float32)\n    c = length // 2\n    cos_a, sin_a = np.cos(np.radians(angle)), np.sin(np.radians(angle))\n    for i in range(length):\n        x = int(c + (i - c) * cos_a)\n        y = int(c + (i - c) * sin_a)\n        if 0 <= x < length and 0 <= y < length:\n            k[y, x] = 1.0\n    if k.sum() > 0:\n        k /= k.sum()\n\n    return np.clip(cv2.filter2D(img01, -1, k, borderType=cv2.BORDER_REPLICATE), 0.0, 1.0)\n\ndef mixed_poisson_gaussian(img01, mode=\"quarter\"):\n    if mode == \"clean\":\n        return img01\n\n    peak_rng, sigma_rng = (\n        (PEAK_RANGE_EXTREME, SIGMA_E_EXTREME) if mode == \"extreme\"\n        else (PEAK_RANGE_QUARTER, SIGMA_E_QUARTER)\n    )\n    peak = random.uniform(*peak_rng)\n    sigma_e = random.uniform(*sigma_rng)\n    noisy_p = np.random.poisson(np.clip(img01 * peak, 0, None)).astype(np.float32) / peak\n    noisy_g = np.random.randn(*img01.shape).astype(np.float32) * sigma_e\n    return np.clip(noisy_p + noisy_g, 0.0, 1.0)\n\ndef _choice_weighted(prob_dict):\n    r = random.random()\n    acc = 0.0\n    for k, p in prob_dict.items():\n        acc += p\n        if r <= acc:\n            return k\n    return list(prob_dict.keys())[-1]\n\ndef _sample_regime_params(prob_dict):\n    reg = _choice_weighted(prob_dict)\n    if reg == \"clean\":\n        return 0, \"clean\", False\n    if reg == \"typical\":\n        return random.choice([1, 3, 5]), \"quarter\", False\n    if reg == \"hard\":\n        return random.choice([3, 5, 8]), \"extreme\", False\n    if reg == \"motion\":\n        return random.choice([1, 3, 5]), random.choice([\"quarter\", \"extreme\"]), True\n    return 3, \"quarter\", False\n\n# ============================================================\n# datasets\n# ============================================================\n\nclass CTDeblur25D(Dataset):\n    def __init__(\n        self,\n        uids,\n        series_root,\n        target_shape=(64, 448, 448),\n        patch_size=112,\n        patches_per_slice=1,\n        regime_probs=None,\n        p_identity=0.20,\n        enable_motion=True,\n        p_motion=0.15,\n        hard_uid_set=None,\n        hard_uid_repeat=1,\n    ):\n        self.uids = list(uids)\n        self.series_root = series_root\n        self.target_shape = target_shape\n        self.patch_size = int(patch_size)\n        self.patches_per_slice = int(patches_per_slice)\n        self.regime_probs = _normalize_probs(regime_probs or REGIME_PROBS_STAGE_A)\n        self.p_identity = float(p_identity)\n        self.enable_motion = bool(enable_motion)\n        self.p_motion = float(p_motion)\n        self.cache = VolumeLRU(max_items=8)\n\n        self.hard_uid_set = set(hard_uid_set or [])\n        self.hard_uid_repeat = int(max(1, hard_uid_repeat))\n\n        items = []\n        for ui in range(len(self.uids)):\n            uid = self.uids[ui]\n            rep = self.hard_uid_repeat if uid in self.hard_uid_set else 1\n            for _ in range(rep):\n                items.extend([(ui, z) for z in range(1, target_shape[0] - 1)])\n        self.items = items\n\n    def __len__(self):\n        return len(self.items)\n\n    def _sample_patch_xy(self, cent, ps):\n        H, W = cent.shape\n        if random.random() < P_RANDOM_PATCH:\n            return np.random.randint(0, H - ps + 1), np.random.randint(0, W - ps + 1)\n\n        for _ in range(ANATOMY_REJECT_TRIES):\n            y = np.random.randint(0, H - ps + 1)\n            x = np.random.randint(0, W - ps + 1)\n            patch = cent[y:y+ps, x:x+ps]\n            if patch.mean() > ANATOMY_MEAN_TH and patch.std() > ANATOMY_STD_TH:\n                return y, x\n\n        return (H - ps) // 2, (W - ps) // 2\n\n    def __getitem__(self, idx):\n        ui, z = self.items[idx]\n        uid = self.uids[ui]\n\n        vol = self.cache.get(uid)\n        if vol is None:\n            vol = load_series_volume(uid, self.series_root, self.target_shape)\n            if vol is None:\n                return self.__getitem__(random.randint(0, len(self.items) - 1))\n            self.cache.put(uid, vol)\n\n        ps = self.patch_size\n        y, x = self._sample_patch_xy(vol[z], ps)\n\n        clean = vol[z][y:y+ps, x:x+ps].copy()\n        pp = vol[z-1][y:y+ps, x:x+ps].copy()\n        cc = vol[z][y:y+ps, x:x+ps].copy()\n        nn_ = vol[z+1][y:y+ps, x:x+ps].copy()\n\n        is_identity = 0.0\n        if random.random() < self.p_identity:\n            blur_level, dose_mode, do_motion, is_identity = 0, \"clean\", False, 1.0\n        else:\n            blur_level, dose_mode, do_motion = _sample_regime_params(self.regime_probs)\n            do_motion = do_motion and self.enable_motion and (random.random() < self.p_motion)\n\n        bp = gaussian_psf_surrogate(pp, blur_level, alpha=DIFFUSION_ALPHA)\n        bc = gaussian_psf_surrogate(cc, blur_level, alpha=DIFFUSION_ALPHA)\n        bn = gaussian_psf_surrogate(nn_, blur_level, alpha=DIFFUSION_ALPHA)\n\n        if do_motion:\n            L = random.choice([3, 5, 7, 9, 11])\n            A = random.uniform(0, 180)\n            bp = motion_artifact_surrogate(bp, L, A)\n            bc = motion_artifact_surrogate(bc, L, A)\n            bn = motion_artifact_surrogate(bn, L, A)\n\n        if dose_mode != \"clean\":\n            bp = mixed_poisson_gaussian(bp, dose_mode)\n            bc = mixed_poisson_gaussian(bc, dose_mode)\n            bn = mixed_poisson_gaussian(bn, dose_mode)\n\n        cp = clean.copy()\n        if random.random() > 0.5:\n            cp, bp, bc, bn = cp[::-1].copy(), bp[::-1].copy(), bc[::-1].copy(), bn[::-1].copy()\n        if random.random() > 0.5:\n            cp, bp, bc, bn = cp[:, ::-1].copy(), bp[:, ::-1].copy(), bc[:, ::-1].copy(), bn[:, ::-1].copy()\n        k = random.randint(0, 3)\n        if k > 0:\n            cp, bp, bc, bn = np.rot90(cp, k).copy(), np.rot90(bp, k).copy(), np.rot90(bc, k).copy(), np.rot90(bn, k).copy()\n\n        t_norm = float(blur_level) / BLUR_LEVEL_MAX if BLUR_LEVEL_MAX > 0 else 0.0\n\n        inp = np.stack([bp, bc, bn, np.full_like(bc, t_norm, dtype=np.float32)], axis=0).astype(np.float32)\n        tgt = cp[np.newaxis, ...].astype(np.float32)\n\n        meta = np.array([\n            t_norm,\n            1.0 if do_motion else 0.0,\n            1.0 if dose_mode == \"clean\" else 0.0,\n            is_identity\n        ], dtype=np.float32)\n\n        return (\n            torch.from_numpy(inp).float(),\n            torch.from_numpy(tgt).float(),\n            torch.from_numpy(meta).float()\n        )\n\nclass UIDBatchSampler(Sampler):\n    def __init__(self, dataset, batch_size, seed=42, batches_per_epoch=None):\n        self.dataset = dataset\n        self.batch_size = int(batch_size)\n        self.rng = random.Random(seed)\n        self.by_ui = {}\n\n        for idx, (ui, z) in enumerate(dataset.items):\n            self.by_ui.setdefault(ui, []).append(idx)\n\n        self.ui_keys = list(self.by_ui.keys())\n        self.batches_per_epoch = int(batches_per_epoch) if batches_per_epoch else len(dataset) // self.batch_size\n\n    def __len__(self):\n        return self.batches_per_epoch\n\n    def __iter__(self):\n        for _ in range(self.batches_per_epoch):\n            ui = self.rng.choice(self.ui_keys)\n            pool = self.by_ui[ui]\n            if len(pool) >= self.batch_size:\n                yield self.rng.sample(pool, self.batch_size)\n            else:\n                yield [self.rng.choice(pool) for _ in range(self.batch_size)]\n\n# ============================================================\n# loss\n# ============================================================\n\ndef charbonnier_loss(pred, target, eps=1e-3):\n    return torch.mean(torch.sqrt((pred - target) ** 2 + eps ** 2))\n\ndef fft_spectrum_loss(pred, target, fft_mask):\n    pred_fft = torch.fft.rfft2(pred.float(), dim=(-2, -1), norm=\"ortho\")\n    tgt_fft = torch.fft.rfft2(target.float(), dim=(-2, -1), norm=\"ortho\")\n    m = fft_mask.view(1, 1, *fft_mask.shape)\n    return charbonnier_loss(torch.abs(pred_fft) * m, torch.abs(tgt_fft) * m)\n\ndef ssim_loss(pred, target, window_size=11):\n    C1, C2 = 0.01 ** 2, 0.03 ** 2\n    pad = window_size // 2\n\n    mu_x = F.avg_pool2d(pred, window_size, stride=1, padding=pad)\n    mu_y = F.avg_pool2d(target, window_size, stride=1, padding=pad)\n\n    sigma_x2 = F.avg_pool2d(pred ** 2, window_size, stride=1, padding=pad) - mu_x ** 2\n    sigma_y2 = F.avg_pool2d(target ** 2, window_size, stride=1, padding=pad) - mu_y ** 2\n    sigma_xy = F.avg_pool2d(pred * target, window_size, stride=1, padding=pad) - mu_x * mu_y\n\n    ssim_map = ((2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)) / (\n        (mu_x ** 2 + mu_y ** 2 + C1) * (sigma_x2 + sigma_y2 + C2) + 1e-8\n    )\n    return 1.0 - ssim_map.mean()\n\ndef _make_fft_mask(H, W, fcut=0.20, device=\"cpu\"):\n    fy = torch.fft.fftfreq(H, d=1.0, device=device).view(H, 1).abs()\n    fx = torch.fft.rfftfreq(W, d=1.0, device=device).view(1, W // 2 + 1).abs()\n    return (torch.sqrt(fx * fx + fy * fy) >= fcut).float()\n\nclass UltimatePhysicsLoss(nn.Module):\n    def __init__(self, patch_size=112):\n        super().__init__()\n        self.sobel_x = torch.tensor([[[-1., 0., 1.], [-2., 0., 2.], [-1., 0., 1.]]]).view(1, 1, 3, 3).to(device)\n        self.sobel_y = torch.tensor([[[-1., -2., -1.], [0., 0., 0.], [1., 2., 1.]]]).view(1, 1, 3, 3).to(device)\n        self.lap = torch.tensor([[[0., 1., 0.], [1., -4., 1.], [0., 1., 0.]]]).view(1, 1, 3, 3).to(device)\n        self.register_buffer(\"fft_mask\", _make_fft_mask(patch_size, patch_size, fcut=FFT_FCUTOFF, device=device))\n\n    def forward(self, pred, target, allow_fft=False):\n        total = W_CHARBONNIER * charbonnier_loss(pred, target)\n        total += W_SSIM * ssim_loss(pred, target)\n\n        p_pad = F.pad(pred, (1, 1, 1, 1), mode=\"replicate\")\n        t_pad = F.pad(target, (1, 1, 1, 1), mode=\"replicate\")\n\n        total += W_SOBEL * (\n            charbonnier_loss(F.conv2d(p_pad, self.sobel_x), F.conv2d(t_pad, self.sobel_x))\n            + charbonnier_loss(F.conv2d(p_pad, self.sobel_y), F.conv2d(t_pad, self.sobel_y))\n        )\n        total += W_LAP * charbonnier_loss(F.conv2d(p_pad, self.lap), F.conv2d(t_pad, self.lap))\n\n        if allow_fft:\n            total += W_FFT * fft_spectrum_loss(pred, target, self.fft_mask)\n        return total\n\ndef alg_humility_penalty(pred, inp, is_id, weight_change_id):\n    center = inp[:, 1:2]\n    per_sample = torch.mean(torch.abs(pred - center), dim=(1, 2, 3))\n    return torch.mean(per_sample * is_id * weight_change_id)\n\ndef low_t_edit_penalty(pred, inp, meta, weight_low_t):\n    center = inp[:, 1:2]\n    t_norm = meta[:, 0]\n    low_mask = (t_norm <= LOW_T_THR).float()\n    per_sample = torch.mean(torch.abs(pred - center), dim=(1, 2, 3))\n    return torch.mean(per_sample * low_mask * weight_low_t)\n\ndef authority_tv_penalty(authority, weight_auth_tv):\n    dy = torch.abs(authority[:, :, 1:, :] - authority[:, :, :-1, :]).mean()\n    dx = torch.abs(authority[:, :, :, 1:] - authority[:, :, :, :-1]).mean()\n    return (dx + dy) * weight_auth_tv\n\n# ============================================================\n# validation\n# ============================================================\n\n@torch.no_grad()\ndef eval_model_psnr(model, val_uids, series_root, target_shape=(64, 448, 448), max_uids=8):\n    model.eval()\n    scores = []\n    cache = VolumeLRU(max_items=2)\n\n    for uid in list(val_uids)[:max_uids]:\n        vol = cache.get(uid)\n        if vol is None:\n            vol = load_series_volume(uid, series_root, target_shape)\n            if vol is None:\n                continue\n            cache.put(uid, vol)\n\n        D = vol.shape[0]\n        for blur_level in [3, 8]:\n            for z in range(1, D - 1, 8):\n                cl = vol[z].astype(np.float32)\n                prev = vol[z - 1].astype(np.float32)\n                cent = vol[z].astype(np.float32)\n                next_ = vol[z + 1].astype(np.float32)\n\n                bp = mixed_poisson_gaussian(gaussian_psf_surrogate(prev, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bc = mixed_poisson_gaussian(gaussian_psf_surrogate(cent, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bn = mixed_poisson_gaussian(gaussian_psf_surrogate(next_, blur_level, DIFFUSION_ALPHA), \"quarter\")\n\n                t_norm = float(blur_level) / BLUR_LEVEL_MAX\n                inp_np = np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32)\n                meta_np = np.array([[t_norm, 0.0, 0.0, 0.0]], dtype=np.float32)\n\n                inp_t = torch.from_numpy(inp_np).unsqueeze(0).to(device)\n                meta_t = torch.from_numpy(meta_np).to(device)\n\n                with AMP_CTX():\n                    pred = model(inp_t, meta=meta_t)[0, 0].float().cpu().numpy()\n\n                mse = float(np.mean((pred - cl) ** 2))\n                scores.append(99.0 if mse <= 0 else 10.0 * math.log10(1.0 / mse))\n\n    return float(np.mean(scores)) if scores else None\n\n# ============================================================\n# hard-case mining\n# only on refine split\n# ============================================================\n\n@torch.no_grad()\ndef mine_hard_uids(model, refine_uids, series_root, target_shape=(64, 448, 448), max_uids=None):\n    model.eval()\n    cache = VolumeLRU(max_items=2)\n    rows = []\n\n    use_uids = list(refine_uids)\n    if max_uids is not None:\n        use_uids = use_uids[:max_uids]\n\n    for uid in use_uids:\n        vol = cache.get(uid)\n        if vol is None:\n            vol = load_series_volume(uid, series_root, target_shape)\n            if vol is None:\n                continue\n            cache.put(uid, vol)\n\n        D = vol.shape[0]\n        psnr_scores = []\n        overedit_scores = []\n\n        for blur_level in MINE_BLUR_LEVELS:\n            for z in range(1, D - 1, MINE_STRIDE_Z):\n                cl = vol[z].astype(np.float32)\n                prev = vol[z - 1].astype(np.float32)\n                cent = vol[z].astype(np.float32)\n                next_ = vol[z + 1].astype(np.float32)\n\n                bp = mixed_poisson_gaussian(gaussian_psf_surrogate(prev, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bc = mixed_poisson_gaussian(gaussian_psf_surrogate(cent, blur_level, DIFFUSION_ALPHA), \"quarter\")\n                bn = mixed_poisson_gaussian(gaussian_psf_surrogate(next_, blur_level, DIFFUSION_ALPHA), \"quarter\")\n\n                t_norm = float(blur_level) / BLUR_LEVEL_MAX\n                inp_np = np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32)\n                meta_np = np.array([[t_norm, 0.0, 0.0, 0.0]], dtype=np.float32)\n\n                inp_t = torch.from_numpy(inp_np).unsqueeze(0).to(device)\n                meta_t = torch.from_numpy(meta_np).to(device)\n\n                with AMP_CTX():\n                    pred, aux = model(inp_t, meta=meta_t, return_aux=True)\n                    pred_np = pred[0, 0].float().cpu().numpy()\n\n                mse = float(np.mean((pred_np - cl) ** 2))\n                psnr = 99.0 if mse <= 0 else 10.0 * math.log10(1.0 / mse)\n                psnr_scores.append(psnr)\n\n                overedit = float(np.mean(np.abs(pred_np - bc)))\n                overedit_scores.append(overedit)\n\n        if len(psnr_scores) == 0:\n            continue\n\n        psnr_mean = float(np.mean(psnr_scores))\n        overedit_mean = float(np.mean(overedit_scores))\n\n        rows.append({\n            \"uid\": uid,\n            \"psnr_mean\": psnr_mean,\n            \"overedit_mean\": overedit_mean,\n        })\n\n    df_hard = pd.DataFrame(rows)\n    if len(df_hard) == 0:\n        return [], df_hard\n\n    # 分数越大越难：低 PSNR + 高 overedit\n    psnr_norm = (df_hard[\"psnr_mean\"].max() - df_hard[\"psnr_mean\"])\n    if psnr_norm.max() > 0:\n        psnr_norm = psnr_norm / (psnr_norm.max() + 1e-8)\n\n    over_norm = df_hard[\"overedit_mean\"]\n    if over_norm.max() > 0:\n        over_norm = over_norm / (over_norm.max() + 1e-8)\n\n    df_hard[\"hard_score\"] = HARD_SCORE_W_PSNR * psnr_norm + HARD_SCORE_W_OVEREDIT * over_norm\n    df_hard = df_hard.sort_values(\"hard_score\", ascending=False).reset_index(drop=True)\n\n    hard_uids = df_hard[\"uid\"].tolist()[:min(N_HARD_UIDS, len(df_hard))]\n    return hard_uids, df_hard\n\n# ============================================================\n# generic training loop\n# ============================================================\n\ndef run_training_stage(\n    stage_name,\n    model,\n    train_loader,\n    val_uids,\n    save_best,\n    save_last,\n    epochs,\n    lr,\n    weight_low_t,\n    weight_auth_tv,\n    extra_metadata=None,\n):\n    opt = torch.optim.AdamW(model.parameters(), lr=lr, weight_decay=WEIGHT_DECAY)\n    sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=epochs)\n    crit = UltimatePhysicsLoss(patch_size=PATCH_SIZE).to(device)\n\n    best_psnr = -1.0\n    hist = []\n\n    print(f\"\\n=== {stage_name} ===\")\n    print(\"lr:\", lr, \"| epochs:\", epochs)\n\n    for ep in range(1, epochs + 1):\n        model.train()\n        losses = []\n        t0 = time.time()\n\n        for b, (inp, tgt, meta) in enumerate(train_loader, 1):\n            inp = inp.to(device, non_blocking=True)\n            tgt = tgt.to(device, non_blocking=True)\n            meta = meta.to(device, non_blocking=True)\n\n            opt.zero_grad(set_to_none=True)\n\n            with AMP_CTX():\n                pred, aux = model(inp, meta=meta, return_aux=True)\n\n                t_scalar = meta[:, 0].mean().item() * BLUR_LEVEL_MAX\n                is_motion = meta[:, 1]\n                is_id = meta[:, 3]\n\n                allow_fft = (t_scalar <= FFT_ONLY_IF_T_LE) and not bool((is_motion > 0.5).any().item())\n\n                loss_main = crit(pred, tgt, allow_fft=allow_fft)\n                loss_id = alg_humility_penalty(pred, inp, is_id, W_CHANGE_ID)\n                loss_lowt = low_t_edit_penalty(pred, inp, meta, weight_low_t)\n                loss_auth = authority_tv_penalty(aux[\"authority\"], weight_auth_tv)\n\n                loss = loss_main + loss_id + loss_lowt + loss_auth\n\n            if USE_AMP:\n                scaler.scale(loss).backward()\n                scaler.unscale_(opt)\n                nn.utils.clip_grad_norm_(model.parameters(), GRAD_CLIP)\n                scaler.step(opt)\n                scaler.update()\n            else:\n                loss.backward()\n                nn.utils.clip_grad_norm_(model.parameters(), GRAD_CLIP)\n                opt.step()\n\n            losses.append(float(loss.item()))\n\n            if b == 1 or b % 50 == 0 or b == len(train_loader):\n                print(\n                    f\"  [{stage_name} | Epoch {ep:02d} | Batch {b:03d}/{len(train_loader)}] \"\n                    f\"Loss={np.mean(losses[-20:]):.4f} \"\n                    f\"(main={float(loss_main.item()):.4f}, id={float(loss_id.item()):.4f}, \"\n                    f\"lowt={float(loss_lowt.item()):.4f}, auth={float(loss_auth.item()):.6f})\"\n                )\n\n        sched.step()\n        epoch_loss = float(np.mean(losses)) if len(losses) else np.nan\n        print(f\"{stage_name} Epoch {ep:02d} done | Loss={epoch_loss:.5f} | Time={(time.time()-t0)/60:.1f} min\")\n\n        row = {\"stage\": stage_name, \"epoch\": ep, \"train_loss\": epoch_loss}\n\n        if ep % 3 == 0 or ep == epochs:\n            psnr_val = eval_model_psnr(model, val_uids, RSNA_DATA_ROOT, (TARGET_D, TARGET_H, TARGET_W), max_uids=8)\n            row[\"val_psnr\"] = psnr_val\n\n            if psnr_val is not None:\n                if psnr_val > best_psnr:\n                    best_psnr = psnr_val\n                    pack = {\n                        \"model\": model.state_dict(),\n                        \"epoch\": ep,\n                        \"best_val_psnr\": best_psnr,\n                        \"stage\": stage_name,\n                        \"extra_metadata\": extra_metadata or {},\n                    }\n                    torch.save(pack, save_best)\n                    print(f\"  [Val] PSNR={psnr_val:.2f} dB ★ NEW BEST\")\n                else:\n                    print(f\"  [Val] PSNR={psnr_val:.2f} dB\")\n\n        hist.append(row)\n\n    final_pack = {\n        \"model\": model.state_dict(),\n        \"epoch\": epochs,\n        \"best_val_psnr\": best_psnr,\n        \"stage\": stage_name,\n        \"extra_metadata\": extra_metadata or {},\n    }\n    torch.save(final_pack, save_last)\n\n    return best_psnr, pd.DataFrame(hist)\n\n# ============================================================\n# 1) split data\n# ============================================================\n\nprint(\"\\n=== Step 1. Build scientific train/refine/val split ===\")\ntrain_uids, refine_uids, val_uids = build_ct_uid_lists_3way(\n    META_CSV,\n    RSNA_DATA_ROOT,\n    TRAIN_LOCALIZERS_CSV,\n    N_TRAIN_UIDS,\n    N_REFINE_UIDS,\n    N_VAL_UIDS,\n    SEED,\n)\n\npd.DataFrame({\"SeriesInstanceUID\": train_uids}).to_csv(OUT_TRAIN_UIDS, index=False)\npd.DataFrame({\"SeriesInstanceUID\": refine_uids}).to_csv(OUT_REFINE_UIDS, index=False)\npd.DataFrame({\"SeriesInstanceUID\": val_uids}).to_csv(OUT_VAL_UIDS, index=False)\n\nprint(f\"train={len(train_uids)} | refine={len(refine_uids)} | val={len(val_uids)}\")\n\n# ============================================================\n# 2) Stage A dataset\n# ============================================================\n\nprint(\"\\n=== Step 2. Stage A dataset ===\")\ntrain_ds_A = CTDeblur25D(\n    train_uids,\n    RSNA_DATA_ROOT,\n    (TARGET_D, TARGET_H, TARGET_W),\n    PATCH_SIZE,\n    PATCHES_PER_SLICE,\n    regime_probs=REGIME_PROBS_STAGE_A,\n    p_identity=P_IDENTITY,\n    enable_motion=ENABLE_MOTION,\n    p_motion=P_MOTION,\n    hard_uid_set=None,\n    hard_uid_repeat=1,\n)\n\ntrain_loader_A = DataLoader(\n    train_ds_A,\n    batch_sampler=UIDBatchSampler(train_ds_A, BATCH_SIZE, SEED, BATCHES_PER_EPOCH_A),\n    num_workers=NUM_WORKERS,\n    pin_memory=True,\n)\n\n# ============================================================\n# 3) init model\n# ============================================================\n\nprint(\"\\n=== Step 3. Initialize model ===\")\nmodel = DeblurUNet25D_Ultimate(\n    in_ch=4,\n    out_ch=1,\n    base=32,\n    res_min=RES_MIN,\n    res_max=RES_MAX,\n    meta_dim=4,\n    film_hidden=64,\n    authority_bias_init=2.0,\n).to(device)\n\nprint(\"Trainable params:\", f\"{sum(p.numel() for p in model.parameters() if p.requires_grad):,}\")\n\n# ============================================================\n# 4) Stage A: base training\n# ============================================================\n\nstageA_meta = {\n    \"seed\": SEED,\n    \"train_uids\": train_uids,\n    \"refine_uids\": refine_uids,\n    \"val_uids\": val_uids,\n    \"regime_probs\": REGIME_PROBS_STAGE_A,\n    \"phase\": \"base_training\",\n}\n\nbestA, histA = run_training_stage(\n    stage_name=\"StageA_Base\",\n    model=model,\n    train_loader=train_loader_A,\n    val_uids=val_uids,\n    save_best=SAVE_STAGEA_BEST,\n    save_last=SAVE_STAGEA_LAST,\n    epochs=EPOCHS_STAGE_A,\n    lr=LR_STAGE_A,\n    weight_low_t=W_LOW_T_EDIT,\n    weight_auth_tv=W_AUTH_TV,\n    extra_metadata=stageA_meta,\n)\n\nhistA.to_csv(os.path.join(OUTDIR, \"history_stageA.csv\"), index=False)\n\nprint(\"\\nStage A best val PSNR:\", bestA)\n\n# ============================================================\n# 5) Failure analysis on REFINE split only\n# ============================================================\n\nprint(\"\\n=== Step 5. Mine hard cases on refine split only ===\")\nif os.path.exists(SAVE_STAGEA_BEST):\n    packA = torch.load(SAVE_STAGEA_BEST, map_location=\"cpu\")\n    model.load_state_dict(packA[\"model\"], strict=True)\n    model = model.to(device).eval()\n\nhard_uids, df_hard = mine_hard_uids(\n    model,\n    refine_uids,\n    RSNA_DATA_ROOT,\n    (TARGET_D, TARGET_H, TARGET_W),\n    max_uids=MINE_EVAL_MAX_UIDS,\n)\n\ndf_hard.to_csv(os.path.join(OUTDIR, \"refine_uid_difficulty.csv\"), index=False)\npd.DataFrame({\"SeriesInstanceUID\": hard_uids}).to_csv(OUT_HARD_UIDS, index=False)\n\nprint(f\"Hard UIDs selected for Stage B: {len(hard_uids)}\")\ndisplay(df_hard.head(20))\n\n# ============================================================\n# 6) Stage B dataset\n# only uses:\n#   - original TRAIN split\n#   - failure-derived hard UID list from REFINE split analysis\n# no external test leakage\n# ============================================================\n\nprint(\"\\n=== Step 6. Stage B refinement dataset ===\")\nstageB_uids = list(train_uids) + list(hard_uids)\n\ntrain_ds_B = CTDeblur25D(\n    stageB_uids,\n    RSNA_DATA_ROOT,\n    (TARGET_D, TARGET_H, TARGET_W),\n    PATCH_SIZE,\n    PATCHES_PER_SLICE,\n    regime_probs=REGIME_PROBS_STAGE_B,\n    p_identity=P_IDENTITY,\n    enable_motion=ENABLE_MOTION,\n    p_motion=P_MOTION,\n    hard_uid_set=set(hard_uids),\n    hard_uid_repeat=HARD_UID_REPEAT,\n)\n\ntrain_loader_B = DataLoader(\n    train_ds_B,\n    batch_sampler=UIDBatchSampler(train_ds_B, BATCH_SIZE, SEED + 101, BATCHES_PER_EPOCH_B),\n    num_workers=NUM_WORKERS,\n    pin_memory=True,\n)\n\nprint(f\"Stage B total UID list = {len(stageB_uids)} (train + hard refine)\")\nprint(f\"Hard UID repeat = {HARD_UID_REPEAT}\")\n\n# ============================================================\n# 7) load Stage A best and do Stage B refinement\n# ============================================================\n\nprint(\"\\n=== Step 7. Stage B refinement ===\")\nif os.path.exists(SAVE_STAGEA_BEST):\n    packA = torch.load(SAVE_STAGEA_BEST, map_location=\"cpu\")\n    model.load_state_dict(packA[\"model\"], strict=True)\n    model = model.to(device)\n\nstageB_meta = {\n    \"seed\": SEED,\n    \"train_uids\": train_uids,\n    \"refine_uids\": refine_uids,\n    \"val_uids\": val_uids,\n    \"hard_uids\": hard_uids,\n    \"regime_probs\": REGIME_PROBS_STAGE_B,\n    \"phase\": \"failure_driven_refinement\",\n    \"note\": \"hard cases mined only from refine split, not from external test set\",\n}\n\nbestB, histB = run_training_stage(\n    stage_name=\"StageB_Refine\",\n    model=model,\n    train_loader=train_loader_B,\n    val_uids=val_uids,\n    save_best=SAVE_STAGEB_BEST,\n    save_last=SAVE_STAGEB_LAST,\n    epochs=EPOCHS_STAGE_B,\n    lr=LR_STAGE_B,\n    weight_low_t=W_LOW_T_EDIT_STAGE_B,\n    weight_auth_tv=W_AUTH_TV_STAGE_B,\n    extra_metadata=stageB_meta,\n)\n\nhistB.to_csv(os.path.join(OUTDIR, \"history_stageB.csv\"), index=False)\n\n# ============================================================\n# 8) final summary\n# ============================================================\n\nsummary = {\n    \"train_uids\": len(train_uids),\n    \"refine_uids\": len(refine_uids),\n    \"val_uids\": len(val_uids),\n    \"hard_uids_stageB\": len(hard_uids),\n    \"best_stageA_val_psnr\": bestA,\n    \"best_stageB_val_psnr\": bestB,\n    \"stageA_best_ckpt\": SAVE_STAGEA_BEST,\n    \"stageB_best_ckpt\": SAVE_STAGEB_BEST,\n}\n\nwith open(os.path.join(OUTDIR, \"training_summary.json\"), \"w\") as f:\n    json.dump(summary, f, indent=2)\n\nprint(\"\\n\" + \"=\" * 90)\nprint(\"Scientific-method training finished\")\nprint(\"=\" * 90)\nfor k, v in summary.items():\n    print(f\"{k}: {v}\")\n\nprint(\"\\nSaved files:\")\nfor p in sorted(os.listdir(OUTDIR)):\n    print(\" -\", os.path.join(OUTDIR, p))\n\nprint(\"\\n✅ Weight File for later evaluation：\")\nprint(\"   \", SAVE_STAGEB_BEST if os.path.exists(SAVE_STAGEB_BEST) else SAVE_STAGEA_BEST)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-14T22:13:58.337425Z","iopub.execute_input":"2026-03-14T22:13:58.337836Z","iopub.status.idle":"2026-03-15T00:00:47.440469Z","shell.execute_reply.started":"2026-03-14T22:13:58.337813Z","shell.execute_reply":"2026-03-15T00:00:47.438665Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nBEST_CKPT_FOR_EVAL = SAVE_STAGEB_BEST if os.path.exists(SAVE_STAGEB_BEST) else SAVE_STAGEA_BEST\nprint(\"BEST_CKPT_FOR_EVAL =\", BEST_CKPT_FOR_EVAL)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T00:00:47.445734Z","iopub.execute_input":"2026-03-15T00:00:47.445976Z","iopub.status.idle":"2026-03-15T00:00:47.450965Z","shell.execute_reply.started":"2026-03-15T00:00:47.445951Z","shell.execute_reply":"2026-03-15T00:00:47.4502Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# Downstream Tests\n\nThese downstream experiments are not included simply to show that the images look cleaner. Their purpose is to test the **methodology** from several complementary angles.\n\nThe central question of this project is:\n\n> **When true paired data are unavailable, can physics-grounded synthetic supervision train a restoration model that is clinically useful, externally transferable, and structurally trustworthy?**\n\nA single metric cannot answer that question.  \nA model may score well on an internal validation set yet fail on real data, or improve pixel-level quality while making unsafe or anatomically meaningless edits.  \nFor that reason, the downstream section is organized into multiple evidence blocks, each addressing a different reviewer concern.\n\n---\n\n## 1. Clinical Rescue Matrix (CRM)\n\n### What it measures\nCRM evaluates whether the restoration model can **recover clinically relevant signal** that was suppressed by degradation on a strict held-out OOD CT/CTA set.\n\nInstead of asking only whether the image becomes smoother or sharper, CRM asks:\n\n> **Can the model restore information that matters for a downstream clinical proxy task?**\n\n### Why it matters\nThis is the most direct test of **functional usefulness**.  \nIf the synthetic supervision is genuinely valuable, the model should learn more than denoising or smoothing; it should recover image content that helps preserve or recover clinically meaningful signal.\n\n### Core questions it answers\n- Does the model help on a clinically motivated downstream proxy?\n- Does it outperform or complement traditional image filters?\n- Does it introduce harmful changes?\n- Is its behavior conservative or overly aggressive?\n\n### Role in the overall argument\nCRM is the first evidence block:\n\n> **It shows that the method is useful.**\n\nIn other words, it addresses **effectiveness**.\n\n---\n\n## 2. Mayo External Paired Evaluation\n\n### What it measures\nThe Mayo experiment evaluates the model on **real external paired low-dose / full-dose CT**.\n\nHere the comparison is made directly between:\n- real Quarter Dose CT,\n- the restored output,\n- and real Full Dose CT as reference.\n\nThe key question is not whether the model performs well on data generated by the same synthetic process used in training, but rather:\n\n> **Does training on synthetic physics-based pairs transfer to real low-dose clinical data?**\n\n### Why it matters\nThis is the strongest evidence for **sim-to-real transfer** in the notebook.  \nThe training supervision was not built from true paired low-dose / high-dose scans of the same patient. It was created through physics-grounded degradation synthesis.  \nTherefore, the critical test is whether this synthetic supervision teaches the model something that remains valid on real low-dose CT.\n\n### Core questions it answers\n- Is the synthetic supervision externally credible?\n- Does the model improve over the original real low-dose input?\n- Does it do so conservatively, rather than by aggressively rewriting the image?\n\n### Role in the overall argument\nMayo is the second major evidence block:\n\n> **It shows that the method transfers to the real world.**\n\nIn other words, it addresses **transferability**.\n\n---\n\n## 3. TotalSegmentator Anatomy Analysis\n\n### What it measures\nThe TotalSegmentator analysis is not included to win a segmentation benchmark.  \nIts purpose is to answer a more basic question:\n\n> **Where does the model make changes?**\n\nSpecifically, it examines whether the restoration edits are:\n- concentrated in anatomically relevant regions,\n- associated with structures and boundaries,\n- or merely distributed as uniform smoothing across the whole image.\n\n### Why it matters\nA major risk in medical image restoration is that a model may improve certain metrics by making broad, indiscriminate changes that are not anatomically meaningful.  \nSuch behavior can look visually attractive while still being scientifically or clinically untrustworthy.\n\nTotalSegmentator provides an **anatomy-anchored view** of the edits:\n- if changes are nearly absent, the model may be very conservative;\n- if changes are measurable and concentrated around structures, that suggests the edits are not random;\n- if segmentation agreement improves relative to Full Dose CT, then the restoration is likely doing something structurally meaningful rather than only smoothing pixels.\n\n### Core questions it answers\n- Are the edits anatomically grounded?\n- Is the model performing structured restoration instead of global smoothing?\n- Does the restored output become more anatomically consistent with the reference?\n\n### Role in the overall argument\nTotalSegmentator is the third evidence block:\n\n> **It shows that the model’s edits are structured, interpretable, and therefore more trustworthy.**\n\nIn other words, it addresses **interpretability**.\n\n---\n\n## 4. Why these tests belong together\n\nThese downstream tests are not redundant.  \nThey form a progressive evidence chain:\n\n### CRM\nshows that the model has **functional value**.\n\n### Mayo\nshows that this value **transfers to real external low-dose CT**.\n\n### TotalSegmentator\nshows that the edits are **anatomically structured rather than arbitrary**.\n\nTogether, they support the broader methodological claim:\n\n> **Physics-grounded synthetic supervision is not merely a way to generate images that resemble low-quality CT. It is a way to train restoration models that are useful, externally transferable, and structurally interpretable even when true paired data are missing.**\n\n---\n\n## How to read this section\n\n- If you care most about **clinical usefulness**, focus on **CRM**.\n- If you care most about **sim-to-real transfer**, focus on **Mayo**.\n- If you care most about **where the model edits and why those edits should be trusted**, focus on **TotalSegmentator**.\n\nAll three sections support the same central idea:\n\n> **The contribution of this project is not only a restoration model, but a method for training medical restoration systems under paired-data scarcity, together with downstream tests that evaluate whether that method deserves to be trusted.**\n","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# Imports & 环境\n# ============================================================\nimport os, sys, gc, math, time, random, hashlib, inspect, shutil, subprocess\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport pydicom\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom contextlib import nullcontext\n\ntry:\n    from IPython.display import display\nexcept Exception:\n    display = print\n\ntry:\n    cv2.setNumThreads(0)\nexcept Exception:\n    pass\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nUSE_AMP = (device.type == \"cuda\")\nAMP_CTX = (lambda: torch.amp.autocast(\"cuda\")) if USE_AMP else (lambda: nullcontext())\n\nprint(\"device:\", device)\nprint(\"USE_AMP:\", USE_AMP)\n# ============================================================\n#  Config\n# ============================================================\nfrom pathlib import Path\nimport os\n\nRSNA_DATA_ROOT = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/series\"\nMETA_CSV = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train.csv\"\nTRAIN_LOCALIZERS_CSV = \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train_localizers.csv\"  # optional\n\n# ------------------------------------------------------------\n# Final model checkpoint: use the best Stage-B scientific model\n# ------------------------------------------------------------\nCKPT_PATH = \"/kaggle/input/datasets/linguoyuemma/deblur-stageb-best-pt/deblur_stageB_best.pt\"\n\n# Optional: split files from the scientific training pipeline\nTRAIN_UIDS_CSV  = \"/kaggle/input/datasets/linguoyuemma/refine-uids/train_uids.csv\"\nREFINE_UIDS_CSV = \"/kaggle/input/datasets/linguoyuemma/refine-uids/refine_uids.csv\"\nVAL_UIDS_CSV    = \"/kaggle/input/datasets/linguoyuemma/refine-uids/val_uids.csv\"\n\n# Clinical Judge (9th place flayer)\nPREDICTION_PY = \"/kaggle/input/datasets/mingzeli2009/rsna-prediction/prediction.py\"\nMODEL_BASE = \"/kaggle/input/models/tom99763/9th-place-models-rsna-iad/pytorch/default/1\"\nFLAYER_DIR = f\"{MODEL_BASE}/flayer/outputs_heatmap_aux_v1_acc2\"\n\n# ------------------------------------------------------------\n# Output directory\n# ------------------------------------------------------------\nOUTDIR = Path(\"/kaggle/working/vultimate_scientific_eval\")\nOUTDIR.mkdir(parents=True, exist_ok=True)\n\n# ------------------------------------------------------------\n# Data\n# ------------------------------------------------------------\nTARGET_D, TARGET_H, TARGET_W = 64, 448, 448\nTARGET_SHAPE = (TARGET_D, TARGET_H, TARGET_W)\n\nHU_MIN, HU_MAX = -1024.0, 3072.0\nHU_RANGE = HU_MAX - HU_MIN\nKEEP_MODALITIES = {\"CT\", \"CTA\"}\n\n# ------------------------------------------------------------\n# Degradation settings for evaluation\n# Keep these aligned with the training-time degradation model\n# ------------------------------------------------------------\nEVAL_T = 8.0\nEVAL_DOSE = \"quarter\"\nDIFFUSION_ALPHA = 0.20\nLAM = DIFFUSION_ALPHA          # backward compatibility for older cells\nBLUR_T_MAX = 8.0\nRESTORE_BATCH = 16\n\n# ------------------------------------------------------------\n# Evaluation scale\n# ------------------------------------------------------------\nSEED_CASES = 2026\n\nN_PILOT_COMPARE = 10      # quick sanity check\nN_OOD_EVAL = 50           # main CRM size; can later increase to 100/200\nN_MC_CASES = 100          # Monte Carlo stability test\nMC_SEEDS = [10, 42, 23, 55, 83, 9999, 7, 11, 19, 29]\n\n# ------------------------------------------------------------\n# TotalSegmentator\n# Set this to the actual number you want to run\n# ------------------------------------------------------------\nN_TOTALSEG_CASES = 40\nTOTALSEG_TASK = \"total\"\n\n\n# ------------------------------------------------------------\n# Diagnostics\n# ------------------------------------------------------------\nprint(\"CKPT_PATH exists:\", os.path.exists(CKPT_PATH))\nprint(\"TRAIN_UIDS_CSV exists:\", os.path.exists(TRAIN_UIDS_CSV))\nprint(\"REFINE_UIDS_CSV exists:\", os.path.exists(REFINE_UIDS_CSV))\nprint(\"VAL_UIDS_CSV exists:\", os.path.exists(VAL_UIDS_CSV))\nprint(\"RSNA_DATA_ROOT exists:\", os.path.exists(RSNA_DATA_ROOT))\nprint(\"META_CSV exists:\", os.path.exists(META_CSV))\nprint(\"OUTDIR:\", OUTDIR)\n# ============================================================\n# Load NEW model checkpoint + Clinical Judge\n# ============================================================\nassert os.path.exists(CKPT_PATH), f\"Checkpoint not found: {CKPT_PATH}\"\n\nckpt = torch.load(CKPT_PATH, map_location=\"cpu\")\nstate_dict = ckpt[\"model\"] if isinstance(ckpt, dict) and \"model\" in ckpt else ckpt\n\n# --- NEW V-Ultimate (FiLM + authority map) ---\nmodel_25d = DeblurUNet25D_Ultimate(\n    in_ch=4,\n    out_ch=1,\n    base=32,\n    res_min=0.02,\n    res_max=0.15,\n    meta_dim=4,\n    film_hidden=64,\n    authority_bias_init=2.0,\n).to(device)\n\nmodel_25d.load_state_dict(state_dict, strict=True)\nmodel_25d.eval()\n\nprint(\"✅ NEW V-Ultimate model loaded.\")\nprint(\"has auth_head?\", hasattr(model_25d, \"auth_head\"))\nprint(\"has film_e1?\", hasattr(model_25d, \"film_e1\"))\nprint(\"train_uids in ckpt:\", len(ckpt.get(\"train_uids\", [])) if isinstance(ckpt, dict) else \"NA\")\nprint(\"best_val_psnr:\", ckpt.get(\"best_val_psnr\", \"NA\") if isinstance(ckpt, dict) else \"NA\")\n\n# ---- Load clinical judge ----\nimport importlib.util\nassert os.path.exists(PREDICTION_PY), f\"prediction.py not found: {PREDICTION_PY}\"\n\nspec = importlib.util.spec_from_file_location(\"prediction\", PREDICTION_PY)\npred_mod = importlib.util.module_from_spec(spec)\nsys.modules[\"prediction\"] = pred_mod\nspec.loader.exec_module(pred_mod)\n\nclassifier = pred_mod.FlayerClassifier(flayer_dir=FLAYER_DIR)\nclassifier.load()\n\n@torch.no_grad()\ndef aneurysm_predict(volume_uint8):\n    return float(classifier.predict(volume_uint8)[\"aneurysm_prob\"])\n\nprint(\"✅ Clinical Judge loaded.\")\n\n# ============================================================\n# Minimal bridge + CRM utilities\n# Run this BEFORE the CRM cell\n# ============================================================\n\nimport os\nimport hashlib\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\n\n# ------------------------------------------------------------\n# 1) Output dirs\n# ------------------------------------------------------------\nCRM_OUTDIR = OUTDIR / \"clinical_rescue_matrix\"\nCRM_OUTDIR.mkdir(parents=True, exist_ok=True)\n\nMAYO_OUTDIR = OUTDIR / \"mayo_external_eval\"\nMAYO_OUTDIR.mkdir(parents=True, exist_ok=True)\n\nTOTALSEG_OUTDIR = OUTDIR / \"totalseg\"\nTOTALSEG_OUTDIR.mkdir(parents=True, exist_ok=True)\n\n# ------------------------------------------------------------\n# 2) Load UID splits -> dev_uid_set\n# ------------------------------------------------------------\ndef read_uid_csv(path):\n    if path is None or not os.path.exists(path):\n        return []\n    df = pd.read_csv(path)\n    if df.shape[1] == 0:\n        return []\n    return df.iloc[:, 0].astype(str).tolist()\n\ntrain_uids = read_uid_csv(TRAIN_UIDS_CSV)\nrefine_uids = read_uid_csv(REFINE_UIDS_CSV)\nval_uids = read_uid_csv(VAL_UIDS_CSV)\n\ntrain_uid_set = set(train_uids)\nrefine_uid_set = set(refine_uids)\nval_uid_set = set(val_uids)\ndev_uid_set = train_uid_set | refine_uid_set | val_uid_set\n\nprint(\"train_uid_set :\", len(train_uid_set))\nprint(\"refine_uid_set:\", len(refine_uid_set))\nprint(\"val_uid_set   :\", len(val_uid_set))\nprint(\"dev_uid_set   :\", len(dev_uid_set))\n\n# ------------------------------------------------------------\n# 3) Small helpers\n# ------------------------------------------------------------\ndef stable_uid_seed(uid: str) -> int:\n    s = str(uid).encode(\"utf-8\")\n    return int(hashlib.md5(s).hexdigest()[:8], 16)\n\ndef uid_tail4(uid: str) -> str:\n    return str(uid)[-4:]\n\ndef _clip01(x):\n    return np.clip(np.asarray(x, dtype=np.float32), 0.0, 1.0).astype(np.float32)\n\ndef vol01_to_flayer_uint8(vol01):\n    vol01 = _clip01(vol01)\n    return (vol01 * 255.0).astype(np.uint8)\n\n# ------------------------------------------------------------\n# 4) CRM metrics\n# ------------------------------------------------------------\ndef calc_abs_gain(p_gt, p_deg, p_rec):\n    \"\"\"\n    Raw signed change in model score after restoration.\n    Positive means p_rec > p_deg, negative means p_rec < p_deg.\n    \"\"\"\n    return float(p_rec - p_deg)\n\ndef calc_target_gain(p_gt, p_deg, p_rec):\n    \"\"\"\n    Improvement toward the ground-truth score.\n    Positive means restoration moved closer to p_gt.\n    \"\"\"\n    return float(abs(p_deg - p_gt) - abs(p_rec - p_gt))\n\ndef is_iatrogenic(p_gt, p_deg, p_rec, eps=1e-8):\n    \"\"\"\n    Restoration is harmful if it moves farther away from p_gt than the degraded input.\n    \"\"\"\n    return bool((abs(p_rec - p_gt) - abs(p_deg - p_gt)) > eps)\n\nprint(\"✅ CRM bridge utilities ready.\")\nprint(\"CRM_OUTDIR:\", CRM_OUTDIR)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T00:43:57.543395Z","iopub.execute_input":"2026-03-15T00:43:57.544096Z","iopub.status.idle":"2026-03-15T00:44:03.655934Z","shell.execute_reply.started":"2026-03-15T00:43:57.544059Z","shell.execute_reply":"2026-03-15T00:44:03.655151Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# Clinical Rescue Matrix — streamlined evaluation-only version\n# STRICT held-out RSNA CT/CTA OOD evaluation\n# Depends on:\n#   - model_25d\n#   - RSNA_DATA_ROOT, META_CSV, CRM_OUTDIR\n#   - dev_uid_set (from bridge cell)\n#   - aneurysm_predict, load_series_volume, stable_uid_seed, uid_tail4\n#   - _clip01, calc_abs_gain, calc_target_gain, is_iatrogenic\n#   - vol01_to_flayer_uint8\n# ============================================================\n\nimport os, gc, time, math, random\nfrom pathlib import Path\nfrom contextlib import nullcontext\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport torch\n\ntry:\n    from IPython.display import display\nexcept Exception:\n    display = print\n\ntry:\n    cv2.setNumThreads(0)\nexcept Exception:\n    pass\n\n# ------------------------------------------------------------\n# 0) Checks / config\n# ------------------------------------------------------------\nrequired_globals = [\n    \"model_25d\", \"device\", \"AMP_CTX\",\n    \"RSNA_DATA_ROOT\", \"META_CSV\", \"CRM_OUTDIR\",\n    \"dev_uid_set\",\n    \"load_series_volume\", \"stable_uid_seed\", \"uid_tail4\",\n    \"_clip01\", \"calc_abs_gain\", \"calc_target_gain\", \"is_iatrogenic\",\n    \"vol01_to_flayer_uint8\", \"aneurysm_predict\",\n]\nmissing = [k for k in required_globals if k not in globals()]\nif missing:\n    raise RuntimeError(f\"Missing required globals/functions: {missing}\")\n\nOUTDIR = Path(CRM_OUTDIR)\nOUTDIR.mkdir(parents=True, exist_ok=True)\n\nTARGET_D = int(globals().get(\"TARGET_D\", 64))\nTARGET_H = int(globals().get(\"TARGET_H\", 448))\nTARGET_W = int(globals().get(\"TARGET_W\", 448))\nTARGET_SHAPE = (TARGET_D, TARGET_H, TARGET_W)\n\nHU_MIN = float(globals().get(\"HU_MIN\", -1024.0))\nHU_MAX = float(globals().get(\"HU_MAX\", 3072.0))\nHU_RANGE = HU_MAX - HU_MIN\n\nDIFFUSION_ALPHA = float(globals().get(\"DIFFUSION_ALPHA\", globals().get(\"LAM\", 0.20)))\nBLUR_T_MAX = float(globals().get(\"BLUR_T_MAX\", 8.0))\nRESTORE_BATCH = int(globals().get(\"RESTORE_BATCH\", 16))\n\nEVAL_T = float(globals().get(\"EVAL_T\", 8.0))\nEVAL_DOSE = str(globals().get(\"EVAL_DOSE\", \"quarter\"))\nN_CRM_CASES = int(globals().get(\"N_OOD_EVAL\", 50))\nSEED_CRM = 2026\n\nGAIN_POS_TH = 0.005\nGAIN_NEG_TH = -0.005\n\nGAUSS_SIGMA = 0.8\nMEDIAN_KSIZE = 3\nBILATERAL_D = 7\nBILATERAL_SIGMACOLOR = 35\nBILATERAL_SIGMASPACE = 35\nUNSHARP_SIGMA = 1.0\nUNSHARP_AMOUNT = 0.8\nUSE_NLM = False\nNLM_H = 7\n\nRAW_CSV = OUTDIR / f\"crm_raw_N{N_CRM_CASES}.csv\"\nSUM_CSV = OUTDIR / f\"crm_summary_N{N_CRM_CASES}.csv\"\nPAIR_CSV = OUTDIR / f\"crm_paired_vs_vultimate_N{N_CRM_CASES}.csv\"\nBUCKET_CSV = OUTDIR / f\"crm_bucket_breakdown_N{N_CRM_CASES}.csv\"\nFIXED_UID_CSV = OUTDIR / f\"crm_eval_uids_fixed_N{N_CRM_CASES}.csv\"\nUSED_UID_CSV = OUTDIR / f\"crm_used_eval_uids_N{N_CRM_CASES}.csv\"\n\nMETHODS = [\n    (\"Degraded\", \"degraded\"),\n    (\"Gaussian\", \"gaussian\"),\n    (\"Median3\", \"median3\"),\n    (\"Bilateral\", \"bilateral\"),\n    (\"Unsharp\", \"unsharp\"),\n]\nif USE_NLM:\n    METHODS.append((\"NLM\", \"nlm\"))\nMETHODS.append((\"V-Ultimate\", \"vultimate\"))\n\nprint(\"Methods:\", [m[0] for m in METHODS])\n\n# ------------------------------------------------------------\n# 1) Helpers\n# ------------------------------------------------------------\ndef psnr01_full(a, b):\n    a = np.asarray(a, dtype=np.float32)\n    b = np.asarray(b, dtype=np.float32)\n    mse = float(np.mean((a - b) ** 2))\n    return 99.0 if mse <= 0 else 10.0 * math.log10(1.0 / mse)\n\ndef crm_outcome_label(p_gt, p_deg, p_rec, abs_gain):\n    if is_iatrogenic(p_gt, p_deg, p_rec):\n        return \"⚠️ Iatrogenic\"\n    if abs_gain > GAIN_POS_TH:\n        return \"✅ Successful Rescue\"\n    if abs_gain < GAIN_NEG_TH:\n        return \"⏬ Minor Deviation\"\n    return \"➖ Algorithmic Humility\"\n\ndef crm_bucket(p_gt, p_deg, tau=0.03):\n    if p_gt >= 0.5 and p_deg < p_gt - tau:\n        return \"harmed_positive\"\n    if p_gt >= 0.5 and p_deg > p_gt + tau:\n        return \"overcall_positive\"\n    return \"neutral_other\"\n\ndef hu01_to_hu(vol01):\n    return np.asarray(vol01, dtype=np.float32) * HU_RANGE + HU_MIN\n\ndef hu_to_01(hu):\n    return np.clip((hu - HU_MIN) / HU_RANGE, 0.0, 1.0).astype(np.float32)\n\ndef window_hu_to_uint8(hu, center=40.0, width=400.0):\n    x = np.clip((hu - (center - width / 2.0)) / (width + 1e-6), 0.0, 1.0)\n    return (x * 255.0).astype(np.uint8)\n\ndef degrade_volume_fixed_uid(vol01, t, uid, dose_mode=\"quarter\"):\n    local_seed = stable_uid_seed(uid)\n    py_state, np_state = random.getstate(), np.random.get_state()\n    random.seed(local_seed)\n    np.random.seed(local_seed % (2**32 - 1))\n\n    sigma = math.sqrt(max(1e-8, 2.0 * DIFFUSION_ALPHA * float(t)))\n    out = np.empty_like(vol01, dtype=np.float32)\n\n    for z in range(vol01.shape[0]):\n        x = cv2.GaussianBlur(\n            vol01[z].astype(np.float32),\n            (0, 0),\n            sigmaX=sigma,\n            sigmaY=sigma,\n            borderType=cv2.BORDER_REPLICATE,\n        )\n        if dose_mode != \"clean\":\n            peak = random.uniform(3000.0, 6000.0)\n            sigma_e = random.uniform(0.01, 0.02)\n            x = np.random.poisson(np.clip(x * peak, 0, None)).astype(np.float32) / peak\n            x = x + np.random.randn(*x.shape).astype(np.float32) * sigma_e\n        out[z] = np.clip(x, 0.0, 1.0).astype(np.float32)\n\n    random.setstate(py_state)\n    np.random.set_state(np_state)\n    return out\n\ndef trad_identity(vol):\n    return np.asarray(vol, dtype=np.float32).copy()\n\ndef trad_gaussian(vol, sigma=GAUSS_SIGMA):\n    out = np.empty_like(vol, dtype=np.float32)\n    for z in range(vol.shape[0]):\n        out[z] = _clip01(cv2.GaussianBlur(\n            vol[z], (0, 0), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REPLICATE\n        ))\n    return out\n\ndef trad_median(vol, ksize=MEDIAN_KSIZE):\n    out = np.empty_like(vol, dtype=np.float32)\n    for z in range(vol.shape[0]):\n        x8 = (np.clip(vol[z], 0, 1) * 255).astype(np.uint8)\n        out[z] = cv2.medianBlur(x8, ksize).astype(np.float32) / 255.0\n    return out\n\ndef trad_bilateral(vol, d=BILATERAL_D, sigmaColor=BILATERAL_SIGMACOLOR, sigmaSpace=BILATERAL_SIGMASPACE):\n    out = np.empty_like(vol, dtype=np.float32)\n    for z in range(vol.shape[0]):\n        x8 = (np.clip(vol[z], 0, 1) * 255).astype(np.uint8)\n        out[z] = cv2.bilateralFilter(\n            x8, d=d, sigmaColor=sigmaColor, sigmaSpace=sigmaSpace\n        ).astype(np.float32) / 255.0\n    return out\n\ndef trad_unsharp(vol, sigma=UNSHARP_SIGMA, amount=UNSHARP_AMOUNT):\n    out = np.empty_like(vol, dtype=np.float32)\n    for z in range(vol.shape[0]):\n        x = vol[z].astype(np.float32)\n        blur = cv2.GaussianBlur(x, (0, 0), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REPLICATE)\n        out[z] = np.clip(x + amount * (x - blur), 0.0, 1.0)\n    return out\n\ndef trad_nlm(vol, h=NLM_H, center=40.0, width=400.0):\n    vol_hu = hu01_to_hu(vol)\n    out = np.empty_like(vol, dtype=np.float32)\n    for z in range(vol.shape[0]):\n        u8 = window_hu_to_uint8(vol_hu[z], center=center, width=width)\n        den = cv2.fastNlMeansDenoising(u8, None, h=h, templateWindowSize=7, searchWindowSize=21)\n        den01 = den.astype(np.float32) / 255.0\n        den_hu = den01 * width + (center - width / 2.0)\n        out[z] = hu_to_01(den_hu)\n    return out\n\nTRAD_METHODS = {\n    \"degraded\": trad_identity,\n    \"gaussian\": trad_gaussian,\n    \"median3\": trad_median,\n    \"bilateral\": trad_bilateral,\n    \"unsharp\": trad_unsharp,\n    \"nlm\": trad_nlm,\n}\n\n@torch.no_grad()\ndef run_vultimate(vol_deg01, t=EVAL_T, restore_batch=RESTORE_BATCH):\n    vol_deg01 = np.asarray(vol_deg01, dtype=np.float32)\n    D = vol_deg01.shape[0]\n    out = vol_deg01.copy()\n\n    t_norm = np.float32(0.0 if t <= 0 else float(t) / float(BLUR_T_MAX))\n    meta_row = np.array([t_norm, 0.0, 0.0, 0.0], dtype=np.float32)\n\n    for s in range(0, D, restore_batch):\n        zs = list(range(s, min(D, s + restore_batch)))\n        inp_batch = []\n        for z in zs:\n            bp = vol_deg01[max(0, z - 1)]\n            bc = vol_deg01[z]\n            bn = vol_deg01[min(D - 1, z + 1)]\n            inp_batch.append(np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32))\n\n        inp_t = torch.from_numpy(np.stack(inp_batch, axis=0)).to(device, non_blocking=True)\n        meta_t = torch.from_numpy(np.repeat(meta_row[None, :], len(zs), axis=0)).to(device, non_blocking=True)\n\n        with AMP_CTX():\n            pred_obj = model_25d(inp_t, meta=meta_t)\n            pred = pred_obj[0] if isinstance(pred_obj, (tuple, list)) else pred_obj\n            pred = pred.float().cpu().numpy()[:, 0]\n\n        for k, z in enumerate(zs):\n            out[z] = _clip01(pred[k])\n\n    return out\n\ndef apply_method(method_key, deg):\n    if method_key == \"vultimate\":\n        return run_vultimate(deg, t=EVAL_T, restore_batch=RESTORE_BATCH)\n    if method_key in TRAD_METHODS:\n        return TRAD_METHODS[method_key](deg)\n    raise ValueError(f\"Unknown method_key: {method_key}\")\n\n# ------------------------------------------------------------\n# 2) Build strict held-out OOD pool\n# ------------------------------------------------------------\nprint(\"=== Clinical Rescue Matrix | STRICT RSNA OOD CT/CTA ===\")\n\nmeta = pd.read_csv(META_CSV)\nct_uids = set(\n    meta.loc[meta[\"Modality\"].astype(str).isin({\"CT\", \"CTA\"}), \"SeriesInstanceUID\"].astype(str).tolist()\n)\n\nall_series_dirs = [\n    u for u in os.listdir(RSNA_DATA_ROOT)\n    if os.path.isdir(os.path.join(RSNA_DATA_ROOT, u))\n]\n\nood_pool = [u for u in all_series_dirs if (u in ct_uids) and (u not in dev_uid_set)]\n\nprint(\"\\n=== Exclusion summary ===\")\nprint(\"dev_uid_set :\", len(dev_uid_set))\nprint(\"ood_pool    :\", len(ood_pool))\nprint(\"overlap     :\", len(set(ood_pool) & dev_uid_set), \"(should be 0)\")\n\nif FIXED_UID_CSV.exists():\n    crm_eval_uids = pd.read_csv(FIXED_UID_CSV)[\"uid_full\"].astype(str).tolist()\n    crm_eval_uids = [u for u in crm_eval_uids if u in ood_pool]\n    print(f\"Loaded fixed eval UID list: {len(crm_eval_uids)}\")\nelse:\n    random.Random(SEED_CRM).shuffle(ood_pool)\n    crm_eval_uids = ood_pool[:N_CRM_CASES]\n    pd.DataFrame({\"uid_full\": crm_eval_uids}).to_csv(FIXED_UID_CSV, index=False)\n    print(f\"Saved fixed eval UID list: {len(crm_eval_uids)}\")\n\n# ------------------------------------------------------------\n# 3) Main loop\n# ------------------------------------------------------------\nrows = []\nused_eval_uids = []\nt0_all = time.time()\n\nfor i, uid in enumerate(crm_eval_uids, 1):\n    case_t0 = time.time()\n\n    try:\n        try:\n            vol = load_series_volume(uid, RSNA_DATA_ROOT, TARGET_SHAPE)\n        except TypeError:\n            try:\n                vol = load_series_volume(uid, RSNA_DATA_ROOT)\n            except TypeError:\n                vol = load_series_volume(uid)\n    except Exception:\n        vol = None\n\n    if vol is None:\n        continue\n\n    gt = vol.astype(np.float32)\n    deg = degrade_volume_fixed_uid(gt, EVAL_T, uid, dose_mode=EVAL_DOSE)\n\n    p_gt = float(aneurysm_predict(vol01_to_flayer_uint8(gt)))\n    p_deg = float(aneurysm_predict(vol01_to_flayer_uint8(deg)))\n    bucket = crm_bucket(p_gt, p_deg, tau=0.03)\n    uid4 = uid_tail4(uid)\n\n    print(f\"\\n[{len(used_eval_uids)+1:03d}/{len(crm_eval_uids)}] UID:{uid4} | GT:{p_gt:.4f} -> Deg:{p_deg:.4f} | bucket={bucket}\")\n\n    case_ok = False\n    for m_name, m_key in METHODS:\n        try:\n            rec = apply_method(m_key, deg)\n            p_rec = float(aneurysm_predict(vol01_to_flayer_uint8(rec)))\n            abs_gain = float(calc_abs_gain(p_gt, p_deg, p_rec))\n            tgt_gain = float(calc_target_gain(p_gt, p_deg, p_rec))\n            iatro = int(is_iatrogenic(p_gt, p_deg, p_rec))\n            psnr_val = float(psnr01_full(rec, gt))\n            outcome = crm_outcome_label(p_gt, p_deg, p_rec, abs_gain)\n\n            rows.append({\n                \"uid_full\": uid,\n                \"uid4\": uid4,\n                \"bucket\": bucket,\n                \"method\": m_name,\n                \"method_key\": m_key,\n                \"p_gt\": p_gt,\n                \"p_deg\": p_deg,\n                \"p_rec\": p_rec,\n                \"abs_gain\": abs_gain,\n                \"target_gain\": tgt_gain,\n                \"iatrogenic\": iatro,\n                \"outcome\": outcome,\n                \"psnr_db\": psnr_val,\n                \"eval_t\": EVAL_T,\n                \"eval_dose\": EVAL_DOSE,\n            })\n\n            tag = (\n                \"✅ rescue\" if tgt_gain > GAIN_POS_TH\n                else \"⚠️ negative\" if tgt_gain < GAIN_NEG_TH\n                else \"➖ identity\"\n            )\n            print(\n                f\"  ├─ {m_name:<10s} | Rec:{p_rec:.4f} | \"\n                f\"TGain:{tgt_gain:+.4f} | AGain:{abs_gain:+.4f} | {outcome} | {tag}\"\n            )\n            case_ok = True\n\n            if m_key != \"degraded\":\n                del rec\n\n        except Exception as e:\n            rows.append({\n                \"uid_full\": uid,\n                \"uid4\": uid4,\n                \"bucket\": bucket,\n                \"method\": m_name,\n                \"method_key\": m_key,\n                \"error\": repr(e),\n            })\n            print(f\"  ├─ {m_name:<10s} | ERROR: {repr(e)}\")\n\n    if case_ok:\n        used_eval_uids.append(uid)\n\n    print(f\"  -> case done in {time.time() - case_t0:.1f}s\")\n\n    del gt, deg, vol\n    gc.collect()\n    if torch.cuda.is_available():\n        torch.cuda.empty_cache()\n\npd.DataFrame({\"uid_full\": used_eval_uids}).to_csv(USED_UID_CSV, index=False)\n\nprint(f\"\\nAll CRM cases done. elapsed={(time.time() - t0_all)/60:.1f} min\")\nprint(f\"Valid cases: {len(used_eval_uids)}\")\n\n# ------------------------------------------------------------\n# 4) Save + summary\n# ------------------------------------------------------------\ndf = pd.DataFrame(rows)\ndf.to_csv(RAW_CSV, index=False)\nif len(df) == 0:\n    raise RuntimeError(\"No CRM results generated.\")\n\ndf_ok = df.dropna(subset=[\"target_gain\", \"abs_gain\", \"psnr_db\"]).copy()\n\nsummary_rows = []\nfor m_name, _ in METHODS:\n    sub = df_ok[df_ok[\"method\"] == m_name].copy()\n    if len(sub) == 0:\n        continue\n\n    oc = sub[\"outcome\"].value_counts()\n    summary_rows.append({\n        \"Method\": m_name,\n        \"N\": int(len(sub)),\n        \"Mean Target Gain\": round(float(sub[\"target_gain\"].mean()), 4),\n        \"Target Win %\": round(float((sub[\"target_gain\"] > GAIN_POS_TH).mean() * 100), 1),\n        \"Target Neg %\": round(float((sub[\"target_gain\"] < GAIN_NEG_TH).mean() * 100), 1),\n        \"Mean Abs Gain\": round(float(sub[\"abs_gain\"].mean()), 4),\n        \"Iatrogenic %\": round(float(sub[\"iatrogenic\"].mean() * 100), 1),\n        \"PSNR (dB)\": round(float(sub[\"psnr_db\"].mean()), 2),\n        \"CRM_Rescue\": f\"{int(oc.get('✅ Successful Rescue', 0))} ({100.0 * oc.get('✅ Successful Rescue', 0) / len(sub):.1f}%)\",\n        \"CRM_Humble\": f\"{int(oc.get('➖ Algorithmic Humility', 0))} ({100.0 * oc.get('➖ Algorithmic Humility', 0) / len(sub):.1f}%)\",\n        \"CRM_Deviate\": f\"{int(oc.get('⏬ Minor Deviation', 0))} ({100.0 * oc.get('⏬ Minor Deviation', 0) / len(sub):.1f}%)\",\n        \"CRM_Iatro\": f\"{int(oc.get('⚠️ Iatrogenic', 0))} ({100.0 * oc.get('⚠️ Iatrogenic', 0) / len(sub):.1f}%)\",\n    })\n\ndf_sum = pd.DataFrame(summary_rows)\ndf_sum.to_csv(SUM_CSV, index=False)\n\nprint(\"\\n=== Clinical Rescue Matrix — Summary ===\")\ndisplay(df_sum)\n\n# bucket breakdown\nbucket_rows = []\nfor bucket in sorted(df_ok[\"bucket\"].dropna().unique()):\n    for m_name, _ in METHODS:\n        sub = df_ok[(df_ok[\"bucket\"] == bucket) & (df_ok[\"method\"] == m_name)].copy()\n        if len(sub) == 0:\n            continue\n        bucket_rows.append({\n            \"bucket\": bucket,\n            \"method\": m_name,\n            \"N\": int(len(sub)),\n            \"mean_target_gain\": round(float(sub[\"target_gain\"].mean()), 4),\n            \"mean_abs_gain\": round(float(sub[\"abs_gain\"].mean()), 4),\n            \"iatrogenic_%\": round(float(sub[\"iatrogenic\"].mean() * 100), 1),\n            \"psnr_db\": round(float(sub[\"psnr_db\"].mean()), 2),\n        })\n\ndf_bucket = pd.DataFrame(bucket_rows)\ndf_bucket.to_csv(BUCKET_CSV, index=False)\n\nprint(\"\\n=== Bucket breakdown ===\")\ndisplay(df_bucket)\n\n# paired vs V-Ultimate\nif \"V-Ultimate\" in set(df_ok[\"method\"].unique()):\n    base = df_ok[df_ok[\"method\"] == \"V-Ultimate\"][[\n        \"uid_full\", \"target_gain\", \"abs_gain\", \"iatrogenic\", \"psnr_db\"\n    ]].rename(columns={\n        \"target_gain\": \"tgain_base\",\n        \"abs_gain\": \"again_base\",\n        \"iatrogenic\": \"iatro_base\",\n        \"psnr_db\": \"psnr_base\",\n    })\n\n    paired_rows = []\n    for m_name, _ in METHODS:\n        if m_name == \"V-Ultimate\":\n            continue\n\n        cmp_df = df_ok[df_ok[\"method\"] == m_name][[\n            \"uid_full\", \"target_gain\", \"abs_gain\", \"iatrogenic\", \"psnr_db\"\n        ]].rename(columns={\n            \"target_gain\": \"tgain_cmp\",\n            \"abs_gain\": \"again_cmp\",\n            \"iatrogenic\": \"iatro_cmp\",\n            \"psnr_db\": \"psnr_cmp\",\n        })\n\n        merged = base.merge(cmp_df, on=\"uid_full\", how=\"inner\")\n        if len(merged) == 0:\n            continue\n\n        dt = merged[\"tgain_base\"] - merged[\"tgain_cmp\"]\n        paired_rows.append({\n            \"vs\": m_name,\n            \"N\": int(len(merged)),\n            \"V wins TGain %\": round(float((dt > 0.001).mean() * 100), 1),\n            \"V loses TGain %\": round(float((dt < -0.001).mean() * 100), 1),\n            \"ΔTGain\": round(float(dt.mean()), 4),\n            \"ΔPSNR\": round(float((merged[\"psnr_base\"] - merged[\"psnr_cmp\"]).mean()), 3),\n        })\n\n    df_pair = pd.DataFrame(paired_rows)\n    df_pair.to_csv(PAIR_CSV, index=False)\n\n    print(\"\\n=== Paired comparison (baseline: V-Ultimate) ===\")\n    display(df_pair)\n\nprint(\"\\nSaved:\")\nprint(\" raw      ->\", RAW_CSV)\nprint(\" summary  ->\", SUM_CSV)\nprint(\" bucket   ->\", BUCKET_CSV)\nprint(\" pair     ->\", PAIR_CSV if PAIR_CSV.exists() else \"(not generated)\")\nprint(\" fixedUID ->\", FIXED_UID_CSV)\nprint(\" usedUID  ->\", USED_UID_CSV)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T00:47:50.408733Z","iopub.execute_input":"2026-03-15T00:47:50.409039Z","iopub.status.idle":"2026-03-15T01:32:07.225123Z","shell.execute_reply.started":"2026-03-15T00:47:50.409007Z","shell.execute_reply":"2026-03-15T01:32:07.224552Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Results Interpretation — Clinical Rescue Matrix (Strict RSNA OOD)\n\n### 1. What question does this experiment actually answer?\n\nThe Clinical Rescue Matrix (CRM) is not designed to test which method produces the smoothest image. Instead, it asks:\n\n> When degradation suppresses clinically relevant signal, can a restoration method **recover the lost target signal** while minimizing harmful deviations (iatrogenic change)?\n\nFor that reason, the key metrics here are not PSNR alone, but:\n\n- **Mean Target Gain**: average clinically relevant recovery\n- **Target Win %**: fraction of cases with clear positive rescue\n- **Iatrogenic %**: fraction of cases with harmful deviation\n- **Bucket breakdown**: which kinds of cases each method actually helps\n\n---\n\n### 2. Overall result: V-Ultimate is the strongest overall method, but not a universal winner on every single case\n\nAcross 50 strict held-out OOD CT/CTA cases, V-Ultimate achieved:\n\n- **Mean Target Gain = 0.0324** (**highest** among all methods)\n- **Target Win % = 66.0%** (**highest**)\n- **Iatrogenic % = 32.0%** (**tied lowest** with Bilateral)\n- **PSNR = 42.74 dB** (**substantially higher** than all traditional baselines)\n- **Successful Rescue = 31/50 = 62.0%** (**highest**)\n\nThis means that a model trained only with synthetic supervision is already capable of producing the **strongest overall clinical proxy benefit** on strict held-out OOD data, while maintaining a relatively low harmful-modification rate.\n\nThe most important point is not that V-Ultimate wins every case, but that:\n\n> **At the population level, V-Ultimate gives the best combination of benefit, win rate, and safety.**\n\n---\n\n### 3. Gaussian is not weak, but it is not the end of the story\n\nGaussian is a strong baseline in this benchmark:\n\n- Mean Target Gain = **0.0125**\n- Target Win % = **50.0%**\n- Iatrogenic % = **44.0%**\n\nThis shows that simple smoothing can help under some degradation patterns. Therefore, these results should **not** be interpreted as:\n\n> “Traditional methods are useless, so deep learning is always necessary.”\n\nA more accurate interpretation is:\n\n> **Simple smoothing can recover signal in some cases, but V-Ultimate achieves higher overall benefit, higher win rate, and better safety, suggesting that it learns something more targeted than generic denoising.**\n\nIn other words, Gaussian serves as a strong control condition. It helps distinguish whether the model is doing:\n\n- ordinary smoothing, or\n- more targeted, controlled restoration\n\n---\n\n### 4. The most important finding: V-Ultimate is strongest where degradation truly harms the target\n\nThe most important CRM bucket is **`harmed_positive`**, meaning:\n\n> the original high-quality scan contains a strong target signal, but degradation clearly suppresses it.\n\nThis bucket is the most clinically meaningful because it represents cases in which degradation actually hides information that should be detected.\n\nIn this key bucket (**N = 36**), V-Ultimate achieved:\n\n- **mean_target_gain = 0.0485** (**highest**)\n- **mean_abs_gain = 0.0664** (**highest**)\n- **iatrogenic = 16.7%** (**lowest**)\n- **psnr = 42.56 dB** (**highest**)\n\nBy comparison:\n\n- Gaussian: target_gain = **0.0194**, iatrogenic = **36.1%**\n- Bilateral: target_gain = **0.0206**, iatrogenic = **22.2%**\n- Median3: target_gain = **0.0082**, iatrogenic = **38.9%**\n\nThis is the central result of the CRM experiment:\n\n> **When degradation truly damages clinically relevant signal, V-Ultimate is the best method at rescuing that signal, and it does so with the lowest harm rate.**\n\nThat is much more meaningful than simply saying the model “denoises better.”\n\n---\n\n### 5. This suggests the model is learning targeted recovery, not blind smoothing\n\nIf a method is merely smoothing, two patterns often appear:\n\n1. it may help in some clearly degraded cases;\n2. but it also tends to over-edit or over-smooth cases that do not need much intervention.\n\nV-Ultimate’s gains are concentrated most strongly in the **`harmed_positive`** bucket, which is exactly where recovery is most needed.\n\nThat matches the intended design of the method:\n\n> **The model is not trying to smooth every image uniformly; it is more effective when signal has actually been damaged and needs rescue.**\n\nThis makes V-Ultimate better described as a **target-oriented restoration model** rather than just a generic denoiser.\n\n---\n\n### 6. Important limitation: neutral and overcall cases remain challenging\n\nThe results also show clearly that V-Ultimate is **not perfect**.\n\nIn the **`neutral_other`** bucket (**N = 12**):\n\n- mean_target_gain = **-0.0123**\n- iatrogenic = **75.0%**\n\nThis suggests that:\n\n> when degradation does not strongly suppress the target signal, the model can still introduce unnecessary edits.\n\nIn the **`overcall_positive`** bucket (**N = 2**), the sample size is too small for strong conclusions, but the current result also indicates instability in cases where the task is to reduce inflated signal rather than recover suppressed signal.\n\nTherefore, the correct conclusion is **not**:\n\n> “V-Ultimate is best in every type of case.”\n\nInstead, it is:\n\n> **V-Ultimate has a very clear strength profile: it is most useful in the clinically important cases where degradation truly harms the target, while still leaving room for improvement in cases that may not require substantial correction.**\n\nThis is exactly why the later Monte Carlo, Mayo, and anatomy-aware analyses are important.\n\n---\n\n### 7. Pairwise comparison further supports that the advantage is systematic, not driven by a few lucky cases\n\nUsing V-Ultimate as the reference in paired comparisons:\n\n- vs Gaussian:  \n  - **V wins TGain % = 64.0%**\n  - ΔTGain = **+0.0199**\n\n- vs Median3:  \n  - **V wins TGain % = 74.0%**\n  - ΔTGain = **+0.0288**\n\n- vs Bilateral:  \n  - **V wins TGain % = 62.0%**\n  - ΔTGain = **+0.0209**\n\n- vs Unsharp:  \n  - **V wins TGain % = 82.0%**\n  - ΔTGain = **+0.0499**\n\nThis shows that V-Ultimate’s advantage is not driven by a handful of outliers. Instead, it wins in the majority of matched case-level comparisons, including against strong traditional baselines such as Gaussian and Bilateral.\n\n---\n\n### 8. Final interpretation of this section\n\nThe CRM results support not a simple “score-chasing” conclusion, but a more important methodological one:\n\n> **A controlled restoration model trained only with physics-based synthetic supervision can already produce real, systematic recovery benefits on strict held-out OOD clinical proxy tasks.**\n\nMore specifically:\n\n1. **It achieves the highest overall target gain and the highest repair win rate;**\n2. **Its strongest advantage appears in the clinically important cases where degradation actually suppresses signal;**\n3. **Its behavior cannot be explained by generic smoothing alone, because its gains are concentrated in rescue-relevant cases;**\n4. **At the same time, it still struggles in neutral or overcall settings, which motivates the later stability analysis and external real-data validation.**\n\nTherefore, the most accurate conclusion is not:\n\n> “V-Ultimate completely crushes all traditional baselines.”\n\nIt is:\n\n> **V-Ultimate demonstrates that physics-based synthetic supervision can train a controlled restoration model that meaningfully rescues clinically relevant signal under strict OOD conditions.**","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# 8) Stage C: Monte Carlo Noise Stability Test (NEW MODEL)\n# Compatible with: FiLM + authority map V-Ultimate\n# ============================================================\n\nimport os, sys, time, math, random, gc\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport torch\nfrom contextlib import nullcontext\n\ntry:\n    from IPython.display import display\nexcept Exception:\n    display = print\n\ntry:\n    cv2.setNumThreads(0)\nexcept Exception:\n    pass\n\n# -----------------------------\n# 0) Prerequisites check\n# -----------------------------\nrequired_any = {\n    \"model\": (\"model_25d\" in globals()),\n    \"judge\": (\"aneurysm_predict\" in globals()),\n    \"loader\": (\"load_series_volume\" in globals()),\n}\nif not all(required_any.values()):\n    missing = [k for k, v in required_any.items() if not v]\n    raise RuntimeError(\n        f\"Missing prerequisite objects/functions: {missing}\\n\"\n        f\"Please run: new model loading + Clinical Judge + DICOM loading cells first.\"\n    )\n\nMODEL_OBJ = globals().get(\"model_25d\", None)\nassert MODEL_OBJ is not None, \"model_25d not found\"\nassert hasattr(MODEL_OBJ, \"auth_head\"), \"Not the new model: missing auth_head\"\nassert hasattr(MODEL_OBJ, \"film_e1\"), \"Not the new model: missing FiLM module\"\n\n# -----------------------------\n# 1) Configuration\n# -----------------------------\nRSNA_DATA_ROOT = globals().get(\"RSNA_DATA_ROOT\", \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/series\")\nMETA_CSV       = globals().get(\"META_CSV\", \"/kaggle/input/competitions/rsna-intracranial-aneurysm-detection/train.csv\")\n\nTARGET_D = int(globals().get(\"TARGET_D\", 64))\nTARGET_H = int(globals().get(\"TARGET_H\", 448))\nTARGET_W = int(globals().get(\"TARGET_W\", 448))\nTARGET_SHAPE = (TARGET_D, TARGET_H, TARGET_W)\n\nHU_MIN   = float(globals().get(\"HU_MIN\", -1024.0))\nHU_MAX   = float(globals().get(\"HU_MAX\", 3072.0))\nHU_RANGE = float(globals().get(\"HU_RANGE\", HU_MAX - HU_MIN))\n\nDIFFUSION_ALPHA = float(globals().get(\"DIFFUSION_ALPHA\", globals().get(\"LAM\", 0.20)))\nBLUR_T_MAX = float(globals().get(\"BLUR_T_MAX\", 8.0))\n\nEVAL_T    = float(globals().get(\"EVAL_T\", 8.0))\nEVAL_DOSE = globals().get(\"EVAL_DOSE\", \"quarter\")\n\nN_MC_TO_RUN = int(globals().get(\"N_MC_CASES\", 100))\nMC_SEEDS    = list(globals().get(\"MC_SEEDS\", list(range(10))))\nRESTORE_BATCH = int(globals().get(\"RESTORE_BATCH\", 16))\n\nfrom pathlib import Path\nOUTDIR = str(Path(CRM_OUTDIR) / \"monte_carlo_stability\")\nos.makedirs(OUTDIR, exist_ok=True)\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nUSE_AMP = (device.type == \"cuda\")\nAMP_CTX = (lambda: torch.amp.autocast(\"cuda\")) if USE_AMP else (lambda: nullcontext())\n\n# Thresholds\nTARGET_POS_TH = 0.005\nTARGET_NEG_TH = -0.005\nABS_POS_TH = 0.005\nABS_NEG_TH = -0.005\n\n# -----------------------------\n# 2) Utility functions\n# -----------------------------\ndef mc_bucket(p_gt, p_deg, tau=0.03):\n    if p_gt >= 0.5 and p_deg < p_gt - tau:\n        return \"harmed_positive\"\n    if p_gt >= 0.5 and p_deg > p_gt + tau:\n        return \"overcall_positive\"\n    return \"neutral_other\"\n\ndef degrade_volume_mc_seeded(vol01, t, dose_mode=\"quarter\", enable_motion=False):\n    \"\"\"\n    Stochastic degradation matching training/CRM style.\n    Randomness is controlled by the external seed set before calling.\n    \"\"\"\n    D = vol01.shape[0]\n    out = np.empty_like(vol01, dtype=np.float32)\n\n    sigma = math.sqrt(max(1e-8, 2.0 * DIFFUSION_ALPHA * float(t)))\n\n    for z in range(D):\n        x = vol01[z].astype(np.float32)\n        x = cv2.GaussianBlur(\n            x, (0, 0),\n            sigmaX=sigma, sigmaY=sigma,\n            borderType=cv2.BORDER_REPLICATE\n        )\n\n        # Motion artefact (disabled by default)\n        if enable_motion:\n            pass\n\n        if dose_mode != \"clean\":\n            if dose_mode == \"extreme\":\n                peak = random.uniform(1000.0, 3000.0)\n                sigma_e = random.uniform(0.02, 0.04)\n            else:\n                peak = random.uniform(3000.0, 6000.0)\n                sigma_e = random.uniform(0.01, 0.02)\n\n            noisy_p = np.random.poisson(np.clip(x * peak, 0, None)).astype(np.float32) / peak\n            x = noisy_p + np.random.randn(*x.shape).astype(np.float32) * sigma_e\n\n        out[z] = np.clip(x, 0.0, 1.0).astype(np.float32)\n\n    return out\n\n@torch.no_grad()\ndef run_vultimate_mc(vol_deg01, t=EVAL_T, restore_batch=RESTORE_BATCH):\n    \"\"\"\n    FiLM + authority map compatible inference.\n    meta = [t_norm, do_motion, dose_clean, is_identity]\n    \"\"\"\n    vol_deg01 = np.asarray(vol_deg01, dtype=np.float32)\n    D = vol_deg01.shape[0]\n    out = vol_deg01.copy()\n\n    t_norm = np.float32(0.0 if t <= 0 else (float(t) / float(BLUR_T_MAX)))\n    meta_row = np.array([t_norm, 0.0, 0.0, 0.0], dtype=np.float32)\n\n    for s in range(0, D, restore_batch):\n        zs = list(range(s, min(D, s + restore_batch)))\n        inp_batch = []\n\n        for z in zs:\n            bp = vol_deg01[max(0, z - 1)]\n            bc = vol_deg01[z]\n            bn = vol_deg01[min(D - 1, z + 1)]\n            inp_batch.append(\n                np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32)\n            )\n\n        inp_t = torch.from_numpy(np.stack(inp_batch, axis=0)).to(device, non_blocking=True)\n        meta_t = torch.from_numpy(np.repeat(meta_row[None, :], len(zs), axis=0)).to(device, non_blocking=True)\n\n        with AMP_CTX():\n            pred_obj = MODEL_OBJ(inp_t, meta=meta_t)\n            if isinstance(pred_obj, (tuple, list)):\n                pred_b = pred_obj[0].float().cpu().numpy()[:, 0]\n            else:\n                pred_b = pred_obj.float().cpu().numpy()[:, 0]\n\n        for k, z in enumerate(zs):\n            out[z] = np.clip(pred_b[k], 0.0, 1.0).astype(np.float32)\n\n    return out\n\n# -----------------------------\n# 3) Build OOD CT/CTA pool\n# -----------------------------\nprint(\"=== Stage C | Monte Carlo Stability Test ===\")\n\nmeta = pd.read_csv(META_CSV)\nct_uids = set(meta[meta[\"Modality\"].astype(str).isin({\"CT\", \"CTA\"})][\"SeriesInstanceUID\"].astype(str).tolist())\n\ntrain_uid_set = set()\nfor p in [\n    \"/kaggle/working/vultimate_sleep_safe/train_uids_ultimate.csv\",\n    \"/kaggle/working/train_uids_ultimate.csv\",\n    \"/kaggle/input/datasets/mingzeli2009/train-uids-ultimate/train_uids_ultimate.csv\",\n    \"/kaggle/working/train_uids_ct_only.csv\",\n]:\n    if os.path.exists(p):\n        try:\n            train_uid_set = set(pd.read_csv(p)[\"SeriesInstanceUID\"].astype(str).tolist())\n            print(f\"[Train exclusion] loaded {len(train_uid_set)} train UIDs from: {p}\")\n            break\n        except Exception:\n            pass\n\nall_series_dirs = [u for u in os.listdir(RSNA_DATA_ROOT) if os.path.isdir(os.path.join(RSNA_DATA_ROOT, u))]\nood_pool = [u for u in all_series_dirs if (u in ct_uids) and (u not in dev_uid_set)]\n\nrandom.Random(2026).shuffle(ood_pool)\n\nmc_candidates = ood_pool[:max(N_MC_TO_RUN * 3, 300)]\n\n# -----------------------------\n# 4) Pre-load cases\n# -----------------------------\nmc_cases = []\nt_prep = time.time()\n\nfor uid in mc_candidates:\n    if len(mc_cases) >= N_MC_TO_RUN:\n        break\n\n    try:\n        vol = load_series_volume(uid, RSNA_DATA_ROOT, TARGET_SHAPE)\n    except TypeError:\n        try:\n            vol = load_series_volume(uid, RSNA_DATA_ROOT)\n        except TypeError:\n            vol = load_series_volume(uid)\n\n    if vol is None:\n        continue\n\n    gt_vol = vol.astype(np.float32)\n    p_gt = float(aneurysm_predict(vol01_to_flayer_uint8(gt_vol)))\n\n    mc_cases.append({\n        \"uid\": uid,\n        \"uid4\": uid_tail4(uid),\n        \"gt_vol\": gt_vol,\n        \"p_gt\": p_gt,\n    })\n\n    if (len(mc_cases) == 1) or (len(mc_cases) % 10 == 0) or (len(mc_cases) == N_MC_TO_RUN):\n        print(f\"[prep {len(mc_cases):03d}/{N_MC_TO_RUN}] UID:{uid_tail4(uid)} | p_gt={p_gt:.4f}\")\n\nprint(f\"MC selected cases: {len(mc_cases)} / {N_MC_TO_RUN}\")\nprint(f\"MC seeds: {MC_SEEDS}\")\nprint(f\"MC preparation elapsed: {(time.time()-t_prep)/60:.1f} min\")\n\n# -----------------------------\n# 5) Monte Carlo main loop\n# -----------------------------\nmc_raw_rows = []\nt_mc = time.time()\n\nfor ci, case in enumerate(mc_cases, 1):\n    uid = case[\"uid\"]\n    uid4 = case[\"uid4\"]\n    gt_vol = case[\"gt_vol\"]\n    p_gt = float(case[\"p_gt\"])\n\n    print(f\"\\n[{ci:03d}/{len(mc_cases)}] UID:{uid4} | p_gt={p_gt:.4f} | running {len(MC_SEEDS)} seeds ...\")\n    case_t0 = time.time()\n\n    case_bucket_votes = []\n\n    for s in MC_SEEDS:\n        random.seed(s)\n        np.random.seed(s)\n        torch.manual_seed(s)\n        if torch.cuda.is_available():\n            torch.cuda.manual_seed_all(s)\n\n        deg_vol = degrade_volume_mc_seeded(gt_vol, EVAL_T, dose_mode=EVAL_DOSE, enable_motion=False)\n        rec_vol = run_vultimate_mc(deg_vol, t=EVAL_T, restore_batch=RESTORE_BATCH)\n\n        p_deg = float(aneurysm_predict(vol01_to_flayer_uint8(deg_vol)))\n        p_rec = float(aneurysm_predict(vol01_to_flayer_uint8(rec_vol)))\n\n        bucket = mc_bucket(p_gt, p_deg, tau=0.03)\n        case_bucket_votes.append(bucket)\n\n        t_gain = float(calc_target_gain(p_gt, p_deg, p_rec))\n        a_gain = float(calc_abs_gain(p_gt, p_deg, p_rec))\n\n        mc_raw_rows.append({\n            \"uid_full\": uid,\n            \"uid4\": uid4,\n            \"seed\": int(s),\n            \"bucket\": bucket,\n            \"p_gt\": p_gt,\n            \"p_deg\": p_deg,\n            \"p_rec\": p_rec,\n            \"target_gain\": t_gain,\n            \"abs_gain\": a_gain,\n            \"target_positive\": int(t_gain > TARGET_POS_TH),\n            \"target_negative\": int(t_gain < TARGET_NEG_TH),\n            \"abs_positive\": int(a_gain > ABS_POS_TH),\n            \"abs_negative\": int(a_gain < ABS_NEG_TH),\n        })\n\n        del deg_vol, rec_vol\n        gc.collect()\n        if torch.cuda.is_available():\n            torch.cuda.empty_cache()\n\n    tmp = pd.DataFrame([r for r in mc_raw_rows if r[\"uid_full\"] == uid])\n    major_bucket = tmp[\"bucket\"].mode().iloc[0] if len(tmp) else \"NA\"\n\n    print(\n        \"  -> case done in {:.1f}s | bucket={} | \"\n        \"TGain mean={:+.4f}, std={:.4f}, pos/neg={:.2f}/{:.2f} | \"\n        \"AGain mean={:+.4f}, pos/neg={:.2f}/{:.2f}\".format(\n            time.time()-case_t0,\n            major_bucket,\n            tmp[\"target_gain\"].mean(), tmp[\"target_gain\"].std(ddof=0),\n            (tmp[\"target_gain\"] > TARGET_POS_TH).mean(), (tmp[\"target_gain\"] < TARGET_NEG_TH).mean(),\n            tmp[\"abs_gain\"].mean(),\n            (tmp[\"abs_gain\"] > ABS_POS_TH).mean(), (tmp[\"abs_gain\"] < ABS_NEG_TH).mean(),\n        )\n    )\n\ndf_mc_raw = pd.DataFrame(mc_raw_rows)\n\n# -----------------------------\n# 6) Aggregation\n# -----------------------------\ndef agg_mc_case(g):\n    t = g[\"target_gain\"].to_numpy(dtype=float)\n    a = g[\"abs_gain\"].to_numpy(dtype=float)\n    pdeg = g[\"p_deg\"].to_numpy(dtype=float)\n    prec = g[\"p_rec\"].to_numpy(dtype=float)\n    pgt = float(g[\"p_gt\"].iloc[0])\n\n    t_pos = t > TARGET_POS_TH\n    t_neg = t < TARGET_NEG_TH\n    a_pos = a > ABS_POS_TH\n    a_neg = a < ABS_NEG_TH\n\n    bucket_mode = g[\"bucket\"].mode().iloc[0] if len(g[\"bucket\"].mode()) else \"NA\"\n\n    return pd.Series({\n        \"n_runs\": int(len(g)),\n        \"bucket_major\": bucket_mode,\n        \"p_gt\": pgt,\n\n        \"target_gain_mean\": float(np.mean(t)),\n        \"target_gain_std\": float(np.std(t, ddof=0)),\n        \"target_gain_min\": float(np.min(t)),\n        \"target_gain_max\": float(np.max(t)),\n        \"target_pos_rate\": float(np.mean(t_pos)),\n        \"target_neg_rate\": float(np.mean(t_neg)),\n        \"target_flip\": bool(np.any(t_pos) and np.any(t_neg)),\n\n        \"abs_gain_mean\": float(np.mean(a)),\n        \"abs_gain_std\": float(np.std(a, ddof=0)),\n        \"abs_gain_min\": float(np.min(a)),\n        \"abs_gain_max\": float(np.max(a)),\n        \"abs_pos_rate\": float(np.mean(a_pos)),\n        \"abs_neg_rate\": float(np.mean(a_neg)),\n        \"abs_flip\": bool(np.any(a_pos) and np.any(a_neg)),\n\n        \"p_deg_std\": float(np.std(pdeg, ddof=0)),\n        \"p_rec_std\": float(np.std(prec, ddof=0)),\n    })\n\ndf_mc_agg = (\n    df_mc_raw.groupby([\"uid_full\", \"uid4\"], as_index=False)\n    .apply(agg_mc_case)\n    .reset_index(drop=True)\n)\n\nmc_raw_path = os.path.join(OUTDIR, \"mc_noise_100cases_10seeds_raw.csv\")\nmc_agg_path = os.path.join(OUTDIR, \"mc_noise_100cases_10seeds_agg.csv\")\ndf_mc_raw.to_csv(mc_raw_path, index=False)\ndf_mc_agg.to_csv(mc_agg_path, index=False)\n\n# -----------------------------\n# 7) Summary statistics\n# -----------------------------\nn = len(df_mc_agg)\n\nstable_pos = ((df_mc_agg[\"target_pos_rate\"] > 0) & (df_mc_agg[\"target_neg_rate\"] == 0)).sum()\nflip       = ((df_mc_agg[\"target_pos_rate\"] > 0) & (df_mc_agg[\"target_neg_rate\"] > 0)).sum()\nstable_neg = ((df_mc_agg[\"target_pos_rate\"] == 0) & (df_mc_agg[\"target_neg_rate\"] > 0)).sum()\nneutral    = ((df_mc_agg[\"target_pos_rate\"] == 0) & (df_mc_agg[\"target_neg_rate\"] == 0)).sum()\n\nmc_summary = {\n    \"n_cases\": n,\n    \"n_seeds_per_case\": len(MC_SEEDS),\n    \"n_total_runs\": len(df_mc_raw),\n\n    \"target_gain_mean(run-level)\": float(df_mc_raw[\"target_gain\"].mean()),\n    \"target_positive_rate(run-level)\": float((df_mc_raw[\"target_gain\"] > TARGET_POS_TH).mean()),\n    \"target_negative_rate(run-level)\": float((df_mc_raw[\"target_gain\"] < TARGET_NEG_TH).mean()),\n\n    \"abs_gain_mean(run-level)\": float(df_mc_raw[\"abs_gain\"].mean()),\n\n    \"cases_with_target_flip\": int(df_mc_agg[\"target_flip\"].sum()),\n    \"stable_positive\": int(stable_pos),\n    \"noise_sensitive_flip\": int(flip),\n    \"stable_negative\": int(stable_neg),\n    \"neutral\": int(neutral),\n\n    \"elapsed_min\": round((time.time() - t_mc) / 60.0, 1),\n}\n\nprint(\"\\n\" + \"=\"*100)\nprint(f\"Monte Carlo Stability Summary ({n} cases × {len(MC_SEEDS)} seeds)\")\nprint(\"=\"*100)\nfor k, v in mc_summary.items():\n    if isinstance(v, float):\n        print(f\"{k:>45}: {v:.4f}\")\n    else:\n        print(f\"{k:>45}: {v}\")\n\nprint(f\"\\nCase classification:\")\nprint(f\"   Stable positive:    {stable_pos}/{n} ({stable_pos/n*100:.0f}%)\")\nprint(f\"   Noise-sensitive:    {flip}/{n} ({flip/n*100:.0f}%)\")\nprint(f\"   Stable negative:    {stable_neg}/{n} ({stable_neg/n*100:.0f}%)\")\nprint(f\"   Neutral:            {neutral}/{n} ({neutral/n*100:.0f}%)\")\n\nprint(\"\\n[Top noise-sensitive cases by target_gain_std]\")\ndisplay(\n    df_mc_agg.sort_values([\"target_gain_std\", \"target_neg_rate\"], ascending=[False, False])\n    .head(10)\n    .reset_index(drop=True)\n)\n\nprint(\"\\n[Top unstable / failure-prone cases by target_neg_rate]\")\ndisplay(\n    df_mc_agg.sort_values([\"target_neg_rate\", \"target_gain_std\"], ascending=[False, False])\n    .head(10)\n    .reset_index(drop=True)\n)\n\nprint(\"\\n[Bucket breakdown]\")\ndisplay(\n    df_mc_agg.groupby(\"bucket_major\", as_index=False)\n    .agg(\n        n_cases=(\"uid_full\", \"size\"),\n        mean_target_gain=(\"target_gain_mean\", \"mean\"),\n        mean_target_std=(\"target_gain_std\", \"mean\"),\n        cases_with_flip=(\"target_flip\", \"sum\"),\n        mean_prec_std=(\"p_rec_std\", \"mean\"),\n    )\n)\n\nprint(\"\\nsaved:\", mc_raw_path)\nprint(\"saved:\", mc_agg_path)\nprint(\"\\n✅ Monte Carlo complete.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T01:52:28.522378Z","iopub.execute_input":"2026-03-15T01:52:28.523054Z","iopub.status.idle":"2026-03-15T06:01:24.882122Z","shell.execute_reply.started":"2026-03-15T01:52:28.523017Z","shell.execute_reply":"2026-03-15T06:01:24.88128Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Results Interpretation — Monte Carlo Stability Test\n\n### 1. What is this experiment testing?\n\nThe Monte Carlo stability test asks a different question from the main CRM table.\n\nThe CRM summary tells us whether V-Ultimate is beneficial **on average** across a fixed evaluation set.  \nThe Monte Carlo analysis asks:\n\n> **If the degradation noise realization changes while the same underlying case is kept fixed, does the model still behave consistently?**\n\nThis is important because a restoration method can look strong in one deterministic evaluation, yet be unstable when the exact corruption pattern varies.\n\nIn other words, this section measures **robustness to stochastic degradation**, not just average performance.\n\n---\n\n### 2. Overall result: the method remains net-positive under repeated stochastic perturbations\n\nAcross **100 cases × 10 noise seeds = 1000 total runs**, the run-level summary is:\n\n- **Mean target gain = 0.0335**\n- **Target-positive rate = 71.7%**\n- **Target-negative rate = 20.3%**\n- **Mean absolute gain = 0.0549**\n\nThis is the most important first conclusion:\n\n> **Even after repeatedly perturbing each case with different stochastic degradation realizations, V-Ultimate remains net-positive overall.**\n\nSo the main CRM conclusion is **not** a one-seed accident.  \nThe model still shows a positive average rescue effect under repeated noisy corruption.\n\n---\n\n### 3. But the stability is not uniform: most cases are noise-sensitive to some degree\n\nThe case-level stability breakdown is:\n\n- **Stable positive:** 36 / 100 (**36%**)\n- **Noise-sensitive:** 62 / 100 (**62%**)\n- **Stable negative:** 2 / 100 (**2%**)\n- **Neutral:** 0 / 100 (**0%**)\n\nThis means the correct interpretation is **not**:\n\n> “The model behaves identically under all noise realizations.”\n\nInstead, the correct interpretation is:\n\n> **The model is beneficial overall, but its case-level response is often sensitive to the exact stochastic corruption pattern.**\n\nThat is why this section is so important. It shows that the method is **robust in expectation**, but **not uniformly deterministic at the individual-case level**.\n\n---\n\n### 4. The most encouraging result: clinically important harmed-positive cases remain positive on average\n\nThe most important bucket is still **`harmed_positive`**, where degradation suppresses a clinically relevant signal that should ideally be recovered.\n\nFor this bucket:\n\n- **n_cases = 80**\n- **mean_target_gain = 0.0442**\n- **mean_target_std = 0.0324**\n- **cases_with_flip = 45**\n- **mean_pred_std = 0.0243**\n\nThis is the core takeaway:\n\n> **In the clinically important harmed-positive cases, the average effect remains clearly positive, even though many individual cases still exhibit seed-to-seed fluctuation.**\n\nSo the model is doing the right thing in the most relevant bucket **on average**, but some cases sit close enough to the decision boundary that different noise realizations can change whether the outcome is counted as positive or negative.\n\nThat is exactly why many cases are labeled **noise-sensitive**, rather than simply “successful” or “failed.”\n\n---\n\n### 5. Neutral and overcall cases remain the main weakness\n\nThe weaker stability appears in the buckets where recovery is less clearly needed:\n\n#### `neutral_other`\n- **n_cases = 18**\n- **mean_target_gain = -0.0093**\n- **mean_target_std = 0.0271**\n- **cases_with_flip = 15**\n\n#### `overcall_positive`\n- **n_cases = 2**\n- **mean_target_gain = -0.0126**\n- **mean_target_std = 0.0309**\n- **cases_with_flip = 2**\n\nThese results are consistent with the earlier CRM interpretation:\n\n> **V-Ultimate is strongest when the target signal has genuinely been damaged, but it is less reliable when the case is neutral or already overcalled.**\n\nThis matters because it shows that the model’s current limitation is not random across all cases.  \nIts instability is concentrated in the scenarios where the restoration objective is inherently more ambiguous.\n\n---\n\n### 6. The meaning of “noise-sensitive” is important\n\nA “noise-sensitive” case does **not** necessarily mean the model is bad.  \nIt means that across different stochastic corruptions of the same underlying case, the direction or strength of the rescue effect can change.\n\nFor example:\n\n- some cases remain positive, but with variable strength;\n- some cases occasionally cross from positive to negative;\n- some cases are near the clinical decision boundary and are therefore especially sensitive.\n\nSo the right interpretation is:\n\n> **Noise-sensitive cases are not all failures; they are cases where restoration outcome depends meaningfully on the exact degradation realization.**\n\nThis is a more nuanced and scientifically useful conclusion than simply labeling everything as success or failure.\n\n---\n\n### 7. The top unstable cases show where the method still needs improvement\n\nThe most unstable examples are informative because they reveal the current edge cases of the method.\n\nTwo patterns stand out:\n\n#### A. Some `harmed_positive` cases are beneficial on average, but still flip across seeds\nExamples such as `3439`, `0800`, `4778`, and `6455` have:\n\n- positive mean target gain,\n- but substantial standard deviation,\n- and seed-level sign flips.\n\nThis suggests that:\n\n> the model often moves in the right direction, but its response is still somewhat sensitive when the corruption severity lands near a difficult boundary.\n\n#### B. Many `neutral_other` cases are consistently fragile\nExamples such as `6332`, `3779`, `3477`, and `3141` show:\n\n- negative or near-zero mean target gain,\n- frequent flips,\n- and clear instability across seeds.\n\nThis reinforces the earlier finding that:\n\n> **the model is currently much better at rescuing damaged signal than at deciding when to leave an already acceptable case alone.**\n\n---\n\n### 8. What does this mean for the overall claim of the project?\n\nThe Monte Carlo result does **not** weaken the project’s central claim.  \nInstead, it sharpens it.\n\nA naive interpretation would be:\n\n> “Because many cases are noise-sensitive, the method is unreliable.”\n\nThat would be too strong.\n\nA better interpretation is:\n\n> **The method is robust at the population level and clearly beneficial in the most clinically relevant harmed-positive bucket, but it is not yet fully stable at the individual-case level under stochastic degradation variation.**\n\nThis is actually a meaningful scientific result, because it shows both:\n\n1. **why the method is promising**, and  \n2. **where its current limits still are.**\n\n---\n\n### 9. Final interpretation of this section\n\nThe Monte Carlo stability test supports the following conclusion:\n\n> **V-Ultimate’s CRM benefit is not a single-seed artifact: under repeated stochastic degradation, the method remains net-positive overall and especially positive in clinically important harmed-positive cases.**\n\nAt the same time, it also reveals an important limitation:\n\n> **Case-level behavior is often noise-sensitive, particularly in neutral or ambiguous cases, indicating that the model still needs stronger mechanisms for deciding when not to edit.**\n\nTherefore, this section strengthens the project in two ways:\n\n- it confirms that the observed rescue effect is **real and reproducible in expectation**;\n- it identifies **precisely where further safety and stability improvements are needed**.\n\nThe correct summary is not:\n\n> “The model is perfectly stable.”\n\nIt is:\n\n> **The model is directionally stable at the population level, strongest where clinically relevant signal is genuinely suppressed, but still sensitive in ambiguous cases under stochastic corruption.**","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# Mayo external paired evaluation — 10 unique matched series\n# - strict matching by relative directory\n# - exactly one center window per matched series\n# - no repeated windows from the same series\n# Depends on:\n#   - model_25d, device, AMP_CTX\n#   - Q_DIR, F_DIR, MAYO_OUTDIR\n# ============================================================\n\nMAYO_ROOT = \"/kaggle/input/datasets/andrewmvd/ct-low-dose-reconstruction/CT_low_dose_reconstruction_dataset/Original Data\"\nQ_DIR = os.path.join(MAYO_ROOT, \"Quarter Dose\")\nF_DIR = os.path.join(MAYO_ROOT, \"Full Dose\")\n\nMAYO_OUTDIR = Path(\"/kaggle/working/vultimate_scientific_eval/mayo_external_eval\")\nMAYO_OUTDIR.mkdir(parents=True, exist_ok=True)\n\nprint(\"MAYO_ROOT exists?\", os.path.exists(MAYO_ROOT))\nprint(\"Q_DIR exists?\", os.path.exists(Q_DIR))\nprint(\"F_DIR exists?\", os.path.exists(F_DIR))\nprint(\"MAYO_OUTDIR:\", MAYO_OUTDIR)\n\nimport os, gc, math, time\nfrom pathlib import Path\nfrom contextlib import nullcontext\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport torch\n\ntry:\n    from IPython.display import display\nexcept Exception:\n    display = print\n\n# ------------------------------------------------------------\n# 0) Config / checks\n# ------------------------------------------------------------\nrequired_globals = [\n    \"model_25d\", \"device\", \"AMP_CTX\",\n    \"Q_DIR\", \"F_DIR\", \"MAYO_OUTDIR\",\n    \"HU_MIN\", \"HU_MAX\", \"HU_RANGE\",\n    \"BLUR_T_MAX\", \"RESTORE_BATCH\",\n]\nmissing = [k for k in required_globals if k not in globals()]\nif missing:\n    raise RuntimeError(f\"Missing required globals: {missing}\")\n\nOUTDIR = Path(MAYO_OUTDIR)\nOUTDIR.mkdir(parents=True, exist_ok=True)\n\nN_CASES_MAIN = int(globals().get(\"N_CASES_MAIN\", 10))   # unique matched series\nVOL_DEPTH = int(globals().get(\"VOL_DEPTH\", 128))\nTARGET_H = int(globals().get(\"TARGET_H\", 448))\nTARGET_W = int(globals().get(\"TARGET_W\", 448))\n\nT_INFER_NORM = 0.05\nT_INFER = T_INFER_NORM * float(BLUR_T_MAX)\nMAYO_USE_GAUSSIAN = True\nMAYO_GAUSSIAN_SIGMA = 0.8\nBOOT_N = 2000\nSEED_MAYO = 2026\n\nprint(\"Q_DIR exists?\", os.path.exists(Q_DIR))\nprint(\"F_DIR exists?\", os.path.exists(F_DIR))\nprint(\"N_CASES_MAIN:\", N_CASES_MAIN)\n\n# ------------------------------------------------------------\n# 1) Helpers\n# ------------------------------------------------------------\ndef is_dicom_like(fname: str):\n    if fname.startswith(\".\"):\n        return False\n    if fname.endswith((\".dcm\", \".DCM\", \".ima\", \".IMA\")):\n        return True\n    if \".\" not in fname:\n        return True\n    return False\n\ndef dcm_sort_key(fp: str):\n    try:\n        ds = pydicom.dcmread(fp, stop_before_pixels=True, force=True)\n\n        inst = getattr(ds, \"InstanceNumber\", None)\n        if inst is not None:\n            try:\n                return (0, int(inst), fp)\n            except Exception:\n                pass\n\n        ipp = getattr(ds, \"ImagePositionPatient\", None)\n        if ipp is not None and len(ipp) >= 3:\n            try:\n                return (1, float(ipp[2]), fp)\n            except Exception:\n                pass\n\n        sl = getattr(ds, \"SliceLocation\", None)\n        if sl is not None:\n            try:\n                return (2, float(sl), fp)\n            except Exception:\n                pass\n    except Exception:\n        pass\n\n    return (9, fp)\n\ndef list_dicom_series(root_dir):\n    series = {}\n    for r, _, fs in os.walk(root_dir):\n        dcm_files = [os.path.join(r, f) for f in fs if is_dicom_like(f)]\n        if len(dcm_files) == 0:\n            continue\n        rel = os.path.relpath(r, root_dir)\n        series[rel] = sorted(dcm_files, key=dcm_sort_key)\n    return series\n\ndef dcm_to_hu(ds):\n    arr = ds.pixel_array.astype(np.float32)\n    slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n    return arr * slope + intercept\n\ndef hu_to_01(hu):\n    return np.clip((hu - HU_MIN) / HU_RANGE, 0.0, 1.0).astype(np.float32)\n\ndef hu01_to_hu(x01):\n    return np.asarray(x01, dtype=np.float32) * HU_RANGE + HU_MIN\n\ndef psnr01(a, b):\n    a = np.asarray(a, dtype=np.float32)\n    b = np.asarray(b, dtype=np.float32)\n    mse = float(np.mean((a - b) ** 2))\n    return 99.0 if mse <= 0 else 10.0 * math.log10(1.0 / mse)\n\ndef mae01(a, b):\n    return float(np.mean(np.abs(np.asarray(a, dtype=np.float32) - np.asarray(b, dtype=np.float32))))\n\ndef ssim_fast01_2d(a, b):\n    a = a.astype(np.float32)\n    b = b.astype(np.float32)\n    C1, C2 = 0.01 ** 2, 0.03 ** 2\n    mu_a = cv2.GaussianBlur(a, (11, 11), 1.5)\n    mu_b = cv2.GaussianBlur(b, (11, 11), 1.5)\n    mu_a2, mu_b2, mu_ab = mu_a * mu_a, mu_b * mu_b, mu_a * mu_b\n    sigma_a2 = cv2.GaussianBlur(a * a, (11, 11), 1.5) - mu_a2\n    sigma_b2 = cv2.GaussianBlur(b * b, (11, 11), 1.5) - mu_b2\n    sigma_ab = cv2.GaussianBlur(a * b, (11, 11), 1.5) - mu_ab\n    ssim_map = ((2 * mu_ab + C1) * (2 * sigma_ab + C2)) / (\n        (mu_a2 + mu_b2 + C1) * (sigma_a2 + sigma_b2 + C2) + 1e-8\n    )\n    return float(np.mean(ssim_map))\n\ndef build_paired_case_manifest_unique(q_root, f_root, depth=128, n_target=10):\n    q_series = list_dicom_series(q_root)\n    f_series = list_dicom_series(f_root)\n    common_rels = sorted(set(q_series.keys()) & set(f_series.keys()))\n\n    cases = []\n    for rel in common_rels:\n        qf = q_series[rel]\n        ff = f_series[rel]\n        usable_len = min(len(qf), len(ff))\n        if usable_len < depth:\n            continue\n\n        center_start = max(0, (usable_len - depth) // 2)\n        cases.append({\n            \"rel_dir\": rel,\n            \"start_idx\": int(center_start),\n            \"usable_len\": int(usable_len),\n            \"q_files\": qf,\n            \"f_files\": ff,\n        })\n\n    if len(cases) == 0:\n        raise RuntimeError(\"No matched Quarter/Full series with enough slices were found.\")\n\n    # prefer longer usable series; one window per unique matched series only\n    cases = sorted(cases, key=lambda x: (-x[\"usable_len\"], x[\"rel_dir\"]))[:n_target]\n\n    out = []\n    for i, m in enumerate(cases, 1):\n        out.append({\n            \"case_id\": i,\n            \"rel_dir\": m[\"rel_dir\"],\n            \"start_idx\": m[\"start_idx\"],\n            \"usable_len\": m[\"usable_len\"],\n            \"q_files\": m[\"q_files\"],\n            \"f_files\": m[\"f_files\"],\n        })\n    return out\n\ndef load_paired_case_window(q_files, f_files, start_idx, depth=128, out_hw=(448, 448)):\n    q_hu_list, f_hu_list = [], []\n    max_len = min(len(q_files), len(f_files))\n\n    for k in range(depth):\n        idx = start_idx + k\n        if idx < 0 or idx >= max_len:\n            break\n\n        ds_q = pydicom.dcmread(q_files[idx], force=True)\n        ds_f = pydicom.dcmread(f_files[idx], force=True)\n\n        hq = dcm_to_hu(ds_q)\n        hf = dcm_to_hu(ds_f)\n\n        if hq.shape != out_hw:\n            hq = cv2.resize(hq, (out_hw[1], out_hw[0]), interpolation=cv2.INTER_LINEAR)\n        if hf.shape != out_hw:\n            hf = cv2.resize(hf, (out_hw[1], out_hw[0]), interpolation=cv2.INTER_LINEAR)\n\n        q_hu_list.append(hq.astype(np.float32))\n        f_hu_list.append(hf.astype(np.float32))\n\n    q_hu = np.stack(q_hu_list, axis=0).astype(np.float32)\n    f_hu = np.stack(f_hu_list, axis=0).astype(np.float32)\n    q01 = hu_to_01(q_hu)\n    f01 = hu_to_01(f_hu)\n    return q_hu, f_hu, q01, f01\n\n@torch.no_grad()\ndef run_vultimate_mayo(vol_deg01, t=T_INFER, restore_batch=RESTORE_BATCH):\n    vol_deg01 = np.asarray(vol_deg01, dtype=np.float32)\n    D = vol_deg01.shape[0]\n    out = vol_deg01.copy()\n\n    t_norm = np.float32(0.0 if t <= 0 else float(t) / float(BLUR_T_MAX))\n    meta_row = np.array([t_norm, 0.0, 0.0, 0.0], dtype=np.float32)\n\n    for s in range(0, D, restore_batch):\n        zs = list(range(s, min(D, s + restore_batch)))\n        inp_batch = []\n\n        for z in zs:\n            bp = vol_deg01[max(0, z - 1)]\n            bc = vol_deg01[z]\n            bn = vol_deg01[min(D - 1, z + 1)]\n            inp_batch.append(np.stack([bp, bc, bn, np.full_like(bc, t_norm)], axis=0).astype(np.float32))\n\n        inp_t = torch.from_numpy(np.stack(inp_batch, axis=0)).to(device, non_blocking=True)\n        meta_t = torch.from_numpy(np.repeat(meta_row[None, :], len(zs), axis=0)).to(device, non_blocking=True)\n\n        with AMP_CTX():\n            pred_obj = model_25d(inp_t, meta=meta_t)\n            pred = pred_obj[0] if isinstance(pred_obj, (tuple, list)) else pred_obj\n            pred = pred.float().cpu().numpy()[:, 0]\n\n        for k, z in enumerate(zs):\n            out[z] = np.clip(pred[k], 0.0, 1.0).astype(np.float32)\n\n    return out\n\ndef gaussian_baseline_3d(vol01, sigma=0.8):\n    out = np.empty_like(vol01, dtype=np.float32)\n    for z in range(vol01.shape[0]):\n        out[z] = cv2.GaussianBlur(\n            vol01[z], (0, 0), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REPLICATE\n        )\n    return np.clip(out, 0.0, 1.0).astype(np.float32)\n\ndef vol_metrics(x01, ref01, q01_for_change):\n    z_idx = list(range(0, x01.shape[0], 8))\n    ssim_vals = [ssim_fast01_2d(x01[z], ref01[z]) for z in z_idx]\n    return {\n        \"psnr\": psnr01(x01, ref01),\n        \"mae\": mae01(x01, ref01),\n        \"ssim\": float(np.mean(ssim_vals)),\n        \"mean_abs_change_vs_quarter\": float(np.mean(np.abs(x01 - q01_for_change))),\n        \"max_abs_change_vs_quarter\": float(np.max(np.abs(x01 - q01_for_change))),\n    }\n\ndef paired_bootstrap_ci(x, y, n_boot=2000, seed=2026):\n    x = np.asarray(x, dtype=np.float32)\n    y = np.asarray(y, dtype=np.float32)\n    d = x - y\n    mean_d = float(np.mean(d))\n\n    rng = np.random.default_rng(seed)\n    boots = []\n    n = len(d)\n    for _ in range(n_boot):\n        idx = rng.integers(0, n, size=n)\n        boots.append(float(np.mean(d[idx])))\n\n    lo = float(np.percentile(boots, 2.5))\n    hi = float(np.percentile(boots, 97.5))\n    return mean_d, lo, hi\n\n# ------------------------------------------------------------\n# 2) Build strict unique-series manifest\n# ------------------------------------------------------------\nmayo_case_manifest = build_paired_case_manifest_unique(\n    Q_DIR, F_DIR, depth=VOL_DEPTH, n_target=N_CASES_MAIN\n)\n\ndf_manifest = pd.DataFrame([\n    {\n        \"case_id\": m[\"case_id\"],\n        \"rel_dir\": m[\"rel_dir\"],\n        \"start_idx\": m[\"start_idx\"],\n        \"usable_len\": m[\"usable_len\"],\n        \"q_nfiles\": len(m[\"q_files\"]),\n        \"f_nfiles\": len(m[\"f_files\"]),\n    }\n    for m in mayo_case_manifest\n])\n\nassert df_manifest[\"rel_dir\"].nunique() == len(df_manifest), \"Expected unique matched series only.\"\n\ndf_manifest.to_csv(OUTDIR / \"mayo_case_manifest.csv\", index=False)\n\nprint(\"\\n=== Mayo case manifest ===\")\ndisplay(df_manifest)\nprint(\"Total selected windows:\", len(df_manifest))\nprint(\"Unique matched series :\", df_manifest[\"rel_dir\"].nunique())\n\n# ------------------------------------------------------------\n# 3) Main experiment\n# ------------------------------------------------------------\nrows = []\nt0 = time.time()\n\nfor i, case_info in enumerate(mayo_case_manifest, 1):\n    cid = int(case_info[\"case_id\"])\n    rel_dir = case_info[\"rel_dir\"]\n    st = int(case_info[\"start_idx\"])\n\n    print(f\"\\n[{i}/{len(mayo_case_manifest)}] case_id={cid:02d} | rel_dir={rel_dir} | start={st}\")\n\n    q_hu, f_hu, q01, f01 = load_paired_case_window(\n        case_info[\"q_files\"], case_info[\"f_files\"], st, depth=VOL_DEPTH, out_hw=(TARGET_H, TARGET_W)\n    )\n\n    m_q = vol_metrics(q01, f01, q01)\n    rows.append({\n        \"case_id\": cid, \"rel_dir\": rel_dir, \"start_idx\": st, \"variant\": \"Quarter\",\n        **m_q,\n    })\n\n    if MAYO_USE_GAUSSIAN:\n        g01 = gaussian_baseline_3d(q01, sigma=MAYO_GAUSSIAN_SIGMA)\n        m_g = vol_metrics(g01, f01, q01)\n        rows.append({\n            \"case_id\": cid, \"rel_dir\": rel_dir, \"start_idx\": st, \"variant\": \"Gaussian\",\n            **m_g,\n        })\n\n    rec01 = run_vultimate_mayo(q01, t=T_INFER, restore_batch=RESTORE_BATCH)\n    m_u = vol_metrics(rec01, f01, q01)\n    rows.append({\n        \"case_id\": cid, \"rel_dir\": rel_dir, \"start_idx\": st, \"variant\": \"V-Ultimate\",\n        **m_u,\n    })\n\n    line = f\"  Quarter PSNR={m_q['psnr']:.2f}\"\n    if MAYO_USE_GAUSSIAN:\n        line += f\" | Gaussian={m_g['psnr']:.2f}\"\n    line += f\" | V-Ultimate={m_u['psnr']:.2f}\"\n    print(line)\n\n    del q_hu, f_hu, q01, f01, rec01\n    if MAYO_USE_GAUSSIAN:\n        del g01\n    gc.collect()\n    if torch.cuda.is_available():\n        torch.cuda.empty_cache()\n\ndf_raw = pd.DataFrame(rows)\ndf_raw.to_csv(OUTDIR / \"mayo_generalization_raw.csv\", index=False)\n\ndf_sum = (\n    df_raw.groupby(\"variant\", as_index=False)\n    .agg(\n        n_case=(\"case_id\", \"nunique\"),\n        psnr_mean=(\"psnr\", \"mean\"),\n        psnr_std=(\"psnr\", \"std\"),\n        mae_mean=(\"mae\", \"mean\"),\n        mae_std=(\"mae\", \"std\"),\n        ssim_mean=(\"ssim\", \"mean\"),\n        ssim_std=(\"ssim\", \"std\"),\n        mean_abs_change_vs_quarter=(\"mean_abs_change_vs_quarter\", \"mean\"),\n        max_abs_change_vs_quarter=(\"max_abs_change_vs_quarter\", \"mean\"),\n    )\n    .sort_values(\"psnr_mean\", ascending=False)\n    .reset_index(drop=True)\n)\ndf_sum.to_csv(OUTDIR / \"mayo_generalization_summary.csv\", index=False)\n\nprint(\"\\n=== Mayo Generalization Summary ===\")\ndisplay(df_sum)\n\n# ------------------------------------------------------------\n# 4) Paired comparison + bootstrap CI\n# ------------------------------------------------------------\nq = df_raw[df_raw[\"variant\"] == \"Quarter\"][[\"case_id\", \"psnr\", \"mae\", \"ssim\"]].rename(\n    columns={\"psnr\": \"psnr_q\", \"mae\": \"mae_q\", \"ssim\": \"ssim_q\"}\n)\n\npaired_rows = []\nboot_rows = []\n\nfor v in sorted(df_raw[\"variant\"].unique()):\n    if v == \"Quarter\":\n        continue\n\n    d = df_raw[df_raw[\"variant\"] == v][[\"case_id\", \"psnr\", \"mae\", \"ssim\"]].rename(\n        columns={\"psnr\": \"psnr_v\", \"mae\": \"mae_v\", \"ssim\": \"ssim_v\"}\n    )\n    m = q.merge(d, on=\"case_id\", how=\"inner\").sort_values(\"case_id\").reset_index(drop=True)\n\n    dpsnr = m[\"psnr_v\"].values - m[\"psnr_q\"].values\n    dmae = m[\"mae_q\"].values - m[\"mae_v\"].values\n    dssim = m[\"ssim_v\"].values - m[\"ssim_q\"].values\n\n    paired_rows.append({\n        \"variant\": v,\n        \"n_case\": len(m),\n        \"ΔPSNR_vs_Quarter\": float(np.mean(dpsnr)),\n        \"ΔMAE_vs_Quarter\": float(np.mean(dmae)),\n        \"ΔSSIM_vs_Quarter\": float(np.mean(dssim)),\n        \"PSNR_win_rate\": float((m[\"psnr_v\"] > m[\"psnr_q\"]).mean()),\n        \"MAE_win_rate\": float((m[\"mae_v\"] < m[\"mae_q\"]).mean()),\n        \"SSIM_win_rate\": float((m[\"ssim_v\"] > m[\"ssim_q\"]).mean()),\n    })\n\n    psnr_mean, psnr_lo, psnr_hi = paired_bootstrap_ci(\n        m[\"psnr_v\"].values, m[\"psnr_q\"].values, n_boot=BOOT_N, seed=SEED_MAYO + 1\n    )\n    mae_mean, mae_lo, mae_hi = paired_bootstrap_ci(\n        m[\"mae_q\"].values, m[\"mae_v\"].values, n_boot=BOOT_N, seed=SEED_MAYO + 2\n    )\n    ssim_mean, ssim_lo, ssim_hi = paired_bootstrap_ci(\n        m[\"ssim_v\"].values, m[\"ssim_q\"].values, n_boot=BOOT_N, seed=SEED_MAYO + 3\n    )\n\n    boot_rows.extend([\n        {\"variant\": v, \"metric\": \"ΔPSNR_vs_Quarter\", \"mean\": psnr_mean, \"ci95_lo\": psnr_lo, \"ci95_hi\": psnr_hi},\n        {\"variant\": v, \"metric\": \"ΔMAE_vs_Quarter\",  \"mean\": mae_mean,  \"ci95_lo\": mae_lo,  \"ci95_hi\": mae_hi},\n        {\"variant\": v, \"metric\": \"ΔSSIM_vs_Quarter\", \"mean\": ssim_mean, \"ci95_lo\": ssim_lo, \"ci95_hi\": ssim_hi},\n    ])\n\ndf_pair = pd.DataFrame(paired_rows).sort_values(\"ΔPSNR_vs_Quarter\", ascending=False).reset_index(drop=True)\ndf_pair.to_csv(OUTDIR / \"mayo_generalization_paired_vs_quarter.csv\", index=False)\n\ndf_boot = pd.DataFrame(boot_rows)\ndf_boot.to_csv(OUTDIR / \"mayo_generalization_bootstrap_vs_quarter.csv\", index=False)\n\nprint(\"\\n=== Paired vs Quarter ===\")\ndisplay(df_pair)\n\nprint(\"\\n=== Bootstrap CI vs Quarter ===\")\ndisplay(df_boot)\n\nprint(\"\\nSaved:\")\nfor p in sorted(OUTDIR.glob(\"*.csv\")):\n    print(\" -\", p)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T06:28:46.843364Z","iopub.execute_input":"2026-03-15T06:28:46.844289Z","iopub.status.idle":"2026-03-15T06:36:29.194515Z","shell.execute_reply.started":"2026-03-15T06:28:46.844255Z","shell.execute_reply":"2026-03-15T06:36:29.193851Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Results Interpretation — Mayo External Paired Evaluation\n\n### 1. What does this experiment test?\n\nThis experiment is the key **sim-to-real transfer test** in the project.\n\nThe model was trained using **synthetic supervision**, but here it is evaluated on **real external paired low-dose CT data** from the Mayo dataset. Because each quarter-dose volume has a matched full-dose reference, this experiment can directly measure whether the restoration learned from synthetic data transfers to real clinical acquisition.\n\nSo the core question here is:\n\n> **Can a restoration model trained with physics-based synthetic supervision improve real external low-dose CT, rather than only working on synthetic test data?**\n\n---\n\n### 2. Dataset construction is strict and conservative\n\nThis evaluation uses:\n\n- **10 unique matched Mayo series**\n- **strict Quarter/Full pairing by relative directory**\n- **one center window per series**\n- **no repeated windows from the same series**\n\nThat matters because it prevents inflating the sample size by extracting many overlapping subvolumes from the same patient or series.  \nIn other words, this is a **small but clean external paired evaluation**, not a repeated-window benchmark.\n\nAnother notable detail is that the 10 cases are split across:\n\n- **Sharp kernel (D45)**\n- **Soft kernel (B30)**\n\nThis is useful because it introduces modest reconstruction-domain diversity within the external benchmark.\n\n---\n\n### 3. Main result: V-Ultimate clearly improves real Mayo quarter-dose CT\n\nCompared with the original quarter-dose input:\n\n- **Quarter mean PSNR = 39.49 dB**\n- **V-Ultimate mean PSNR = 41.84 dB**\n- **ΔPSNR = +2.35 dB**\n- **ΔMAE = +0.00203**\n- **ΔSSIM = +0.04099**\n- **PSNR / MAE / SSIM win rate = 100% / 100% / 100%**\n\nThis is the strongest single conclusion from the Mayo experiment:\n\n> **A model trained using synthetic supervision produces consistent, measurable improvement on every paired Mayo case in PSNR, MAE, and SSIM.**\n\nThat is a very strong transfer result, especially because this benchmark is external to the training pipeline.\n\n---\n\n### 4. V-Ultimate also outperforms Gaussian in mean PSNR, but that is not the main point\n\nGaussian is a strong classical baseline here:\n\n- **Gaussian mean PSNR = 41.18 dB**\n- **ΔPSNR vs Quarter = +1.69 dB**\n- **MAE and SSIM also improve substantially**\n- **PSNR win rate vs Quarter = 90%**\n\nSo Gaussian is clearly useful. This is important, because it shows that the external Mayo benchmark is not trivial: simple smoothing already helps.\n\nHowever, V-Ultimate still achieves:\n\n- **higher mean PSNR than Gaussian**\n- **larger PSNR gain over Quarter**\n- **perfect 10/10 PSNR win rate over Quarter**\n- **lower mean MAE than Gaussian**\n\nSo the result is not merely that “deep learning beats a weak baseline.”  \nRather, the more precise statement is:\n\n> **Even against a strong and relevant smoothing baseline, V-Ultimate achieves the best overall quantitative fidelity on real external low-dose CT.**\n\nThat said, the difference between V-Ultimate and Gaussian should still be interpreted carefully, because the real conceptual contribution here is not just superiority over Gaussian, but successful **external transfer from synthetic supervision to real paired data**.\n\n---\n\n### 5. The most important qualitative finding: V-Ultimate is more conservative than Gaussian\n\nA particularly important result appears in the change-versus-input statistics:\n\n- **Gaussian mean absolute change vs Quarter = 0.00690**\n- **V-Ultimate mean absolute change vs Quarter = 0.00316**\n\nand\n\n- **Gaussian max absolute change vs Quarter = 0.15326**\n- **V-Ultimate max absolute change vs Quarter = 0.02463**\n\nThis is one of the most important findings in the entire external evaluation.\n\nIt means that:\n\n> **V-Ultimate achieves higher mean PSNR than Gaussian while making much smaller modifications to the original quarter-dose image.**\n\nIn other words, the model is not simply “winning by editing more aggressively.”  \nIt is actually:\n\n- changing the image **less** on average,\n- changing the image **far less** in its worst-case magnitude,\n- and still obtaining **better reconstruction fidelity**.\n\nThis strongly supports the project’s **controlled restoration / do-no-harm** framing.\n\nA useful way to summarize this is:\n\n> **Gaussian is a stronger global smoother, but V-Ultimate is a more selective and conservative restorer.**\n\n---\n\n### 6. The gain is especially strong in soft-kernel cases, while sharp-kernel cases are more mixed\n\nAt the case level, a pattern is visible:\n\n#### Soft kernel (B30) cases\nV-Ultimate consistently produces large gains, for example:\n\n- L096: **43.68 → 45.69**\n- L333: **42.88 → 44.63**\n- L192: **44.20 → 45.77**\n- L143: **42.32 → 44.46**\n\nThese are strong improvements over both Quarter and Gaussian.\n\n#### Sharp kernel (D45) cases\nThe picture is more mixed:\n\n- in some sharp-kernel cases, V-Ultimate clearly improves over Quarter and slightly exceeds Gaussian;\n- in others, Gaussian is very close or slightly better.\n\nThis is a meaningful pattern rather than a weakness to hide.\n\nA reasonable interpretation is:\n\n> **The model transfers particularly well to the smoother soft-kernel reconstruction regime, while the sharper reconstruction regime remains somewhat harder and more variable.**\n\nThis is plausible, because sharp-kernel reconstruction tends to preserve more high-frequency texture and noise simultaneously, which makes restoration less straightforward.\n\n---\n\n### 7. Bootstrap confidence intervals confirm the gains are statistically stable\n\nThe bootstrap results strengthen the conclusion:\n\n#### V-Ultimate vs Quarter\n- **ΔPSNR = +2.35 dB**\n- **95% CI: [2.03, 2.67]**\n- **ΔMAE = +0.00203**\n- **95% CI: [0.00138, 0.00276]**\n- **ΔSSIM = +0.04099**\n- **95% CI: [0.02542, 0.05881]**\n\n#### Gaussian vs Quarter\n- **ΔPSNR = +1.69 dB**\n- **95% CI: [0.86, 2.54]**\n\nThe key point is that the confidence intervals for V-Ultimate’s improvement over Quarter are **consistently positive**, which supports the conclusion that the transfer effect is not a coincidence driven by a few favorable cases.\n\nSo the correct interpretation is:\n\n> **The external Mayo improvement is not only visible at the case level, but also statistically stable under paired bootstrap resampling.**\n\n---\n\n### 8. What does this mean for the overall project claim?\n\nThis Mayo experiment is arguably the strongest evidence that the synthetic supervision strategy is genuinely meaningful.\n\nWhy? Because the model was **not** trained on real paired Mayo data.  \nYet it still improves real external quarter-dose CT against matched full-dose reference.\n\nThat directly supports the central methodological claim:\n\n> **Physics-based synthetic supervision can train a restoration model that transfers to real clinical low-dose CT.**\n\nThis is more important than any single comparison against Gaussian.\n\nThe main contribution demonstrated here is not simply:\n\n> “Our model has higher PSNR than a classical filter.”\n\nIt is:\n\n> **A controlled restoration model trained with synthetic supervision can generalize beyond synthetic data and improve real external paired low-dose CT.**\n\n---\n\n### 9. Final interpretation of this section\n\nThe Mayo external paired evaluation supports four key conclusions:\n\n1. **V-Ultimate consistently improves real quarter-dose CT on every evaluated case.**\n2. **It achieves the best overall quantitative fidelity among the tested methods.**\n3. **It does so while making much smaller changes to the input than Gaussian, which supports the controlled-restoration claim.**\n4. **Most importantly, it demonstrates successful sim-to-real transfer from physics-based synthetic supervision to real external clinical data.**\n\nTherefore, the most accurate summary is not just:\n\n> “V-Ultimate beats Gaussian on Mayo.”\n\nIt is:\n\n> **V-Ultimate shows that synthetic supervision is not merely a synthetic-domain trick; it can train a restoration model that transfers to real external low-dose CT, improves fidelity, and does so in a more conservative and controlled way than classical smoothing.**","metadata":{}},{"cell_type":"code","source":"# Sanity check: verify changed voxel counts under different thresholds\ntest_case = selected_manifest[0]\n\nq_hu, f_hu, q01, f01 = load_paired_case_window(\n    test_case[\"q_files\"], test_case[\"f_files\"], test_case[\"start_idx\"],\n    depth=VOL_DEPTH, out_hw=(TARGET_H, TARGET_W)\n)\nrec01 = run_vultimate_mayo(q01)\n\nchange_map = np.abs(rec01 - q01).astype(np.float32)\n\nprint(\"case_id:\", test_case[\"case_id\"])\nprint(\"mean_abs_change:\", float(change_map.mean()))\nprint(\"max_abs_change :\", float(change_map.max()))\n\nfor thr in [0.05, 0.02, 0.01, 0.005, 0.002, 0.001]:\n    n = int((change_map >= thr).sum())\n    print(f\"thr={thr:0.3f} -> changed_voxels={n}\")\n\ndel q_hu, f_hu, q01, f01, rec01, change_map\ngc.collect()\nif torch.cuda.is_available():\n    torch.cuda.empty_cache()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T07:35:48.259079Z","iopub.execute_input":"2026-03-15T07:35:48.259817Z","iopub.status.idle":"2026-03-15T07:35:52.120601Z","shell.execute_reply.started":"2026-03-15T07:35:48.259784Z","shell.execute_reply":"2026-03-15T07:35:52.119925Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# TotalSegmentator on Mayo external paired evaluation (CLEAN / REUSABLE / LOW-PRINT)\n# Purpose:\n#   1) Reuse the same matched Mayo series from the Mayo cell\n#   2) Reuse existing TotalSegmentator outputs when available\n#   3) Quantify downstream anatomy benefit (Dice vs Full)\n#   4) Quantify anatomical overlap of Quarter->Ultimate edits\n#\n# Depends on the Mayo cell having already run:\n#   - mayo_case_manifest\n#   - load_paired_case_window\n#   - run_vultimate_mayo\n#   - hu01_to_hu\n#   - VOL_DEPTH, TARGET_H, TARGET_W\n#   - TOTALSEG_OUTDIR\n# ============================================================\n\nimport os, sys, gc, time, shutil, subprocess\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\ntry:\n    import nibabel as nib\nexcept Exception:\n    subprocess.run([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"nibabel\"], check=False)\n    import nibabel as nib\n\ntry:\n    from IPython.display import display\nexcept Exception:\n    display = print\n\n# ------------------------------------------------------------\n# 0) Checks\n# ------------------------------------------------------------\nrequired_globals = [\n    \"mayo_case_manifest\",\n    \"TOTALSEG_OUTDIR\",\n    \"load_paired_case_window\",\n    \"run_vultimate_mayo\",\n    \"hu01_to_hu\",\n    \"VOL_DEPTH\", \"TARGET_H\", \"TARGET_W\",\n]\nmissing = [k for k in required_globals if k not in globals()]\nif missing:\n    raise RuntimeError(f\"Run the Mayo cell first. Missing: {missing}\")\n\nOUTDIR = Path(TOTALSEG_OUTDIR)\nOUTDIR.mkdir(parents=True, exist_ok=True)\n\n# ------------------------------------------------------------\n# 1) Config\n# ------------------------------------------------------------\nN_TOTALSEG_CASES = min(int(globals().get(\"N_TOTALSEG_CASES\", 10)), len(mayo_case_manifest))\nTOTALSEG_CASE_IDS = globals().get(\"TOTALSEG_CASE_IDS\", None)   # optional, e.g. [1,2,3]\nTOTALSEG_TASK = str(globals().get(\"TOTALSEG_TASK\", \"total\"))\nTOTALSEG_CHANGE_THR = 0.01\nTOTALSEG_FAST = bool(globals().get(\"TOTALSEG_FAST\", True))\nTOTALSEG_REUSE_EXISTING = bool(globals().get(\"TOTALSEG_REUSE_EXISTING\", True))\nTOTALSEG_FORCE_CPU = bool(globals().get(\"TOTALSEG_FORCE_CPU\", True))\nTOTALSEG_SHOW_TOP_K = int(globals().get(\"TOTALSEG_SHOW_TOP_K\", 20))\n\n# ------------------------------------------------------------\n# 2) Select cases\n# ------------------------------------------------------------\nif TOTALSEG_CASE_IDS is not None:\n    selected_manifest = [\n        m for m in mayo_case_manifest\n        if int(m[\"case_id\"]) in set(map(int, TOTALSEG_CASE_IDS))\n    ]\nelse:\n    selected_manifest = mayo_case_manifest[:N_TOTALSEG_CASES]\n\nif len(selected_manifest) == 0:\n    raise RuntimeError(\"No Mayo cases selected for TotalSegmentator.\")\n\n# ------------------------------------------------------------\n# 3) Find / install TotalSegmentator CLI\n# ------------------------------------------------------------\ncli = None\nfor c in [\"TotalSegmentator\", \"totalsegmentator\"]:\n    if shutil.which(c):\n        cli = c\n        break\n\nif cli is None:\n    subprocess.run([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"TotalSegmentator\"], check=False)\n    for c in [\"TotalSegmentator\", \"totalsegmentator\"]:\n        if shutil.which(c):\n            cli = c\n            break\n\nif cli is None:\n    raise RuntimeError(\"TotalSegmentator CLI not available.\")\n\n# ------------------------------------------------------------\n# 4) Helpers\n# ------------------------------------------------------------\ndef save_nifti_hu(vol_hu_dhw, out_path):\n    affine = np.eye(4, dtype=np.float32)\n    vol_hwd = np.transpose(vol_hu_dhw.astype(np.float32), (1, 2, 0))\n    nib.save(nib.Nifti1Image(vol_hwd, affine), str(out_path))\n\ndef run_totalseg(in_nii, out_dir):\n    out_dir = Path(out_dir)\n    out_dir.mkdir(parents=True, exist_ok=True)\n\n    existing = list(out_dir.rglob(\"*.nii.gz\"))\n    if TOTALSEG_REUSE_EXISTING and len(existing) > 0:\n        return 0, \"cached\"\n\n    cmd = [cli, \"-i\", str(in_nii), \"-o\", str(out_dir), \"--task\", TOTALSEG_TASK]\n    if TOTALSEG_FAST:\n        cmd.append(\"--fast\")\n\n    env = os.environ.copy()\n    if TOTALSEG_FORCE_CPU:\n        env[\"CUDA_VISIBLE_DEVICES\"] = \"\"\n\n    p = subprocess.run(\n        cmd,\n        stdout=subprocess.PIPE,\n        stderr=subprocess.STDOUT,\n        text=True,\n        env=env,\n    )\n    return p.returncode, p.stdout\n\ndef dice_bin(a, b):\n    a = (a > 0.5)\n    b = (b > 0.5)\n    denom = int(a.sum()) + int(b.sum())\n    if denom == 0:\n        return 1.0\n    inter = int((a & b).sum())\n    return 2.0 * inter / denom\n\ndef resolve_mask(seg_dir, filename):\n    p = Path(seg_dir) / filename\n    if p.exists():\n        return p\n    cands = list(Path(seg_dir).rglob(filename))\n    return cands[0] if len(cands) else None\n\n# ------------------------------------------------------------\n# 5) Main loop\n# ------------------------------------------------------------\noverlap_rows = []\ndice_rows = []\ncase_rows = []\n\nt0_all = time.time()\n\nfor j, case_info in enumerate(selected_manifest, 1):\n    cid = int(case_info[\"case_id\"])\n    rel_dir = case_info[\"rel_dir\"]\n    st = int(case_info[\"start_idx\"])\n\n    print(f\"[{j:02d}/{len(selected_manifest)}] case_id={cid:02d}\")\n\n    q_hu, f_hu, q01, f01 = load_paired_case_window(\n        case_info[\"q_files\"], case_info[\"f_files\"], st,\n        depth=VOL_DEPTH, out_hw=(TARGET_H, TARGET_W)\n    )\n    rec01 = run_vultimate_mayo(q01)\n\n    # Quarter -> Ultimate change map\n    change_map = np.abs(rec01 - q01).astype(np.float32)\n    change_mask = (change_map >= TOTALSEG_CHANGE_THR).astype(np.uint8)\n    total_changed_vox = int(change_mask.sum())\n\n    # Save NIfTI\n    safe_rel = rel_dir.replace(\"/\", \"__\").replace(\"\\\\\", \"__\")\n    case_dir = OUTDIR / f\"case_{cid:02d}_{safe_rel}_st{st:04d}\"\n    case_dir.mkdir(parents=True, exist_ok=True)\n\n    nii_map = {\n        \"full\": case_dir / \"full.nii.gz\",\n        \"quarter\": case_dir / \"quarter.nii.gz\",\n        \"ultimate\": case_dir / \"ultimate.nii.gz\",\n    }\n    save_nifti_hu(f_hu, nii_map[\"full\"])\n    save_nifti_hu(q_hu, nii_map[\"quarter\"])\n    save_nifti_hu(hu01_to_hu(rec01), nii_map[\"ultimate\"])\n\n    # Run / reuse TotalSeg\n    seg_dirs = {}\n    failed = False\n    for nm, nii_p in nii_map.items():\n        rc, outtxt = run_totalseg(nii_p, case_dir / f\"seg_{nm}\")\n        if rc != 0:\n            failed = True\n            seg_dirs[nm] = None\n            print(f\"   {nm}: FAILED\")\n            tail = \"\\n\".join(str(outtxt).splitlines()[-8:])\n            if tail:\n                print(tail)\n        else:\n            seg_dirs[nm] = case_dir / f\"seg_{nm}\"\n\n    if failed or seg_dirs[\"full\"] is None or seg_dirs[\"quarter\"] is None or seg_dirs[\"ultimate\"] is None:\n        print(\"   skipped (missing segmentation outputs)\")\n        del q_hu, f_hu, q01, f01, rec01, change_map, change_mask\n        gc.collect()\n        if torch.cuda.is_available():\n            torch.cuda.empty_cache()\n        continue\n\n    # overlap + dice\n    full_masks = sorted(list(Path(seg_dirs[\"full\"]).rglob(\"*.nii.gz\")))\n    change_mask_hwd = np.transpose(change_mask, (1, 2, 0))\n\n    best_mask = None\n    best_share = -1.0\n    nz_masks = 0\n\n    for mf in full_masks:\n        try:\n            mask_name = mf.name.replace(\".nii.gz\", \"\")\n            full_arr = nib.load(str(mf)).get_fdata()\n            full_bin = (full_arr > 0.5).astype(np.uint8)\n            seg_vox = int(full_bin.sum())\n            if seg_vox == 0:\n                continue\n\n            inter = int((full_bin * change_mask_hwd).sum())\n            if inter > 0:\n                nz_masks += 1\n\n            change_in_seg_ratio = (inter / seg_vox) if seg_vox > 0 else np.nan\n            seg_share_of_change = (inter / total_changed_vox) if total_changed_vox > 0 else np.nan\n\n            if not np.isnan(seg_share_of_change) and seg_share_of_change > best_share:\n                best_share = seg_share_of_change\n                best_mask = mask_name\n\n            overlap_rows.append({\n                \"case_id\": cid,\n                \"rel_dir\": rel_dir,\n                \"start_idx\": st,\n                \"mask_name\": mask_name,\n                \"seg_vox_full\": seg_vox,\n                \"changed_vox_total\": total_changed_vox,\n                \"intersect_change_vox\": inter,\n                \"change_in_seg_ratio\": change_in_seg_ratio,\n                \"seg_share_of_change\": seg_share_of_change,\n            })\n\n            qmask_path = resolve_mask(seg_dirs[\"quarter\"], mf.name)\n            umask_path = resolve_mask(seg_dirs[\"ultimate\"], mf.name)\n\n            if qmask_path is not None and qmask_path.exists():\n                q_arr = nib.load(str(qmask_path)).get_fdata()\n                dice_rows.append({\n                    \"case_id\": cid,\n                    \"rel_dir\": rel_dir,\n                    \"start_idx\": st,\n                    \"mask_name\": mask_name,\n                    \"variant\": \"quarter\",\n                    \"dice_vs_full\": float(dice_bin(full_arr, q_arr)),\n                })\n\n            if umask_path is not None and umask_path.exists():\n                u_arr = nib.load(str(umask_path)).get_fdata()\n                dice_rows.append({\n                    \"case_id\": cid,\n                    \"rel_dir\": rel_dir,\n                    \"start_idx\": st,\n                    \"mask_name\": mask_name,\n                    \"variant\": \"ultimate\",\n                    \"dice_vs_full\": float(dice_bin(full_arr, u_arr)),\n                })\n\n        except Exception:\n            continue\n\n    case_rows.append({\n        \"case_id\": cid,\n        \"rel_dir\": rel_dir,\n        \"start_idx\": st,\n        \"changed_vox_total\": total_changed_vox,\n        \"n_masks_total\": len(full_masks),\n        \"n_masks_nonzero\": nz_masks,\n        \"top_mask\": best_mask,\n        \"top_mask_seg_share\": best_share if best_share >= 0 else np.nan,\n        \"mean_abs_change_q_to_u\": float(change_map.mean()),\n        \"max_abs_change_q_to_u\": float(change_map.max()),\n    })\n\n    del q_hu, f_hu, q01, f01, rec01, change_map, change_mask, change_mask_hwd\n    gc.collect()\n    if torch.cuda.is_available():\n        torch.cuda.empty_cache()\n\n# ------------------------------------------------------------\n# 6) Save outputs\n# ------------------------------------------------------------\nif len(overlap_rows) == 0:\n    raise RuntimeError(\"No TotalSeg overlap rows generated.\")\n\ndf_case = pd.DataFrame(case_rows)\ndf_case.to_csv(OUTDIR / \"totalseg_case_summary_mayo10.csv\", index=False)\n\ndf_ov = pd.DataFrame(overlap_rows)\ndf_ov.to_csv(OUTDIR / \"totalseg_overlap_raw_mayo10.csv\", index=False)\n\ndf_d = pd.DataFrame(dice_rows)\ndf_d.to_csv(OUTDIR / \"totalseg_dice_raw_mayo10.csv\", index=False)\n\nov_sum = (\n    df_ov.groupby(\"mask_name\", as_index=False)\n    .agg(\n        n_case=(\"case_id\", \"nunique\"),\n        n_case_nonzero=(\"intersect_change_vox\", lambda s: int((np.asarray(s) > 0).sum())),\n        pct_cases_nonzero=(\"intersect_change_vox\", lambda s: float((np.asarray(s) > 0).mean() * 100.0)),\n        mean_change_in_seg_ratio=(\"change_in_seg_ratio\", \"mean\"),\n        median_change_in_seg_ratio=(\"change_in_seg_ratio\", \"median\"),\n        mean_seg_share_of_change=(\"seg_share_of_change\", \"mean\"),\n        median_seg_share_of_change=(\"seg_share_of_change\", \"median\"),\n        max_seg_share_of_change=(\"seg_share_of_change\", \"max\"),\n    )\n    .sort_values([\"n_case_nonzero\", \"mean_seg_share_of_change\"], ascending=[False, False])\n    .reset_index(drop=True)\n)\nov_sum.to_csv(OUTDIR / \"totalseg_overlap_summary_mayo10.csv\", index=False)\n\nif len(df_d) == 0:\n    raise RuntimeError(\"No TotalSeg dice rows generated.\")\n\nd_sum = (\n    df_d.groupby(\"variant\", as_index=False)\n    .agg(\n        mean_dice_vs_full=(\"dice_vs_full\", \"mean\"),\n        std_dice_vs_full=(\"dice_vs_full\", \"std\"),\n        median_dice_vs_full=(\"dice_vs_full\", \"median\"),\n        n_rows=(\"dice_vs_full\", \"size\"),\n    )\n    .sort_values(\"mean_dice_vs_full\", ascending=False)\n    .reset_index(drop=True)\n)\nd_sum.to_csv(OUTDIR / \"totalseg_dice_summary_mayo10.csv\", index=False)\n\nqd = df_d[df_d[\"variant\"] == \"quarter\"][[\"case_id\", \"mask_name\", \"dice_vs_full\"]].rename(\n    columns={\"dice_vs_full\": \"dice_q\"}\n)\nud = df_d[df_d[\"variant\"] == \"ultimate\"][[\"case_id\", \"mask_name\", \"dice_vs_full\"]].rename(\n    columns={\"dice_vs_full\": \"dice_u\"}\n)\ndm = qd.merge(ud, on=[\"case_id\", \"mask_name\"], how=\"inner\")\ndm[\"delta_u_minus_q\"] = dm[\"dice_u\"] - dm[\"dice_q\"]\n\nd_mask = (\n    dm.groupby(\"mask_name\", as_index=False)\n    .agg(\n        n_case=(\"case_id\", \"nunique\"),\n        mean_delta_u_minus_q=(\"delta_u_minus_q\", \"mean\"),\n        median_delta_u_minus_q=(\"delta_u_minus_q\", \"median\"),\n        pct_cases_improved=(\"delta_u_minus_q\", lambda s: float((np.asarray(s) > 0).mean() * 100.0)),\n    )\n    .sort_values([\"mean_delta_u_minus_q\", \"pct_cases_improved\"], ascending=[False, False])\n    .reset_index(drop=True)\n)\nd_mask.to_csv(OUTDIR / \"totalseg_dice_delta_by_mask_mayo10.csv\", index=False)\n\n# ------------------------------------------------------------\n# 7) Compact display\n# ------------------------------------------------------------\nelapsed_min = (time.time() - t0_all) / 60.0\nprint(f\"\\nDone in {elapsed_min:.1f} min | change_thr={TOTALSEG_CHANGE_THR}\")\n\nprint(\"\\n=== Case summary ===\")\ndisplay(df_case)\n\nprint(\"\\n=== Dice summary ===\")\ndisplay(d_sum)\n\nprint(f\"\\n=== Top-{TOTALSEG_SHOW_TOP_K} overlap summary ===\")\ndisplay(ov_sum.head(TOTALSEG_SHOW_TOP_K))\n\nprint(f\"\\n=== Top-{TOTALSEG_SHOW_TOP_K} dice improvement by mask ===\")\ndisplay(d_mask.head(TOTALSEG_SHOW_TOP_K))\n\nprint(\"\\nSaved:\")\nfor p in sorted(OUTDIR.glob(\"*mayo10*.csv\")):\n    print(\" -\", p)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T07:36:13.050864Z","iopub.execute_input":"2026-03-15T07:36:13.051216Z","iopub.status.idle":"2026-03-15T07:44:36.787082Z","shell.execute_reply.started":"2026-03-15T07:36:13.051186Z","shell.execute_reply":"2026-03-15T07:44:36.786442Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Results Interpretation — TotalSegmentator Anatomy Analysis\n\n### 1. What does this experiment test?\n\nThis section asks a different question from CRM and Mayo.\n\nCRM shows whether the model can recover clinically relevant signal.  \nMayo shows whether synthetic-supervision training transfers to real external low-dose CT.  \nThe TotalSegmentator analysis asks:\n\n> **Where is the model actually editing the image, and are those edits anatomically meaningful rather than globally smoothing the whole scan?**\n\nTo answer that, we used the restored Mayo volumes and examined:\n\n- **change localization**: where V-Ultimate differs from the original quarter-dose input,\n- **anatomical overlap**: which organs and structures those changes fall into,\n- **segmentation consistency**: whether the restored image becomes easier to segment relative to the full-dose reference.\n\nThis makes TotalSeg an **anatomy-aware interpretability analysis**, not just another image-quality benchmark.\n\n---\n\n### 2. The key result: very small voxelwise edits produce a large segmentation benefit\n\nThe most important global result is the Dice summary:\n\n- **Quarter mean Dice vs Full = 0.6137**\n- **V-Ultimate mean Dice vs Full = 0.7395**\n- **Absolute Dice improvement = +0.1258**\n\nThis is a large gain.\n\nAt the same time, the per-case intensity change statistics remain modest:\n\n- **mean absolute change (Quarter → Ultimate)** is only about  \n  - **0.0015–0.0023** in the soft-kernel cases,\n  - **0.0037–0.0065** in the sharp-kernel cases.\n- **max absolute change** stays around **0.020–0.026**.\n\nSo the strongest first conclusion is:\n\n> **V-Ultimate does not need large image rewrites to create a substantial anatomy-level improvement.**\n\nThis is exactly the pattern expected from a controlled restoration model:\nsmall but targeted edits can have a large downstream structural effect.\n\n---\n\n### 3. The edits are not globally uniform; they are anatomically concentrated\n\nThe overlap summary shows that changes consistently fall into specific anatomical structures rather than spreading uniformly across the whole scan.\n\nThe most frequently involved masks are:\n\n- **small bowel**\n- **colon**\n- **autochthonous back muscles**\n- **liver**\n- **heart**\n- **iliopsoas**\n- **kidney**\n- **aorta**\n- **vertebrae**\n\nSeveral of these appear in **100% of applicable cases**, with nonzero overlap between the model’s change map and the segmentation mask.\n\nThis supports the interpretation that:\n\n> **The model is not merely denoising the image everywhere equally; it is modifying regions that correspond to real anatomical structures, especially boundary-rich and heterogeneous tissues.**\n\nThat is an important qualitative distinction.  \nA purely global smoother would be expected to reduce noise broadly, but not necessarily concentrate its effect in anatomically meaningful areas.\n\n---\n\n### 4. Bowel, liver, and muscular structures dominate the change map\n\nThe top overlap structures are especially informative:\n\n- **small bowel**\n- **colon**\n- **autochthon_right**\n- **liver**\n- **iliopsoas_right**\n- **heart**\n\nThese are all structures that tend to have:\n\n- complex edges,\n- heterogeneous internal texture,\n- partial-volume effects,\n- and sensitivity to reconstruction noise.\n\nThis suggests that V-Ultimate is especially active in regions where low-dose degradation most strongly interferes with anatomical delineation.\n\nA useful interpretation is:\n\n> **The model’s edits cluster in anatomically difficult, boundary-sensitive regions, rather than in homogeneous background or irrelevant areas.**\n\nThat is consistent with a restoration mechanism that is doing more than simple smoothing.\n\n---\n\n### 5. Sharp-kernel cases require much more editing than soft-kernel cases\n\nThe case summary shows a striking pattern:\n\n- **Sharp-kernel (D45) cases** have very large `changed_vox_total`\n- **Soft-kernel (B30) cases** have much smaller `changed_vox_total`\n\nFor example:\n\n- Sharp kernel cases: millions of changed voxels\n- Soft kernel cases: tens to hundreds of thousands\n\nAt the same time:\n\n- **mean absolute change remains small in both groups**\n- but is consistently **larger in sharp-kernel cases**\n\nThis is a very sensible result.\n\nSharp-kernel reconstructions preserve more high-frequency structure, but they also amplify noise and edge instability. That means the restoration model needs to intervene more broadly to regularize the image while preserving structure. Soft-kernel cases are already smoother, so fewer corrections are needed.\n\nSo this section supports another useful interpretation:\n\n> **The model adapts its editing behavior to reconstruction regime: it edits much more in sharp-kernel data, where the restoration problem is harder, and much less in soft-kernel data, where the input is already more stable.**\n\nThat is a sign of conditional, context-sensitive restoration rather than one fixed smoothing behavior.\n\n---\n\n### 6. Liver and bowel are especially important because they account for both overlap frequency and change share\n\nSome structures appear often, but only occupy a tiny fraction of the total change volume.  \nOthers contribute meaningfully to both:\n\n- how often they are touched,\n- and how much of the total change map they account for.\n\nTwo especially notable examples are:\n\n- **liver**\n- **small bowel / colon**\n\nFor instance, the liver reaches very large `seg_share_of_change` values in some cases, especially in the L143 examples, meaning a substantial fraction of the model’s total edits fall inside that organ. Bowel structures also appear consistently across all cases.\n\nThis suggests that:\n\n> **The model’s corrections are not only frequent in these structures, but sometimes dominant there, indicating that large portions of the restoration effort are being spent in anatomically meaningful regions.**\n\nThis is much more informative than simply saying “the model changes a lot of voxels.”\n\n---\n\n### 7. The strongest segmentation improvements occur in boundary-sensitive and fine structures\n\nThe mask-wise Dice improvement table is one of the most revealing parts of the analysis.\n\nThe largest gains appear in structures such as:\n\n- **gluteal muscles**\n- **vertebrae**\n- **ribs**\n- **inferior vena cava**\n- **pancreas**\n- **esophagus**\n- **iliopsoas**\n- **costal cartilages**\n- **portal / splenic veins**\n\nThese are structures that are often difficult to segment under low-dose degradation because they depend heavily on:\n\n- clean edge definition,\n- preservation of thin boundaries,\n- local contrast stability,\n- and suppression of noisy texture confusion.\n\nThis pattern strongly supports the interpretation that:\n\n> **V-Ultimate is improving anatomical usability at the level of fine structures and boundaries, not merely making images look smoother.**\n\nThat is a much stronger claim than a generic denoising statement.\n\n---\n\n### 8. Some structures show high overlap ratio but low total change share — this is not contradictory\n\nA few masks, such as:\n\n- **spinal cord**\n- **aorta**\n- **vertebrae**\n\nshow relatively high `change_in_seg_ratio` but low `seg_share_of_change`.\n\nThis is not a contradiction. It simply means:\n\n- a substantial portion of that structure is edited,\n- but the structure itself is relatively small compared with the full body volume.\n\nSo these results should be interpreted as:\n\n> **The model is highly active within certain compact, clinically meaningful structures, even if those structures do not dominate the total edited voxel count.**\n\nThis is actually another sign of targeted behavior.\n\n---\n\n### 9. The threshold matters: `change_thr = 0.01` is permissive, yet the edits still remain structured\n\nThis analysis used:\n\n- **change threshold = 0.01**\n\nThat is a relatively permissive threshold, meaning even fairly small voxelwise changes are counted as edits. Under such a threshold, a purely indiscriminate method might appear to modify everything.\n\nBut even with this permissive setting, the overlap results are still anatomically organized and the downstream Dice improves strongly.\n\nThat makes the result more convincing:\n\n> **Even when small edits are counted, the change map still aligns with anatomical structures rather than dissolving into uniform global modification.**\n\n---\n\n### 10. What does this mean for the overall project claim?\n\nThis TotalSeg result does not primarily claim that V-Ultimate is “better at segmentation” as an end goal.  \nInstead, it provides anatomical evidence for **how** the restoration is behaving.\n\nThe overall interpretation is:\n\n1. **The model makes relatively small voxelwise edits;**\n2. **those edits are concentrated in real anatomical structures;**\n3. **the strongest downstream gains appear in boundary-sensitive tissues and fine structures;**\n4. **and the restored image becomes substantially more consistent with the full-dose reference at the segmentation level.**\n\nSo the most important conclusion is:\n\n> **V-Ultimate is not acting like a blind global smoother. It is performing anatomically anchored restoration.**\n\nThat is exactly the kind of evidence needed to support the broader claim that the model is doing controlled, meaningful recovery rather than arbitrary image manipulation.\n\n---\n\n### 11. Final interpretation of this section\n\nThe TotalSegmentator analysis strengthens the project in three major ways:\n\n- **Interpretability:** it shows where the model is editing;\n- **Anatomical plausibility:** it shows that edits align with real organs and structures;\n- **Functional relevance:** it shows that these edits improve segmentation consistency relative to full-dose CT.\n\nTherefore, the best summary of this section is:\n\n> **V-Ultimate produces small but anatomically structured edits, and those edits translate into a large downstream segmentation improvement. This supports the claim that the model performs controlled, anatomy-aware restoration rather than simple global smoothing.**","metadata":{}},{"cell_type":"markdown","source":"## Future Work\n\nWhile the current notebook already supports the main methodological claim — that physics-grounded synthetic supervision can train a clinically useful, externally transferable, and anatomically structured CT restoration model — several important next steps remain.\n\n### 1. Scaling behavior of synthetic supervision\nA natural next question is whether the benefit of synthetic supervision continues to grow as more training cases are added.\n\nFuture work should evaluate checkpoints trained with different numbers of synthetic training UIDs (for example: 10, 30, 60, 100, 200, 300) and measure how performance changes on:\n\n- **strict OOD CRM**\n- **external paired Mayo CT**\n- **anatomy-aware downstream metrics such as TotalSegmentator Dice**\n\nThis would answer a higher-level methodological question:\n\n> Does physics-grounded synthetic supervision merely work, or does it exhibit a meaningful scaling law?\n\nA clear positive trend would strengthen the argument that the synthetic supervision contains real training value rather than only a one-time benefit.\n\n---\n\n### 2. Stronger control in neutral and ambiguous cases\nThe current results show that V-Ultimate is strongest in **harmed-positive** cases, where degradation truly suppresses clinically relevant signal.  \nHowever, **neutral** and **overcall-like** cases remain more fragile and noise-sensitive.\n\nFuture work should therefore focus on improving the model’s ability to decide:\n\n- when restoration is necessary,\n- when restoration should be minimal,\n- and when the safest action is to leave the input largely unchanged.\n\nThis could be explored through stronger edit-gating, uncertainty-aware penalties, or more explicit “do-not-edit” supervision.\n\n---\n\n### 3. Broader external validation\nThe Mayo experiment provides strong external paired evidence, but it is still limited in size and protocol diversity.\n\nFuture work should extend external evaluation to:\n\n- more patients,\n- more scanner protocols,\n- more reconstruction kernels,\n- and ideally additional institutions.\n\nThis would test whether the current sim-to-real transfer result remains stable across broader clinical variation.\n\n---\n\n### 4. Richer anatomy-aware downstream evaluation\nThe current TotalSegmentator analysis shows that the edits are small but anatomically structured.  \nA next step would be to evaluate whether restoration consistently improves downstream performance in additional tasks, such as:\n\n- lesion localization,\n- vessel analysis,\n- organ-specific segmentation,\n- or detection/classification systems beyond the current clinical proxy model.\n\nThis would help clarify whether the restoration is useful across multiple clinically relevant endpoints, not just one.\n\n---\n\n### 5. Better characterization of noise sensitivity\nThe Monte Carlo results suggest that the model is robust **in expectation**, but still noise-sensitive in a substantial fraction of cases.\n\nFuture work should analyze this more systematically by identifying:\n\n- which anatomy types are most noise-sensitive,\n- whether instability correlates with sharp kernels or borderline scores,\n- and whether specific degradation regimes are consistently harder.\n\nThis could help distinguish true model weakness from intrinsic ambiguity in the corrupted input.\n\n---\n\n### 6. Comparison with more modern restoration families\nThis notebook intentionally focuses on validating the synthetic-supervision methodology rather than maximizing architecture complexity.\n\nA useful future direction would be to test whether the same synthetic-supervision framework can train other restoration backbones, such as:\n\n- transformer-based restorers,\n- diffusion-style denoisers,\n- or uncertainty-aware reconstruction models.\n\nThat would help separate the contribution of the **training methodology** from the contribution of the current implementation vehicle.\n\n---\n\n## Summary\nThe present work supports the claim that physics-grounded synthetic supervision is a credible training signal under paired-data scarcity.  \nFuture work should now test how far this idea can scale, how safely it behaves in ambiguous cases, and how broadly it transfers across models, protocols, and downstream tasks.","metadata":{}},{"cell_type":"markdown","source":"# Judge Summary — Framing the Project as a Methodology, Not Just a Model\n\n## What did I actually build?\nI built and tested a **physics-grounded synthetic-supervision pipeline for CT restoration**.\n\nThe restoration model (**V-Ultimate**) is the implementation vehicle, but the broader contribution is **not only the network itself**. The broader contribution is the **training-and-validation methodology** for situations where true paired data are missing.\n\nIn simple terms:\n\n1. start from higher-quality CT,\n2. generate paired supervision through physics-grounded degradation,\n3. train a restoration model under conservative edit constraints,\n4. and then test whether that learned behavior remains useful on strict OOD data and real external low-dose CT.\n\n---\n\n## What is the scientific problem?\nIn medical imaging, we rarely have perfect **clean / degraded pairs from the exact same anatomy**. That makes supervised restoration fundamentally difficult.\n\nSo the real scientific question is not:\n\n> “Can I make a restoration network with better numbers?”\n\nIt is:\n\n> **When true paired data are unavailable, can physics-grounded synthetic degradations provide supervision that is reliable enough to train a clinically relevant and conservative CT restoration model?**\n\n---\n\n## What is my answer?\nMy answer is:\n\n> **Yes — but only if the synthetic supervision survives a multi-level evidence chain.**\n\nIn this project, that evidence chain is:\n\n1. **Clinical-proxy usefulness on strict held-out OOD data (CRM)**  \n   The model shows meaningful rescue behavior rather than merely improving visual appearance.\n\n2. **External transfer to real paired Mayo low-dose CT**  \n   The model improves real quarter-dose CT against matched full-dose reference.\n\n3. **Anatomy-aware structure (TotalSegmentator)**  \n   The edits are small but anatomically organized, rather than random global smoothing.\n\nA Monte Carlo stability analysis was also used as a robustness probe, showing that performance is positive in expectation but still noise-sensitive in some ambiguous cases.\n\n---\n\n## What is the key conclusion?\n> **My contribution is not just a restoration model. It is a methodology for constructing synthetic supervision from high-quality CT, then validating that supervision through clinical-proxy behavior, external transfer, and anatomy-aware analysis.**\n\nThis is the central point of the project.\n\n---\n\n## Why does that matter?\nBecause many medical AI problems are limited not by model ideas, but by the absence of reliable paired supervision.\n\nIf a synthetic-supervision pipeline can be shown to train models that are:\n\n- useful,\n- conservative,\n- externally transferable,\n- and anatomically meaningful,\n\nthen it becomes a practical way to train restoration systems in settings where direct paired ground truth is difficult or impossible to collect.\n\n---\n\n## Careful version of the takeaway\nThis project does **not** claim that the simulator is a perfect replica of every scanner protocol.\n\nIt makes a narrower and more defensible claim:\n\n> **The synthetic supervision is realistic enough to act as a useful training signal, because the resulting restoration policy behaves beneficially on strict OOD data, transfers to real external low-dose CT, and makes small but anatomically structured edits.**\n\n---\n\n# Judge Q&A — Likely Questions and Strong Answers\n\n### Q1. Is your contribution the restoration model, or the synthetic data methodology?\n**Answer:**  \nThe restoration model is the implementation vehicle, but the broader contribution is the methodology. The main question of the project is how to train restoration when real paired data are missing. The model matters because it tests whether the synthetic supervision actually works, but the deeper contribution is the pipeline for generating and validating that supervision.\n\n---\n\n### Q2. Isn’t generating degraded images easy? Why is that scientifically interesting?\n**Answer:**  \nGenerating degraded images is easy if the goal is only to make images look worse. The scientific difficulty is showing that the resulting synthetic pairs are useful enough to train a model that works on real data. My contribution is not “I can blur an image.” It is showing that physics-grounded degradation can serve as a credible training signal when evaluated through OOD behavior, external transfer, and anatomy-aware analysis.\n\n---\n\n### Q3. Did you prove that your simulator exactly matches real scanner physics?\n**Answer:**  \nNo, and that is not my claim. I do not claim distributional identity to every real scanner protocol. My claim is more practical and more defensible: the synthetic supervision is realistic enough to train a restoration policy that behaves usefully and conservatively on real external data.\n\n---\n\n### Q4. If the simulator is not perfect, why should anyone trust the training signal?\n**Answer:**  \nBecause I do not ask the reader to trust the simulator only on theory. I test the consequence of using it. If a model trained on synthetic pairs improves strict held-out OOD cases, transfers to external paired Mayo low-dose CT, and makes small but anatomically structured edits, then the supervision is empirically supported as a useful training signal.\n\n---\n\n### Q5. Why not just use Gaussian smoothing or another classical filter?\n**Answer:**  \nClassical filters are strong baselines, and I treat them seriously. In fact, part of the point is not to deny that smoothing can help. The difference is that my method aims for **controlled restoration**, not just generic denoising. On Mayo, for example, V-Ultimate achieved stronger reconstruction quality while making smaller average and maximum changes than Gaussian smoothing. That supports the claim that it is doing more selective, conservative restoration rather than simply editing more aggressively.\n\n---\n\n### Q6. Did your model beat every baseline on every metric?\n**Answer:**  \nNo, and I do not claim that. This project is not framed as “my model universally wins everything.” The stronger claim is that the synthetic-supervision pipeline produces a model with a good balance of usefulness, restraint, external transfer, and anatomical structure. That is a more meaningful methodological result than a one-metric win.\n\n---\n\n### Q7. If Bilateral or Gaussian is strong in some settings, why is your project still valuable?\n**Answer:**  \nBecause the project is not only about outperforming classical filters in every situation. It is about whether synthetic supervision can train a restoration policy that transfers beyond the synthetic domain. The strongest evidence for that is the external Mayo paired evaluation and the anatomy-aware TotalSegmentator analysis, not just one internal comparison table.\n\n---\n\n### Q8. Where is the model strongest, and where is it weakest?\n**Answer:**  \nIt is strongest in cases where degradation truly harms clinically relevant signal — the “harmed-positive” regime. That is where recovery matters most, and that is where the model shows its clearest benefit. It is weaker in more neutral or ambiguous cases, where the safest behavior is often to do very little. That limitation is important and honest, and it points directly to future work on stronger edit gating and uncertainty-aware restoration.\n\n---\n\n### Q9. What does the Mayo experiment prove that CRM does not?\n**Answer:**  \nCRM shows usefulness on strict held-out OOD data within the internal problem setting. Mayo is different: it is an **external paired real-data transfer test**. Because the model was trained with synthetic supervision but evaluated against real quarter-dose / full-dose pairs, Mayo is the strongest evidence that the learned restoration behavior transfers beyond the synthetic training setup.\n\n---\n\n### Q10. What does TotalSegmentator add beyond PSNR and SSIM?\n**Answer:**  \nPSNR and SSIM tell us whether the restored image is numerically closer to the full-dose reference. TotalSegmentator asks a different question: **where** is the model editing, and do those edits make anatomical structure more usable? In this project, the answer is yes: the edits are very small in magnitude, but they are anatomically concentrated and they improve segmentation Dice relative to full-dose reference.\n\n---\n\n### Q11. Why is “small edits” such an important point?\n**Answer:**  \nIn medical restoration, a method can look strong simply by rewriting the image aggressively. That is risky. A more clinically defensible behavior is to make the smallest edits needed to improve usefulness. In this project, the Mayo and TotalSeg results show that V-Ultimate achieves meaningful gains while changing the image much less than Gaussian smoothing. That supports the do-no-harm framing.\n\n---\n\n### Q12. Does Monte Carlo stability show the method is perfectly robust?\n**Answer:**  \nNo. The Monte Carlo analysis shows something more nuanced: the model is positive in expectation, especially in clinically important harmed-positive cases, but some individual cases are still sensitive to the exact stochastic corruption pattern. That does not invalidate the method; it identifies where stability is already good and where future safety improvements are still needed.\n\n---\n\n### Q13. What is the single most important scientific takeaway?\n**Answer:**  \nThe most important takeaway is that **physics-grounded synthetic supervision can be a credible training signal for CT restoration when real paired data are unavailable**, provided that the resulting model is judged not only by image quality, but also by clinical-proxy usefulness, external transfer, and anatomy-aware structure.\n\n---\n\n### Q14. If you had more time, what would you do next?\n**Answer:**  \nThe most important next step would be a **scaling study**: train or evaluate checkpoints across multiple synthetic training sizes and measure how performance changes on CRM, Mayo, and anatomy-aware downstream metrics. That would test whether the value of physics-grounded synthetic supervision continues to grow with more training cases, rather than only working at one scale.","metadata":{}},{"cell_type":"markdown","source":"# Science Fair abstract (short version)\n\nThis project studies a practical problem in medical AI: paired “clean / degraded” training data for CT restoration are usually unavailable. Instead of relying on arbitrary image corruption, I use **physics-grounded degradation synthesis** to generate paired supervision from higher-quality CT volumes. A restoration model is then trained under a **do-no-harm** principle so that improvements remain conservative rather than overly aggressive.\n\nThe main contribution is methodological: not merely a network architecture, but a framework for **creating and validating synthetic supervision**. The framework is tested through four layers of evidence: internal usefulness on controlled synthetic degradations, Monte Carlo stability across stochastic noise realizations, anatomy-overlap analysis using TotalSegmentator, and cross-domain evaluation on real paired Mayo low-dose CT. Together, these experiments support the conclusion that physics-based synthetic supervision can be trustworthy for training medical image restoration models even when true paired clinical data are missing.\n","metadata":{}},{"cell_type":"markdown","source":"# 90-second oral pitch\n\nMy project asks a simple question: **if hospitals do not have perfect paired training scans, can we still train restoration AI responsibly?**\n\nIn CT imaging, we usually do not have the exact same patient scanned once in a perfect state and once in a realistically degraded state. So instead of waiting for impossible labels, I use imaging physics to **simulate realistic degradations** from higher-quality CT data and create paired supervision for training.\n\nBut generating synthetic data is not enough. The important question is whether that supervision is actually trustworthy. So I validate it in four ways: first, whether the trained model helps on difficult held-out cases; second, whether the result is stable across different random noise realizations; third, whether the model’s edits are concentrated in anatomical regions instead of changing the whole scan randomly; and fourth, whether training on synthetic data transfers to **real paired Mayo low-dose CT**.\n\nSo the main contribution is not just one neural network. It is a **methodology for creating and validating synthetic supervision** when real paired medical data are missing.\n","metadata":{}},{"cell_type":"code","source":"!zip -r output.zip /kaggle/working","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T08:00:38.410136Z","iopub.execute_input":"2026-03-15T08:00:38.410993Z","iopub.status.idle":"2026-03-15T08:02:16.240448Z","shell.execute_reply.started":"2026-03-15T08:00:38.410949Z","shell.execute_reply":"2026-03-15T08:02:16.239638Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!echo '{\"title\":\"temp-output\",\"id\":\"andrewrenlin/temp-output\",\"licenses\":[{\"name\":\"CC0-1.0\"}]}' > dataset-metadata.json","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T08:02:31.920648Z","iopub.execute_input":"2026-03-15T08:02:31.921535Z","iopub.status.idle":"2026-03-15T08:02:32.104415Z","shell.execute_reply.started":"2026-03-15T08:02:31.92146Z","shell.execute_reply":"2026-03-15T08:02:32.103243Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!kaggle datasets create -p .","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-15T08:02:44.378314Z","iopub.execute_input":"2026-03-15T08:02:44.378684Z","iopub.status.idle":"2026-03-15T08:03:05.490217Z","shell.execute_reply.started":"2026-03-15T08:02:44.378634Z","shell.execute_reply":"2026-03-15T08:03:05.489425Z"}},"outputs":[],"execution_count":null}]}