{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# =====================================================\n# OSIC CT 共用前段预处理脚本\n# 目标：把每个病人一堆数量不等的 DICOM 切片，统一处理成\n#       一份「标准化、co-registered 的肺三维体积 + 肺掩膜」，\n#       存成中间文件，供【放射组学】和【CNN】两条线共用。\n#\n# 两位同学从本脚本产出的 .npz 接力，谁都不要重新读 DICOM，\n# 这样两条线起点完全相同，最终对比只差「提取器」本身。\n#\n# 流程：读DICOM → 转HU → 按物理位置排序 → lungmask肺分割\n#       → 重采样到各向同性 → 裁剪到肺包围盒 → 存 npz\n#\n# 在 Kaggle Notebook 上运行；需联网（lungmask 首次会下载权重）；\n# 有 GPU 会自动用 GPU 加速分割。\n# =====================================================\n\nimport os\nimport glob\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# ---------- 依赖安装（Kaggle 通常已有 pydicom；lungmask 需要装）----------\n# 若已安装可注释掉。lungmask 会一并带上 torch 依赖。\nos.system('pip -q install lungmask pydicom 2>/dev/null')\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom scipy.ndimage import zoom\nfrom lungmask import LMInferer\nimport torch\n\n# =====================================================\n# 配置区（两人务必用完全相同的配置，不要各改各的）\n# =====================================================\nBASE       = '/kaggle/input/competitions/osic-pulmonary-fibrosis-progression'\nOUT_DIR    = '/kaggle/working/ct_preprocessed'   # 中间文件输出目录\nSPLITS     = ['train', 'test']                   # train 和 test 都要处理\nTARGET_SPACING = (1.0, 1.0, 1.0)                 # 各向同性间距(mm)；磁盘紧可改 1.5 / 2.0\nHU_CLIP    = (-1024, 0)                           # HU 截断窗口（共用，锁死）\nMAX_PATIENTS = None                               # 先填 3 试跑，确认无误后改回 None 跑全量\n\nos.makedirs(OUT_DIR, exist_ok=True)\n\n# GPU 自动检测\nUSE_GPU = torch.cuda.is_available()\nprint(f\"分割使用设备: {'GPU' if USE_GPU else 'CPU'}\")\ninferer = LMInferer(force_cpu=not USE_GPU)  # 默认 U-net(R231) 模型\n\n\n# =====================================================\n# 1. 读取一个病人的所有切片，按物理位置排序\n# =====================================================\ndef load_scan(patient_dir):\n    files = glob.glob(os.path.join(patient_dir, '*.dcm'))\n    slices = []\n    for f in files:\n        try:\n            d = pydicom.dcmread(f)\n            # 必须能拿到像素和位置信息，否则跳过这张\n            _ = d.pixel_array\n            slices.append(d)\n        except Exception:\n            continue\n    if len(slices) == 0:\n        raise RuntimeError(\"无可用切片\")\n\n    # 优先用 ImagePositionPatient 的 z 坐标排序（最可靠），退化用 InstanceNumber\n    def sort_key(s):\n        try:\n            return float(s.ImagePositionPatient[2])\n        except Exception:\n            return float(getattr(s, 'InstanceNumber', 0))\n    slices.sort(key=sort_key)\n    return slices\n\n\n# =====================================================\n# 2. 转换为 HU，并堆叠成 3D 体积 (z, y, x)\n# =====================================================\ndef get_hu_volume(slices):\n    imgs = []\n    for s in slices:\n        arr = s.pixel_array.astype(np.float32)\n        slope = float(getattr(s, 'RescaleSlope', 1.0))\n        inter = float(getattr(s, 'RescaleIntercept', 0.0))\n        hu = arr * slope + inter\n        imgs.append(hu)\n    volume = np.stack(imgs).astype(np.int16)  # (z, y, x)\n    return volume\n\n\n# =====================================================\n# 3. 估计原始体素间距 (z, y, x)，含兜底\n# =====================================================\ndef get_spacing(slices):\n    # 平面内间距\n    try:\n        py, px = map(float, slices[0].PixelSpacing)\n    except Exception:\n        py, px = 1.0, 1.0\n\n    # 层间距：优先用相邻切片 z 坐标差，退化用 SliceThickness\n    pz = None\n    try:\n        z0 = float(slices[0].ImagePositionPatient[2])\n        z1 = float(slices[1].ImagePositionPatient[2])\n        pz = abs(z1 - z0)\n    except Exception:\n        pz = None\n    if not pz or pz == 0:\n        try:\n            pz = float(slices[0].SliceThickness)\n        except Exception:\n            pz = None\n    if not pz or pz == 0:\n        pz = 1.0  # 最终兜底，避免除零\n\n    return np.array([pz, py, px], dtype=np.float32)\n\n\n# =====================================================\n# 4. 重采样到各向同性间距\n#    volume 用线性插值(order=1)，mask 用最近邻(order=0)\n# =====================================================\ndef resample(arr, spacing, new_spacing, order):\n    factor = np.array(spacing) / np.array(new_spacing)\n    return zoom(arr, factor, order=order)\n\n\n# =====================================================\n# 5. 裁剪到肺掩膜的包围盒（去掉大片空气，省磁盘）\n# =====================================================\ndef crop_to_mask(volume, mask, margin=5):\n    coords = np.argwhere(mask > 0)\n    if len(coords) == 0:\n        return volume, mask  # 分割失败兜底，不裁剪\n    lo = np.maximum(coords.min(axis=0) - margin, 0)\n    hi = np.minimum(coords.max(axis=0) + margin + 1, np.array(volume.shape))\n    sl = tuple(slice(lo[i], hi[i]) for i in range(3))\n    return volume[sl], mask[sl]\n\n\n# =====================================================\n# 主流程：遍历所有病人，断点续跑\n# =====================================================\nmanifest = []   # 记录每个病人的处理结果，方便排查 & 写报告\n\nfor split in SPLITS:\n    split_dir = os.path.join(BASE, split)\n    patients = sorted(os.listdir(split_dir))\n    if MAX_PATIENTS:\n        patients = patients[:MAX_PATIENTS]\n\n    print(f\"\\n===== 处理 {split}：共 {len(patients)} 个病人 =====\")\n\n    for i, pid in enumerate(patients):\n        out_path = os.path.join(OUT_DIR, f'{split}__{pid}.npz')\n\n        # 断点续跑：已处理过的跳过\n        if os.path.exists(out_path):\n            print(f\"[{i+1}/{len(patients)}] {pid} 已存在，跳过\")\n            continue\n\n        try:\n            slices = load_scan(os.path.join(split_dir, pid))\n            n_slices = len(slices)\n            volume = get_hu_volume(slices)          # (z,y,x) HU\n            spacing = get_spacing(slices)           # 原始间距\n\n            # 在原始分辨率上做肺分割（lungmask 内部自行处理，质量最好）\n            mask = inferer.apply(volume)            # 返回 (z,y,x)，肺区为正\n            mask = (mask > 0).astype(np.uint8)\n\n            # 体积与掩膜一起重采样到目标各向同性间距\n            volume_rs = resample(volume, spacing, TARGET_SPACING, order=1).astype(np.int16)\n            mask_rs   = resample(mask,   spacing, TARGET_SPACING, order=0).astype(np.uint8)\n\n            # 形状对齐（插值后可能差 1 个像素）\n            z = min(volume_rs.shape[0], mask_rs.shape[0])\n            y = min(volume_rs.shape[1], mask_rs.shape[1])\n            x = min(volume_rs.shape[2], mask_rs.shape[2])\n            volume_rs = volume_rs[:z, :y, :x]\n            mask_rs   = mask_rs[:z, :y, :x]\n\n            # 统一 HU 截断\n            volume_rs = np.clip(volume_rs, HU_CLIP[0], HU_CLIP[1]).astype(np.int16)\n\n            # 裁剪到肺包围盒\n            volume_c, mask_c = crop_to_mask(volume_rs, mask_rs)\n\n            lung_voxels = int(mask_c.sum())\n\n            # 压缩保存：HU 体积 + 肺掩膜 + 元信息\n            np.savez_compressed(\n                out_path,\n                volume=volume_c,            # int16, HU, 已截断已裁剪\n                mask=mask_c,                # uint8, 肺=1\n                spacing=np.array(TARGET_SPACING, dtype=np.float32),\n            )\n\n            manifest.append({\n                'split': split, 'Patient': pid, 'status': 'ok',\n                'n_slices': n_slices,\n                'orig_spacing_z': float(spacing[0]),\n                'final_shape': str(volume_c.shape),\n                'lung_voxels': lung_voxels,\n            })\n            print(f\"[{i+1}/{len(patients)}] {pid}  切片{n_slices}→形状{volume_c.shape}  肺体素{lung_voxels}\")\n\n        except Exception as e:\n            manifest.append({\n                'split': split, 'Patient': pid, 'status': f'FAIL: {e}',\n                'n_slices': -1, 'orig_spacing_z': -1,\n                'final_shape': '', 'lung_voxels': -1,\n            })\n            print(f\"[{i+1}/{len(patients)}] {pid}  ✗ 失败: {e}\")\n\n# 保存处理清单\nmani_df = pd.DataFrame(manifest)\nmani_df.to_csv(os.path.join(OUT_DIR, 'manifest.csv'), index=False)\n\nprint(\"\\n===== 全部完成 =====\")\nprint(f\"成功: {(mani_df['status']=='ok').sum()}  失败: {(mani_df['status']!='ok').sum()}\")\nfails = mani_df[mani_df['status'] != 'ok']\nif len(fails):\n    print(\"失败病人（需人工关注）：\")\n    print(fails[['split', 'Patient', 'status']].to_string(index=False))\nprint(f\"\\n中间文件目录: {OUT_DIR}\")\nprint(\"两位同学从这里的 *.npz 接力：np.load(path) 取 ['volume'] 和 ['mask']\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-07T13:41:27.289905Z","iopub.execute_input":"2026-06-07T13:41:27.290335Z","iopub.status.idle":"2026-06-07T13:43:31.753040Z","shell.execute_reply.started":"2026-06-07T13:41:27.290306Z","shell.execute_reply":"2026-06-07T13:43:31.752056Z"}},"outputs":[],"execution_count":null}]}