{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":99552,"databundleVersionId":13762876,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"%%capture output\n\n!pip install dicom","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:15:38.657093Z","iopub.execute_input":"2025-09-21T11:15:38.657444Z","iopub.status.idle":"2025-09-21T11:15:42.454763Z","shell.execute_reply.started":"2025-09-21T11:15:38.657419Z","shell.execute_reply":"2025-09-21T11:15:42.453612Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!apt install ffmpeg","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:16:05.410491Z","iopub.execute_input":"2025-09-21T11:16:05.411Z","iopub.status.idle":"2025-09-21T11:16:08.950738Z","shell.execute_reply.started":"2025-09-21T11:16:05.410959Z","shell.execute_reply":"2025-09-21T11:16:08.949579Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\nimport os\nimport random\n\nfrom PIL import Image\n\nfrom types import SimpleNamespace\n\nimport polars as pl\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\n\nplt.style.use('ggplot')\n\nimport pydicom as dicom\nimport pydicom","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:20:16.990296Z","iopub.execute_input":"2025-09-21T11:20:16.990585Z","iopub.status.idle":"2025-09-21T11:20:16.997146Z","shell.execute_reply.started":"2025-09-21T11:20:16.990564Z","shell.execute_reply":"2025-09-21T11:20:16.995963Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cfg = SimpleNamespace()\ncfg.INPUT = Path(\"/kaggle/input/rsna-intracranial-aneurysm-detection\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.239076Z","iopub.execute_input":"2025-09-21T11:04:56.239431Z","iopub.status.idle":"2025-09-21T11:04:56.257633Z","shell.execute_reply.started":"2025-09-21T11:04:56.239399Z","shell.execute_reply":"2025-09-21T11:04:56.256555Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = pl.read_csv(cfg.INPUT / \"train.csv\")\ntrain_localizers = pl.read_csv(cfg.INPUT / \"train_localizers.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.258916Z","iopub.execute_input":"2025-09-21T11:04:56.259221Z","iopub.status.idle":"2025-09-21T11:04:56.284646Z","shell.execute_reply.started":"2025-09-21T11:04:56.259192Z","shell.execute_reply":"2025-09-21T11:04:56.283525Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"# of training data: {len(train)}\")\nprint(f\"# of train localizers: {len(train_localizers['SeriesInstanceUID'].unique())}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.287359Z","iopub.execute_input":"2025-09-21T11:04:56.287664Z","iopub.status.idle":"2025-09-21T11:04:56.295988Z","shell.execute_reply.started":"2025-09-21T11:04:56.287641Z","shell.execute_reply":"2025-09-21T11:04:56.294897Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Age distribution in Aneurysm patients","metadata":{}},{"cell_type":"code","source":"absent, present = train.partition_by(\"Aneurysm Present\")\n\nfig, axes = plt.subplots(1,2, figsize=(8,4))\naxes[0].hist(\n    absent['PatientAge'],\n    bins=30,\n)\naxes[0].set_title(\"Absent\")\naxes[1].hist(\n    present['PatientAge'],\n    bins=30,\n)\naxes[1].set_title(\"Present\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.296931Z","iopub.execute_input":"2025-09-21T11:04:56.297249Z","iopub.status.idle":"2025-09-21T11:04:56.720691Z","shell.execute_reply.started":"2025-09-21T11:04:56.297219Z","shell.execute_reply":"2025-09-21T11:04:56.71972Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"present.group_by(\"PatientSex\").agg(count=pl.col('Aneurysm Present').count())\n\nplt.bar(\n    present.group_by(\"PatientSex\").agg(count=pl.col('Aneurysm Present').count())['PatientSex'],\n    present.group_by(\"PatientSex\").agg(count=pl.col('Aneurysm Present').count())['count']\n)\nplt.title(\"Aneurysm Positive\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.72176Z","iopub.execute_input":"2025-09-21T11:04:56.722027Z","iopub.status.idle":"2025-09-21T11:04:56.866005Z","shell.execute_reply.started":"2025-09-21T11:04:56.722006Z","shell.execute_reply":"2025-09-21T11:04:56.864897Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Show a random image","metadata":{}},{"cell_type":"code","source":"series_sop_df = train_localizers.group_by('SeriesInstanceUID').agg(\n    SOPInstanceUIDs = pl.col('SOPInstanceUID'),\n    SOPInstanceUIDCount = pl.col('SOPInstanceUID').count(),\n    coordinates = pl.col('coordinates')\n).sort('SOPInstanceUIDCount')\n\nsop = series_sop_df.row(by_predicate=(\n        pl.col('SeriesInstanceUID') == random.choice(series_sop_df['SeriesInstanceUID'])\n))\n\nfig, axes = plt.subplots(1, len(sop[1]), figsize=(8,8))\n\nif len(sop[1]) > 1:\n    for i, ax in enumerate(axes.ravel()):\n        image = dicom.dcmread(cfg.INPUT / \"series\" / f\"{sop[0]}\" / f\"{sop[1][i]}.dcm\").pixel_array\n        \n        # there are some 3d data\n        if(len(image.shape) > 2):\n            image = image[1, :, :]\n            \n        ax.imshow(\n            image,\n            cmap=plt.cm.bone\n        )\n        ax.plot(\n            eval(sop[3][i])['x'], eval(sop[3][i])['y'],\n            marker='x',\n            color='red'\n        )\nelif len(sop[1]) == 1:\n    image = dicom.dcmread(cfg.INPUT / \"series\" / f\"{sop[0]}\" / f\"{sop[1][0]}.dcm\").pixel_array\n    \n    # there are some 3d data\n    if(len(image.shape) > 2):\n        image = image[1, :, :]\n        \n    axes.imshow(\n        image,\n        cmap=plt.cm.bone\n    )\n    axes.plot(\n        eval(sop[3][0])['x'], eval(sop[3][0])['y'],\n        marker='x',\n        color='red'\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:04:56.86739Z","iopub.execute_input":"2025-09-21T11:04:56.867718Z","iopub.status.idle":"2025-09-21T11:04:57.247994Z","shell.execute_reply.started":"2025-09-21T11:04:56.867694Z","shell.execute_reply":"2025-09-21T11:04:57.246976Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Create animation","metadata":{}},{"cell_type":"code","source":"train_localizers = train_localizers.with_columns(\n    pl.format(\n        \"/kaggle/input/rsna-intracranial-aneurysm-detection/series/{}/{}.dcm\",\n        pl.col(\"SeriesInstanceUID\"),\n        pl.col(\"SOPInstanceUID\")\n    ).alias(\"filename\")\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:06:21.61573Z","iopub.execute_input":"2025-09-21T11:06:21.61608Z","iopub.status.idle":"2025-09-21T11:06:21.625028Z","shell.execute_reply.started":"2025-09-21T11:06:21.616034Z","shell.execute_reply":"2025-09-21T11:06:21.624021Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# params\npaths = [p for p in train_localizers.get_column(\"filename\").to_list() if p and os.path.exists(p)]\nwant = 5\nout = \"five_multiframe_sidebyside.gif\"\nfps = 8\n\ndef load_multiframe(path):\n    ds = pydicom.dcmread(path, force=True)\n    arr = ds.pixel_array.astype(np.float32)\n    # make sure frames dim exists\n    if arr.ndim == 2:\n        arr = arr[np.newaxis, ...]\n    elif arr.ndim == 3:\n        # assume (F,H,W) if first dim >1 (multi-frame)\n        if arr.shape[0] <= 4 and arr.shape[0] != arr.shape[1]:\n            # sometimes (H,W,C) -> transpose to (H,W,C) handled later; treat as single frame\n            arr = arr[np.newaxis, ...]\n    # rescale\n    slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n    if slope != 1.0 or intercept != 0.0:\n        arr = arr * slope + intercept\n    # photometric inversion for MONOCHROME1\n    if getattr(ds, \"PhotometricInterpretation\", \"\").upper() == \"MONOCHROME1\":\n        arr = np.max(arr) - arr\n    return arr  # shape (F, H, W) ideally\n\ndef to_uint8(img2d, p1=1, p99=99):\n    lo, hi = np.percentile(img2d, (p1, p99))\n    if hi <= lo:\n        lo, hi = img2d.min(), img2d.max()\n    if hi == lo:\n        return np.zeros_like(img2d, dtype=np.uint8)\n    clipped = np.clip(img2d, lo, hi)\n    return ((clipped - lo) / (hi - lo) * 255.0).astype(np.uint8)\n\ndef ensure_size(img, shape):\n    # img: 2D uint8, shape: (H, W)\n    if img.shape == shape:\n        return img\n    img_pil = Image.fromarray(img).resize((shape[1], shape[0]), Image.LANCZOS)\n    return np.asarray(img_pil)\n\ndef make_rgb(img2d):\n    return np.stack([img2d]*3, axis=-1)  # (H,W,3)\n\n# find first 5 multi-frame files\nmulti = []\nfor p in paths:\n    try:\n        arr = pydicom.dcmread(p, force=True).pixel_array\n        if getattr(arr, \"ndim\", 2) == 3 and arr.shape[0] > 1:\n            multi.append(p)\n    except Exception:\n        continue\n    if len(multi) >= want:\n        break\n\nif len(multi) < want:\n    raise RuntimeError(f\"Found only {len(multi)} multi-frame DICOM(s); need {want}.\")\n\n# load arrays for chosen files\narrs = [load_multiframe(p) for p in multi[:want]]\nnframes = min(a.shape[0] for a in arrs)  # sync length by minimum\n\n# choose reference shape from first file's single frame (H, W)\nref_shape = arrs[0].shape[1], arrs[0].shape[2]\n\n# build first combined frame for fig setup\nfirst_pieces = []\nfor a in arrs:\n    img = to_uint8(a[0])\n    img = ensure_size(img, ref_shape)\n    first_pieces.append(make_rgb(img))\nfirst_frame = np.concatenate(first_pieces, axis=1)\n\nfig, ax = plt.subplots(figsize=(first_frame.shape[1]/100, first_frame.shape[0]/100), dpi=100)\nax.axis(\"off\")\nim = ax.imshow(first_frame, animated=True)\n\ndef gen():\n    for i in range(nframes):\n        pieces = []\n        for a in arrs:\n            img = a[i]\n            # if (C,H,W) channel-first unlikely here; handle only (H,W)\n            if img.ndim == 3 and img.shape[0] in (3,4):\n                img = np.transpose(img, (1,2,0))  # to H,W,C\n                # convert color -> grayscale by mean\n                img = img.mean(axis=2)\n            img8 = to_uint8(img)\n            img8 = ensure_size(img8, ref_shape)\n            pieces.append(make_rgb(img8))\n        yield np.concatenate(pieces, axis=1)\n\ndef update(frame):\n    im.set_array(frame)\n    return (im,)\n\nani = animation.FuncAnimation(fig, update, frames=gen(), interval=1000/fps, blit=True)\nani.save(out, writer=\"pillow\", fps=fps)\nplt.close(fig)\nprint(\"Saved:\", out)\nfrom IPython.display import HTML\nHTML(ani.to_jshtml())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T11:22:59.484503Z","iopub.execute_input":"2025-09-21T11:22:59.484963Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, math, ast\nimport numpy as np\nimport pydicom\nfrom PIL import Image\nfrom tqdm import tqdm\n\nfilenames = [f for f in train_localizers.get_column(\"filename\").to_list() if f and os.path.exists(f)]\ncoords_raw = train_localizers.get_column(\"coordinates\").to_list()\nout_gif = \"train_localizers_with_cross.gif\"\nfps = 8\ncross_half = 6\ncross_thickness = 2\nsample_frames_for_window = 8\n\ndef parse_coord(raw):\n    if raw is None:\n        return None\n    if isinstance(raw, dict):\n        return raw\n    try:\n        return ast.literal_eval(raw)\n    except Exception:\n        try:\n            import json\n            return json.loads(raw)\n        except Exception:\n            return None\n\ndef read_pixel_array(path):\n    ds = pydicom.dcmread(path, force=True)\n    arr = ds.pixel_array.astype(np.float32)\n    # If multi-frame (F,H,W), pick middle frame (you can change to mean/max if desired)\n    if arr.ndim == 3:\n        arr = arr[arr.shape[0] // 2]\n    # rescale\n    slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n    if slope != 1.0 or intercept != 0.0:\n        arr = arr * slope + intercept\n    # photometric inversion for MONOCHROME1\n    if getattr(ds, \"PhotometricInterpretation\", \"\").upper() == \"MONOCHROME1\":\n        arr = np.max(arr) - arr\n    return arr\n\nn = len(filenames)\nindices = np.linspace(0, n-1, min(sample_frames_for_window, n), dtype=int)\nsamples = []\nfor i in indices:\n    try:\n        samples.append(read_pixel_array(filenames[i]).ravel())\n    except Exception:\n        continue\nif not samples:\n    raise RuntimeError(\"Failed to read sample frames for window estimation.\")\nall_sample_pixels = np.concatenate(samples)\np1, p99 = np.percentile(all_sample_pixels, (1, 99))\nif p99 <= p1:\n    p1, p99 = all_sample_pixels.min(), all_sample_pixels.max()\nif p99 == p1:\n    p99 = p1 + 1.0\n\ndef to_uint8_with_window(arr, lo=p1, hi=p99):\n    arr = np.clip(arr, lo, hi)\n    scaled = ((arr - lo) / (hi - lo) * 255.0).astype(np.uint8)\n    return scaled\n\nimages = []\nfor i, path in enumerate(tqdm(filenames, desc=\"Preparing frames\")):\n    try:\n        arr = read_pixel_array(path)\n    except Exception:\n        # fallback blank image\n        arr = np.zeros_like(read_pixel_array(filenames[0]))\n    arr8 = to_uint8_with_window(arr)            # 2D uint8\n    # convert to RGB (H, W, 3)\n    rgb = np.stack([arr8, arr8, arr8], axis=-1)\n    # parse coord and draw cross by setting pixel blocks (fast)\n    coord = parse_coord(coords_raw[i]) if i < len(coords_raw) else None\n    if coord and \"x\" in coord and \"y\" in coord:\n        try:\n            x = int(round(float(coord[\"x\"])))\n            y = int(round(float(coord[\"y\"])))\n            H, W = arr8.shape\n            # clamp center\n            x = max(0, min(W - 1, x))\n            y = max(0, min(H - 1, y))\n            # horizontal bar\n            x0 = max(0, x - cross_half)\n            x1 = min(W, x + cross_half + 1)\n            y0 = max(0, y - cross_thickness // 2)\n            y1 = min(H, y + (cross_thickness + 1)//2)\n            rgb[y0:y1, x0:x1, 0] = 255  # R channel\n            rgb[y0:y1, x0:x1, 1] = 0\n            rgb[y0:y1, x0:x1, 2] = 0\n            # vertical bar\n            x0_v = max(0, x - cross_thickness // 2)\n            x1_v = min(W, x + (cross_thickness + 1)//2)\n            y0_v = max(0, y - cross_half)\n            y1_v = min(H, y + cross_half + 1)\n            rgb[y0_v:y1_v, x0_v:x1_v, 0] = 255\n            rgb[y0_v:y1_v, x0_v:x1_v, 1] = 0\n            rgb[y0_v:y1_v, x0_v:x1_v, 2] = 0\n        except Exception:\n            pass\n    # convert to PIL Image and append\n    images.append(Image.fromarray(rgb))\n\nduration_ms = int(1000 / fps)\nimages[0].save(out_gif, save_all=True, append_images=images[1:], duration=duration_ms, loop=0)\nprint(\"Saved:\", out_gif)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%capture output\n\n# requires ffmpeg installed and on PATH\n!ffmpeg -y -i train_localizers_with_cross.gif -movflags +faststart -pix_fmt yuv420p -vcodec libx264 -crf 23 train_localizers_with_cross.mp4","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from IPython.display import Video, display\n\n# local file\ndisplay(Video(\"train_localizers_with_cross.mp4\", embed=True, width=640, height=360))","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}