{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":99552,"databundleVersionId":13441085}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport glob\nimport pickle\nfrom collections import defaultdict, namedtuple\n\nimport numpy as np\nimport pydicom\nfrom pydicom.dataset import Dataset\nfrom tqdm import tqdm\n\n# =========================\n# Configuration\n# =========================\nROOT_DIR = \"/kaggle/input/rsna-aneurysm/train\"  # set to dataset root\nOUT_DIR = \"./preproc_rsna_aneurysm\"             # outputs: pickles + npy volumes\nos.makedirs(OUT_DIR, exist_ok=True)\n\n# Options\nLOAD_PIXEL_DATA = True           # set False if only indexing/metadata needed\nAPPLY_HU_CONVERSION = True       # convert to HU using RescaleSlope/Intercept\nAPPLY_WINDOWING = False          # typical CTA window, toggle as needed\nWINDOW_CENTER = 40.0\nWINDOW_WIDTH = 400.0\nCLIP_TO_WINDOW = True            # clip to [WL-0.5, WH-0.5] after windowing\nFLOAT32_OUTPUT = True            # store arrays as float32\nSAVE_PER_SERIES_VOL = True       # save each 3D volume as .npy\nSAVE_NEIGHBORS = True            # map neighboring slice SOPInstanceUIDs\n\n# =========================\n# Helpers\n# =========================\nSeriesKey = namedtuple(\"SeriesKey\", [\"study_uid\", \"series_uid\"])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:55:42.103146Z","iopub.execute_input":"2025-08-23T17:55:42.103521Z","iopub.status.idle":"2025-08-23T17:55:42.87792Z","shell.execute_reply.started":"2025-08-23T17:55:42.103488Z","shell.execute_reply":"2025-08-23T17:55:42.876939Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def window_image(img, center, width, clip=True):\n    # Linear windowing\n    low = center - width / 2.0\n    high = center + width / 2.0\n    img = (img - low) / (high - low)\n    img = img * 255.0\n    if clip:\n        img = np.clip(img, 0, 255)\n    return img","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:55:52.482257Z","iopub.execute_input":"2025-08-23T17:55:52.482578Z","iopub.status.idle":"2025-08-23T17:55:52.488026Z","shell.execute_reply.started":"2025-08-23T17:55:52.482545Z","shell.execute_reply":"2025-08-23T17:55:52.486963Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def safe_float(x, default=np.nan):\n    try:\n        return float(x)\n    except Exception:\n        return default\n\ndef is_multiframe(ds: Dataset):\n    # Enhanced CT can be multi-frame with PerFrameFunctionalGroupsSequence\n    return hasattr(ds, \"NumberOfFrames\") and ds.NumberOfFrames is not None and ds.NumberOfFrames > 1","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:55:56.08439Z","iopub.execute_input":"2025-08-23T17:55:56.084798Z","iopub.status.idle":"2025-08-23T17:55:56.090178Z","shell.execute_reply.started":"2025-08-23T17:55:56.084768Z","shell.execute_reply":"2025-08-23T17:55:56.089363Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def extract_orientation_normal(iop):\n    # iop length 6: row (0..2), col (3..5)\n    r = np.array(iop[0:3], dtype=float)\n    c = np.array(iop[3:6], dtype=float)\n    n = np.cross(r, c)\n    return n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:00.393919Z","iopub.execute_input":"2025-08-23T17:56:00.394692Z","iopub.status.idle":"2025-08-23T17:56:00.399034Z","shell.execute_reply.started":"2025-08-23T17:56:00.394662Z","shell.execute_reply":"2025-08-23T17:56:00.398079Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def sort_by_geometry(file_list):\n    # Robust geometric sorting using IOP+IPP; fallback to InstanceNumber; fallback to as-is\n    headers = []\n    for f in file_list:\n        ds = pydicom.dcmread(f, stop_before_pixels=True, force=True)\n        ipp = getattr(ds, \"ImagePositionPatient\", None)\n        iop = getattr(ds, \"ImageOrientationPatient\", None)\n        inst = getattr(ds, \"InstanceNumber\", None)\n        headers.append((f, ipp, iop, inst))\n\n    # Try geometric sort\n    if all(h[1] is not None and h[2] is not None for h in headers):\n        n = extract_orientation_normal(headers[2])\n        positions = []\n        for _, ipp, _, _ in headers:\n            positions.append(np.dot(np.array(ipp, dtype=float), n))\n        order = np.argsort(positions)\n        return [headers[i] for i in order]\n\n    # Fallback: InstanceNumber\n    if all(h[3] is not None for h in headers):\n        order = np.argsort([h[3] for h in headers])\n        return [headers[i] for i in order]\n\n    # Last resort: filesystem order\n    return [h for h in headers]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:04.009322Z","iopub.execute_input":"2025-08-23T17:56:04.009663Z","iopub.status.idle":"2025-08-23T17:56:04.017769Z","shell.execute_reply.started":"2025-08-23T17:56:04.009636Z","shell.execute_reply":"2025-08-23T17:56:04.016865Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def index_dicoms(root_dir):\n    # Recursively find all .dcm and group by Study/Series\n    series_map = defaultdict(list)\n    for f in glob.glob(os.path.join(root_dir, \"**\", \"*.dcm\"), recursive=True):\n        try:\n            ds = pydicom.dcmread(f, stop_before_pixels=True, force=True)\n            study_uid = getattr(ds, \"StudyInstanceUID\", None)\n            series_uid = getattr(ds, \"SeriesInstanceUID\", None)\n            if study_uid and series_uid:\n                series_map[SeriesKey(study_uid, series_uid)].append(f)\n        except Exception:\n            # Ignore unreadable\n            continue\n    return series_map","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:08.894368Z","iopub.execute_input":"2025-08-23T17:56:08.894671Z","iopub.status.idle":"2025-08-23T17:56:08.900597Z","shell.execute_reply.started":"2025-08-23T17:56:08.894647Z","shell.execute_reply":"2025-08-23T17:56:08.899705Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def read_multiframe_volume(ds: Dataset):\n    # Extract 3D array and per-frame geometry from multi-frame CT\n    n_frames = int(ds.NumberOfFrames)\n    rows = int(ds.Rows)\n    cols = int(ds.Columns)\n\n    # Pixel data read\n    vol = ds.pixel_array  # shape: (frames, rows, cols) or (rows, cols, frames) depending on handler\n    if vol.shape == rows and vol.shape[1] == cols and vol.shape[-1] == n_frames:\n        # Make it (Z, Y, X)\n        vol = np.moveaxis(vol, -1, 0)\n    elif vol.shape == n_frames:\n        # already (Z, Y, X)\n        pass\n    else:\n        # Attempt to coerce\n        vol = vol.reshape((n_frames, rows, cols))\n\n    # Spacing\n    pixsp = getattr(ds, \"PixelSpacing\", None)\n    dz = safe_float(getattr(ds, \"SpacingBetweenSlices\", getattr(ds, \"SliceThickness\", np.nan)))\n    dy = safe_float(pixsp) if pixsp is not None else np.nan\n    dx = safe_float(pixsp[1]) if pixsp is not None else np.nan\n\n    slope = safe_float(getattr(ds, \"RescaleSlope\", 1.0), 1.0)\n    intercept = safe_float(getattr(ds, \"RescaleIntercept\", 0.0), 0.0)\n\n    meta = dict(\n        modality=getattr(ds, \"Modality\", None),\n        manufacturer=getattr(ds, \"Manufacturer\", None),\n        convolution_kernel=str(getattr(ds, \"ConvolutionKernel\", \"\")),\n        kVp=safe_float(getattr(ds, \"KVP\", np.nan)),\n        exposure=safe_float(getattr(ds, \"Exposure\", np.nan)),\n        pixel_spacing=(dy, dx, dz),\n        rescale_slope=slope,\n        rescale_intercept=intercept,\n        series_description=str(getattr(ds, \"SeriesDescription\", \"\")),\n        study_uid=getattr(ds, \"StudyInstanceUID\", None),\n        series_uid=getattr(ds, \"SeriesInstanceUID\", None)\n    )\n    return vol, meta","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:12.90536Z","iopub.execute_input":"2025-08-23T17:56:12.905677Z","iopub.status.idle":"2025-08-23T17:56:12.915314Z","shell.execute_reply.started":"2025-08-23T17:56:12.905652Z","shell.execute_reply":"2025-08-23T17:56:12.914067Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def read_singleframe_volume(sorted_files):\n    # Load slices into a 3D volume (Z, Y, X)\n    slices = []\n    metas = []\n    for f in sorted_files:\n        ds = pydicom.dcmread(f, force=True)\n        arr = ds.pixel_array\n        slices.append(arr)\n        pixsp = getattr(ds, \"PixelSpacing\", None)\n        dy = safe_float(pixsp) if pixsp is not None else np.nan\n        dx = safe_float(pixsp[1]) if pixsp is not None else np.nan\n        dz = safe_float(getattr(ds, \"SpacingBetweenSlices\", getattr(ds, \"SliceThickness\", np.nan)))\n        slope = safe_float(getattr(ds, \"RescaleSlope\", 1.0), 1.0)\n        intercept = safe_float(getattr(ds, \"RescaleIntercept\", 0.0), 0.0)\n        meta = dict(\n            sop_uid=str(getattr(ds, \"SOPInstanceUID\", \"\")),\n            instance_number=getattr(ds, \"InstanceNumber\", None),\n            image_position=getattr(ds, \"ImagePositionPatient\", None),\n            image_orientation=getattr(ds, \"ImageOrientationPatient\", None),\n            pixel_spacing=(dy, dx),\n            slice_thickness=safe_float(getattr(ds, \"SliceThickness\", np.nan)),\n            spacing_between_slices=safe_float(getattr(ds, \"SpacingBetweenSlices\", np.nan)),\n            rescale_slope=slope,\n            rescale_intercept=intercept,\n            convolution_kernel=str(getattr(ds, \"ConvolutionKernel\", \"\")),\n            kVp=safe_float(getattr(ds, \"KVP\", np.nan)),\n            exposure=safe_float(getattr(ds, \"Exposure\", np.nan)),\n        )\n        metas.append(meta)\n\n    vol = np.stack(slices, axis=0)  # (Z, Y, X)\n\n    # Attempt consistent spacing\n    dy = metas[\"pixel_spacing\"] if metas and metas[\"pixel_spacing\"] == metas[\"pixel_spacing\"] else np.nan\n    dx = metas[\"pixel_spacing\"][1] if metas and metas[\"pixel_spacing\"][1] == metas[\"pixel_spacing\"][1] else np.nan\n    # Prefer SpacingBetweenSlices; fallback to SliceThickness; compute from IPP if needed\n    dz = metas[\"spacing_between_slices\"]\n    if np.isnan(dz) or dz == 0:\n        dz = metas[\"slice_thickness\"]\n    # IPP-based spacing if available\n    z_positions = []\n    has_ipp = all(m[\"image_position\"] is not None for m in metas)\n    has_iop = all(m[\"image_orientation\"] is not None for m in metas)\n    if has_ipp and has_iop:\n        n = extract_orientation_normal(metas[\"image_orientation\"])\n        for m in metas:\n            pos = np.dot(np.array(m[\"image_position\"], dtype=float), n)\n            z_positions.append(pos)\n        z_positions = np.array(z_positions)\n        diffs = np.diff(np.sort(z_positions))\n        if diffs.size > 0:\n            dz_est = float(np.median(np.abs(diffs)))\n            if not np.isnan(dz_est) and dz_est > 0:\n                dz = dz_est\n\n    slope = metas[\"rescale_slope\"] if metas else 1.0\n    intercept = metas[\"rescale_intercept\"] if metas else 0.0\n\n    series_meta = dict(\n        pixel_spacing=(dy, dx, dz),\n        rescale_slope=slope,\n        rescale_intercept=intercept,\n        convolution_kernel=metas[\"convolution_kernel\"] if metas else \"\",\n        kVp=metas[\"kVp\"] if metas else np.nan,\n        exposure=metas[\"exposure\"] if metas else np.nan,\n        per_slice=metas\n    )\n    return vol, series_meta","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:18.029316Z","iopub.execute_input":"2025-08-23T17:56:18.029642Z","iopub.status.idle":"2025-08-23T17:56:18.042815Z","shell.execute_reply.started":"2025-08-23T17:56:18.029617Z","shell.execute_reply":"2025-08-23T17:56:18.041821Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def to_hounsfield(vol, slope, intercept):\n    # Convert to HU in float32\n    vol = vol.astype(np.float32)\n    vol = vol * (slope if slope is not None and not np.isnan(slope) else 1.0)\n    vol = vol + (intercept if intercept is not None and not np.isnan(intercept) else 0.0)\n    return vol\n\ndef build_neighbors(sorted_files):\n    # Return dict sop_uid -> (prev_uid, next_uid, filepath)\n    uids = []\n    for f in sorted_files:\n        ds = pydicom.dcmread(f, stop_before_pixels=True, force=True)\n        uids.append((str(getattr(ds, \"SOPInstanceUID\", \"\")), f))\n    n = len(uids)\n    nb = {}\n    for i, (uid, f) in enumerate(uids):\n        prev_uid = uids[i-1] if i > 0 else uid\n        next_uid = uids[i+1] if i < n-1 else uid\n        nb[uid] = dict(prev_uid=prev_uid, next_uid=next_uid, path=f)\n    return nb","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T17:56:25.719116Z","iopub.execute_input":"2025-08-23T17:56:25.719662Z","iopub.status.idle":"2025-08-23T17:56:25.726473Z","shell.execute_reply.started":"2025-08-23T17:56:25.719626Z","shell.execute_reply":"2025-08-23T17:56:25.725775Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =========================\n# Main pipeline\n# =========================\ndef main():\n    series_map = index_dicoms(ROOT_DIR)\n    print(f\"Found {len(series_map)} series\")\n\n    series_index = {}     # key: (study_uid, series_uid) -> metadata summary\n    image_index = {}      # key: sop_uid -> metadata (neighbors, z proxy, exposure, thickness, path)\n    saved_volumes = []    # paths to saved npy volumes\n\n    for key, files in tqdm(series_map.items(), total=len(series_map)):\n        # Detect multi-frame vs single-frame\n        first_ds = pydicom.dcmread(files, stop_before_pixels=True, force=True)\n        study_uid = getattr(first_ds, \"StudyInstanceUID\", None)\n        series_uid = getattr(first_ds, \"SeriesInstanceUID\", None)\n        series_desc = str(getattr(first_ds, \"SeriesDescription\", \"\"))\n\n        if is_multiframe(first_ds):\n            if not LOAD_PIXEL_DATA:\n                # Index-only\n                series_index[key] = dict(\n                    study_uid=study_uid,\n                    series_uid=series_uid,\n                    series_description=series_desc,\n                    n_slices=int(first_ds.NumberOfFrames),\n                    pixel_spacing=None,\n                )\n            else:\n                ds_full = pydicom.dcmread(files, force=True)  # need PixelData\n                vol, meta = read_multiframe_volume(ds_full)\n                if APPLY_HU_CONVERSION:\n                    vol = to_hounsfield(vol, meta[\"rescale_slope\"], meta[\"rescale_intercept\"])\n                if APPLY_WINDOWING:\n                    vol = window_image(vol, WINDOW_CENTER, WINDOW_WIDTH, clip=CLIP_TO_WINDOW)\n                if FLOAT32_OUTPUT and vol.dtype != np.float32:\n                    vol = vol.astype(np.float32)\n\n                save_path = os.path.join(OUT_DIR, f\"{study_uid}__{series_uid}.npy\")\n                if SAVE_PER_SERIES_VOL:\n                    np.save(save_path, vol)\n                    saved_volumes.append(save_path)\n\n                # Build a minimal image_index using frame numbers\n                n = vol.shape\n                for i in range(n):\n                    sop_uid = f\"{series_uid}_frame_{i+1}\"\n                    prev_uid = f\"{series_uid}_frame_{max(1, i)}\"\n                    next_uid = f\"{series_uid}_frame_{min(n, i+2)}\"\n                    image_index[sop_uid] = dict(\n                        study_id=study_uid,\n                        series_id=series_uid,\n                        image_minus1=prev_uid,\n                        image_plus1=next_uid,\n                        path=f\"{files}#frame={i+1}\",\n                        z_pos=np.nan,\n                        exposure=meta.get(\"exposure\", np.nan),\n                        thickness=meta.get(\"pixel_spacing\", (np.nan, np.nan, np.nan))[2]\n                    )\n\n                series_index[key] = dict(\n                    study_uid=study_uid,\n                    series_uid=series_uid,\n                    series_description=series_desc,\n                    n_slices=n,\n                    pixel_spacing=meta[\"pixel_spacing\"],\n                    saved_path=save_path if SAVE_PER_SERIES_VOL else None\n                )\n\n        else:\n            # Multi-file single-frame series\n            sorted_files = sort_by_geometry(files)\n            if SAVE_NEIGHBORS:\n                nb = build_neighbors(sorted_files)\n                image_index.update({\n                    k: dict(\n                        study_id=study_uid,\n                        series_id=series_uid,\n                        image_minus1=v[\"prev_uid\"],\n                        image_plus1=v[\"next_uid\"],\n                        path=v[\"path\"],\n                        z_pos=np.nan,  # filled later if IPP available\n                        exposure=np.nan,\n                        thickness=np.nan\n                    ) for k, v in nb.items()\n                })\n\n            if not LOAD_PIXEL_DATA:\n                series_index[key] = dict(\n                    study_uid=study_uid,\n                    series_uid=series_uid,\n                    series_description=series_desc,\n                    n_slices=len(sorted_files),\n                    pixel_spacing=None,\n                )\n                continue\n\n            vol, meta = read_singleframe_volume(sorted_files)\n            if APPLY_HU_CONVERSION:\n                vol = to_hounsfield(vol, meta[\"rescale_slope\"], meta[\"rescale_intercept\"])\n            if APPLY_WINDOWING:\n                vol = window_image(vol, WINDOW_CENTER, WINDOW_WIDTH, clip=CLIP_TO_WINDOW)\n            if FLOAT32_OUTPUT and vol.dtype != np.float32:\n                vol = vol.astype(np.float32)\n\n            # Update z_pos, exposure, thickness in image_index if neighbors computed\n            if SAVE_NEIGHBORS and \"per_slice\" in meta:\n                # derive z from IPP along normal if available\n                zvals = []\n                has_ipp = all(m.get(\"image_position\") is not None for m in meta[\"per_slice\"])\n                has_iop = all(m.get(\"image_orientation\") is not None for m in meta[\"per_slice\"])\n                nvec = None\n                if has_ipp and has_iop:\n                    nvec = extract_orientation_normal(meta[\"per_slice\"][0][\"image_orientation\"])\n                    for m in meta[\"per_slice\"]:\n                        zvals.append(float(np.dot(np.array(m[\"image_position\"], dtype=float), nvec)))\n                else:\n                    zvals = [np.nan] * len(meta[\"per_slice\"])\n\n                for i, m in enumerate(meta[\"per_slice\"]):\n                    sop_uid = m[\"sop_uid\"]\n                    if sop_uid in image_index:\n                        image_index[sop_uid][\"z_pos\"] = zvals[i]\n                        image_index[sop_uid][\"exposure\"] = m[\"exposure\"]\n                        # Prefer spacing_between_slices; fallback to slice_thickness\n                        thick = m[\"spacing_between_slices\"]\n                        if np.isnan(thick) or thick == 0:\n                            thick = m[\"slice_thickness\"]\n                        image_index[sop_uid][\"thickness\"] = thick\n\n            save_path = os.path.join(OUT_DIR, f\"{study_uid}__{series_uid}.npy\")\n            if SAVE_PER_SERIES_VOL:\n                np.save(save_path, vol)\n                saved_volumes.append(save_path)\n\n            dz = meta[\"pixel_spacing\"][2] if meta[\"pixel_spacing\"] else np.nan\n            series_index[key] = dict(\n                study_uid=study_uid,\n                series_uid=series_uid,\n                series_description=series_desc,\n                n_slices=vol.shape,\n                pixel_spacing=meta[\"pixel_spacing\"],\n                saved_path=save_path if SAVE_PER_SERIES_VOL else None\n            )\n\n    # Persist indices\n    with open(os.path.join(OUT_DIR, \"series_index.pickle\"), \"wb\") as f:\n        pickle.dump(series_index, f, protocol=pickle.HIGHEST_PROTOCOL)\n    with open(os.path.join(OUT_DIR, \"image_index.pickle\"), \"wb\") as f:\n        pickle.dump(image_index, f, protocol=pickle.HIGHEST_PROTOCOL)\n    with open(os.path.join(OUT_DIR, \"saved_volumes.pickle\"), \"wb\") as f:\n        pickle.dump(saved_volumes, f, protocol=pickle.HIGHEST_PROTOCOL)\n\n    print(f\"Series indexed: {len(series_index)}\")\n    print(f\"Images indexed: {len(image_index)}\")\n    print(f\"Volumes saved: {len(saved_volumes)}\")\n\nif __name__ == \"__main__\":\n    main()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}