{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Creating images from axial CT scans of the abdomen\nIn this notebook, we transform medical images stored in the DICOM format into PNG format.","metadata":{}},{"cell_type":"markdown","source":"<a id = \"table-of-content\"></a>\n## Table of Contents\n- [Setup](#setup)\n- [Transforming DICOM files using Multiprocessing](#core)\n- [Samples](#result)","metadata":{"execution":{"iopub.status.busy":"2023-10-13T18:50:03.253442Z","iopub.execute_input":"2023-10-13T18:50:03.254254Z","iopub.status.idle":"2023-10-13T18:50:03.267168Z","shell.execute_reply.started":"2023-10-13T18:50:03.254216Z","shell.execute_reply":"2023-10-13T18:50:03.264629Z"}}},{"cell_type":"markdown","source":"<a id = \"setup\"></a>\n# Setup","metadata":{}},{"cell_type":"markdown","source":"## Import Statements","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nfrom PIL import Image\nfrom pathlib import Path\nimport gc\nfrom tqdm import tqdm\nfrom concurrent.futures import ProcessPoolExecutor, as_completed","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Setting global variables","metadata":{}},{"cell_type":"code","source":"DATA_ROOT = \"/kaggle/input/rsna-2023-abdominal-trauma-detection\"\nTRAIN_DIR = DATA_ROOT + \"/train_images\"\nTEST_DIR = DATA_ROOT + \"/test_images\"\nTRAIN_CSV = DATA_ROOT + \"/train.csv\"\nTRAIN_TAGS = DATA_ROOT + \"/train_dicom_tags.parquet\"\nTEST_TAGS = DATA_ROOT + \"/test_dicom_tags.parquet\"\nSEGMEN_DIR = DATA_ROOT + \"/segmentations\"\nSERIES_CSV = DATA_ROOT + \"/train_series_meta.csv\"\nSLICES_CSV = DATA_ROOT + \"/image_level_labels.csv\"\n###############################################################################\nOUT_IMG_DIR = \"/kaggle/working/images\"\nOUT_SIZE = (128, 128)\nOUT_EXT = 'png'","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_mappers():\n    series2patient, patient2series = {}, {}\n    for patient_id in Path(TRAIN_DIR).glob(\"*\"):\n        series_ids = (Path(TRAIN_DIR) / patient_id).glob(\"*\")\n        patient_id = int(patient_id.name)\n        series_ids = sorted([int(sid.name) for sid in series_ids])\n        patient2series[patient_id] = series_ids\n        for series_id in series_ids:\n            series2patient[series_id] = patient_id\n    return series2patient, patient2series","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Map SeriesID to PatientID and PatientID to SeriesIDs\nseries2patient, patient2series = create_mappers()\nseries_ids = sorted(series2patient.keys())\npatient_ids = sorted(patient2series.keys())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Extracting images of healthy individuals\nIn this context, healthy means not having any abdominal injury targeted in the competition (i.e. `bowel`, `extravasation`, `kidney`, `liver` and `spleen`).","metadata":{}},{"cell_type":"code","source":"def get_info_from_tags(parquet_path=TRAIN_TAGS):\n    \n    dcmtags = pd.read_parquet(parquet_path)\n    spatial_dims = ['X', 'Y', 'Z']\n    uid_tags = ['SeriesInstanceUID']\n    intscalar_tags = ['PatientID', 'Rows', 'Columns']\n    floattuple_tags = ['PixelSpacing', 'ImagePositionPatient']\n    selected_tags = intscalar_tags\n\n    # Extract SeriesID from SeriesInstanceUID\n    dcmtags['SeriesID'] = dcmtags['SeriesInstanceUID'].str.split('.').str[-1].astype(int)\n    \n    # Transfrom Integer tags into integers\n    for col in intscalar_tags:\n        dcmtags[col] = dcmtags[col].astype(int)\n\n    # Transfrom FloatTuple tags into seperate float columns\n    for col in floattuple_tags:\n        tmp_series = dcmtags[col].str[1:-1].str.split(',', expand=True).astype(float)\n        for idx, d in enumerate(spatial_dims[:tmp_series.shape[-1]]):\n            dcmtags[col+d] = tmp_series.iloc[:, idx]\n            selected_tags += [col+d]\n    \n    # Rename columns\n    cols_dict = {'Rows': 'Height', 'Columns': 'Width'}\n    dcmtags.rename(columns=cols_dict, inplace=True)\n    selected_tags = list(map(cols_dict.get, selected_tags, selected_tags))\n    \n    # Create two columns for PixelSpacingZ and Depth\n    zcol = 'ImagePositionPatientZ'\n    zspacing = dcmtags[['SeriesID', zcol]].sort_values(zcol)\n    zspacing = zspacing.groupby('SeriesID').diff(axis='index')[zcol]\n    zspacing = zspacing.groupby(level=0).mean().abs()\n    depth = dcmtags.groupby('SeriesID').count()['PatientID']\n    dcmtags = dcmtags.groupby('SeriesID').max()\n    dcmtags['Depth'] = depth\n    dcmtags['PixelSpacingZ'] = zspacing.round(6)\n    selected_tags += ['Depth', 'PixelSpacingZ']\n    \n    # Return only Spatial Information\n    dcmtags = dcmtags[selected_tags]\n\n    return dcmtags","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get SeriesIDs corresponding to healthy 512x512xDepth volumes \nlabels_df = pd.read_csv(TRAIN_CSV)\nlabels_df['series_id'] = labels_df['patient_id'].apply(patient2series.get)\nlabels_df = labels_df.explode('series_id', ignore_index=True)\nlabels_df.set_index('series_id', inplace=True)\nlabels_df = labels_df[['patient_id', 'any_injury']]\n\ntrain_dcmtags = get_info_from_tags()\nlabels_df = labels_df.join(train_dcmtags)\n\nlabels_df = labels_df[~labels_df.any_injury.astype(bool)]\nlabels_df = labels_df[labels_df.Width == 512]\nhealthy_ids = sorted(labels_df.index)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Helper functions for Preprocessing DICOMs","metadata":{}},{"cell_type":"code","source":"def correct_uint14(x):\n    vmax_uint14 = 16383\n    if x.max() == vmax_uint14:\n        x[x > vmax_uint14 // 2] = x[x > vmax_uint14 // 2] - (vmax_uint14 + 1)\n    return x\n\ndef load_dicom_slice(path, hounsfield=True):\n    dcm = pydicom.dcmread(path)\n    img = dcm.pixel_array\n    # Convert to Hounsfield Units\n    if hounsfield:\n        slope = int(dcm.RescaleSlope)\n        intercept = int(dcm.RescaleIntercept)\n        img = slope * img + intercept\n    img = correct_uint14(img)\n    return img\n\ndef make_uint8(x, vmin, vmax):\n    # Normalize\n    x = (x - vmin) / (vmax - vmin)\n    # Rescale\n    x = x * 255\n    # Cast\n    x = x.astype('uint8')\n    return x","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id = \"core\"></a>\n# Transforming DICOM files using Multiprocessing","metadata":{}},{"cell_type":"code","source":"!tree {TRAIN_DIR}/10004/21057 | head -n 10 # PatientID/SeriesID/InstanceID.dcm","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Instead of evenly sampling slices from the 3D volumes, we oversample from the top and the bottom in order to have more instances of arms and legs. This will help us locate the torso region in a single axial slice and in the 3D volume.","metadata":{}},{"cell_type":"code","source":"def even_sampling(slices, trim_depth):\n    depth = len(slices)\n    if trim_depth and depth//trim_depth:\n        slices = slices[::depth//trim_depth]\n        slices = slices[:trim_depth]\n    return slices\n\ndef custom_sampling(slices, trim_depth):\n    depth = len(slices)   \n    if trim_depth and depth//trim_depth > 1 :\n        top_section = slices[:depth//4]\n        middle_section = slices[depth//4:3*depth//4]\n        bottom_section = slices[3*depth//4:]\n        n_top = n_bottom = 7 * trim_depth // 16\n        n_middle = trim_depth // 8\n        top_section = top_section[::len(top_section) // n_top]\n        middle_section = middle_section[::len(middle_section) // n_middle]\n        bottom_section = bottom_section[::len(bottom_section) // n_top]\n        top_section = top_section[:n_top]\n        middle_section = middle_section[:n_middle]\n        bottom_section = bottom_section[:n_bottom]\n        slices = top_section + middle_section + bottom_section\n    return slices","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_slices(series_id, trim_depth=64):\n    patient_id = series2patient[series_id]\n    slices = (Path(TRAIN_DIR) / str(patient_id) / str(series_id)).glob(\"*.dcm\") \n    slices = sorted(slices, key=lambda p: int(p.stem))\n    slices = custom_sampling(slices, trim_depth)\n    return slices\n\ndef get_slices_concurrent(series_ids):\n    with ProcessPoolExecutor() as pool:\n        jobs = [pool.submit(get_slices, sid) for sid in series_ids]\n        slices = []\n        for job in tqdm(as_completed(jobs),\n                        total=len(jobs),\n                        desc='Searching for DICOM files...'):\n            slices.extend(job.result())\n    gc.collect()\n    return slices","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def save_axial(path):\n    patient_id = int(path.parent.parent.stem)\n    series_id = int(path.parent.stem)\n    instance_id = int(path.stem)\n    fname = f\"{patient_id:06d}-{series_id:06d}-{instance_id:06d}.{OUT_EXT}\"\n    save_dir = Path(OUT_IMG_DIR)\n    save_dir.mkdir(parents=True, exist_ok=True)\n    img = load_dicom_slice(path)\n    vmin, vmax = -1024, 1023\n    img = np.clip(img, vmin, vmax)\n    img = make_uint8(img, vmin, vmax)\n    pil_img = Image.fromarray(img)\n    pil_img = pil_img.resize(OUT_SIZE, Image.LANCZOS)\n    pil_img.save(save_dir / fname)\n    del img, pil_img\n    gc.collect()\n    \ndef save_axial_concurrent(paths):\n    with ProcessPoolExecutor() as pool:\n        jobs = [pool.submit(save_axial, path) for path in paths]\n        for _ in tqdm(as_completed(jobs), \n                      total=len(jobs),\n                      desc='Processing axial slices...'):\n            pass\n    gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices = get_slices_concurrent(healthy_ids)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"save_axial_concurrent(slices)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id = \"result\"></a>\n# Samples","metadata":{}},{"cell_type":"code","source":"def stack_3dvolume(img, trim_depth=0, depth_last=True, nrows=4):\n    \"\"\"Convert a 3D volume into 2D image\"\"\"\n    assert img.ndim == 3, f\"Expected 3D-array, got f{img.ndim}-array\"\n    if depth_last:\n        img = np.transpose(img, (2, 0, 1))\n    d, h, w = img.shape\n    if trim_depth and d//trim_depth:\n        img = img[::d//trim_depth]\n        img = img[:trim_depth]\n        d = img.shape[0]\n    assert d % nrows == 0, \"Depth should be a multiple of nrows\"\n    img = img.reshape(nrows, d//nrows, h, w).swapaxes(1, 2)\n    img = img.reshape(nrows*h, (d//nrows)*w)    \n    return img","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nrows, ncols = 10, 10\nsample_generator = Path(OUT_IMG_DIR).glob(\"*\")\nimgs = []\nfor idx in range(nrows*ncols):\n    img = Image.open(next(sample_generator))\n    imgs.append(np.array(img))\nimg_grid = np.stack(imgs, axis=-1)\nimg_grid = stack_3dvolume(img_grid, nrows=nrows)\nImage.fromarray(img_grid)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}