{"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":"# RSNA Fracture Detection - DICOM to 3D Numpy Volume\n\n![CT Image](https://storage.googleapis.com/kaggle-datasets-images/674071/1185670/0d449dc88ac1318321ae8f7de974f0fa/dataset-cover.jpg?t=2020-05-25-13-21-33)\n\n- **Author:** *Mariusz Wiśniewski*\n- **Competition:** *[RSNA 2022 Cervical Spine Fracture Detection](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection)*\n- **Date created:** *November 27th, 2022*\n- **Last modified:** *November 27th, 2022*\n\n## Overview\n\nIn this notebook we will prepare 3D numpy volumes from DICOM files.\n\n### Libraries Used\n\n- [SciPy 🔬](https://scipy.org)\n- [OpenCV 🖼️](https://opencv.org)\n- [Pydicom 🩻](https://pydicom.github.io)\n\n### References\n\n- [Uniformizing Techniques to Process CT scans with 3D CNNs for Tuberculosis Prediction 📃](https://arxiv.org/abs/2007.13224)\n- [Grad-CAM: Visual Explanations from Deep Networks via Gradient-based Localization 📃](https://arxiv.org/abs/1610.02391)\n- [3D image classification from CT scans 📝](https://keras.io/examples/vision/3D_image_classification/)\n- [A Comprehensive Introduction to Different Types of Convolutions in Deep Learning](https://towardsdatascience.com/a-comprehensive-introduction-to-different-types-of-convolutions-in-deep-learning-669281e58215)\n- [[RSNA_22] Dicom to NumPy 3D 📓](https://www.kaggle.com/code/vmuzhichenko/rsna-22-dicom-to-numpy-3d)","metadata":{}},{"cell_type":"markdown","source":"# Notebook Steup","metadata":{}},{"cell_type":"markdown","source":"## Required Packages","metadata":{}},{"cell_type":"code","source":"# install pydicom requirements\n!conda install '/kaggle/input/pydicom-conda-helper/libjpeg-turbo-2.1.0-h7f98852_0.tar.bz2' --offline -yq\n!conda install '/kaggle/input/pydicom-conda-helper/libgcc-ng-9.3.0-h2828fa1_19.tar.bz2' --offline -yq\n!cp ../input/gdcm-conda-install/gdcm.tar .\n!tar -xvzf gdcm.tar\n!conda install --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2 -q\n!conda install '/kaggle/input/pydicom-conda-helper/conda-4.10.1-py37h89c1867_0.tar.bz2' --offline -yq\n!conda install '/kaggle/input/pydicom-conda-helper/certifi-2020.12.5-py37h89c1867_1.tar.bz2' --offline -yq\n!conda install '/kaggle/input/pydicom-conda-helper/openssl-1.1.1k-h7f98852_0.tar.bz2' --offline -yq\n!rm -rf gdcm/ gdcm.tar","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":false,"_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2022-11-29T00:30:11.60983Z","iopub.execute_input":"2022-11-29T00:30:11.611357Z","iopub.status.idle":"2022-11-29T00:31:15.035948Z","shell.execute_reply.started":"2022-11-29T00:30:11.611236Z","shell.execute_reply":"2022-11-29T00:31:15.034294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Import Statements","metadata":{}},{"cell_type":"code","source":"import gc\nimport glob\nimport os\n\nimport cv2\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom joblib import Parallel, delayed\nfrom matplotlib import animation, rc\nfrom scipy.ndimage import zoom\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.039047Z","iopub.execute_input":"2022-11-29T00:31:15.039555Z","iopub.status.idle":"2022-11-29T00:31:15.740944Z","shell.execute_reply.started":"2022-11-29T00:31:15.039509Z","shell.execute_reply":"2022-11-29T00:31:15.739978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Project Configuration","metadata":{}},{"cell_type":"code","source":"config = {\n    'img_size': 256,\n    'depth': 128,\n    'ess': False,\n    'hu': True\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.744584Z","iopub.execute_input":"2022-11-29T00:31:15.745247Z","iopub.status.idle":"2022-11-29T00:31:15.750506Z","shell.execute_reply.started":"2022-11-29T00:31:15.745212Z","shell.execute_reply":"2022-11-29T00:31:15.748996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Processing DICOM Files","metadata":{}},{"cell_type":"markdown","source":"## Loading Data and Preprocessing","metadata":{}},{"cell_type":"code","source":"TRAIN_IMG_PATH = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/'\nTEST_IMG_PATH = '../input/rsna-2022-cervical-spine-fracture-detection/test_images/'\nTRAIN_CSV_PATH = '../input/rsna-2022-cervical-spine-fracture-detection/train.csv'\nTEST_CSV_PATH = '../input/rsna-2022-cervical-spine-fracture-detection/test.csv'\nTRAIN_BBOX_CSV_PATH = (\n    '../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv'\n)\n\ntrain_images = os.listdir(TRAIN_IMG_PATH)\ntest_images = os.listdir(TEST_IMG_PATH)\n\ntrain = pd.read_csv(TRAIN_CSV_PATH)\ntest = pd.read_csv(TEST_CSV_PATH)\ntrain_bbox = pd.read_csv(TRAIN_BBOX_CSV_PATH)","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.753194Z","iopub.execute_input":"2022-11-29T00:31:15.753526Z","iopub.status.idle":"2022-11-29T00:31:15.88628Z","shell.execute_reply.started":"2022-11-29T00:31:15.753495Z","shell.execute_reply":"2022-11-29T00:31:15.88533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Creating Output Directories","metadata":{}},{"cell_type":"code","source":"TRAIN_OUTPUT_PATH = (\n    f'./train_volumes_{config[\"img_size\"]}x{config[\"img_size\"]}x{config[\"depth\"]}/'\n)\nTEST_OUTPUT_PATH = (\n    f'./test_volumes_{config[\"img_size\"]}x{config[\"img_size\"]}x{config[\"depth\"]}/'\n)\n\nif not os.path.exists(TRAIN_OUTPUT_PATH):\n    os.mkdir(TRAIN_OUTPUT_PATH)\nif not os.path.exists(TEST_OUTPUT_PATH):\n    os.mkdir(TEST_OUTPUT_PATH)","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.887532Z","iopub.execute_input":"2022-11-29T00:31:15.887902Z","iopub.status.idle":"2022-11-29T00:31:15.895025Z","shell.execute_reply.started":"2022-11-29T00:31:15.887871Z","shell.execute_reply":"2022-11-29T00:31:15.893883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocessing Functions\n\nHere we define several helper functions to process the data. It is valuable to note that the `zoom` function used to resize the volume employs spline interpolation of the requested order. According to [this paper](https://arxiv.org/abs/2007.13224), this operation is known as *spline interpolated zoom (SIZ)*. Additionally, we add a flag to perform *Even Slice Selection (ESS)* instead.\n\nCT scans store raw voxel intensity in *Hounsfield units (HU)*. They range from -2000 to above 3000 in this dataset. Because bones with varying radio intensities exist over 400, this is utilized as an upper bound. A threshold between -1000 and 400 is commonly used to normalize CT scans.","metadata":{}},{"cell_type":"code","source":"def load_dicom_slice(path, ess=False, hu=True):\n    \"\"\"Read and load a single slice\"\"\"\n    dicom = pydicom.read_file(path)\n    data = dicom.pixel_array\n\n    if ess:\n        data = cv2.resize(\n            data, (config['img_size'], config['img_size']), interpolation=cv2.INTER_AREA\n        )\n        \n    if hu:\n        # convert to Hounsfield units (HU)\n        data = data * dicom.RescaleSlope + dicom.RescaleIntercept\n    return np.array(data).astype(np.float32)\n\n\ndef load_dicom_volume(path, ess=False, hu=True):\n    \"\"\"Load dicom files corresponding to one volume in parallel\"\"\"\n    t_paths = sorted(\n        glob.glob(os.path.join(path, '*')),\n        key=lambda x: int(x.split('/')[-1].split('.')[0]),\n    )\n\n    if ess:\n        # instead of zooming whole dicom series, perform Even Slice Selection (ESS)\n        n_scans = len(os.listdir(path))\n        indices = (\n            np.quantile(list(range(n_scans)), np.linspace(0.0, 1.0, config['depth']))\n            .round()\n            .astype(int)\n        )\n        t_paths = [t_paths[i] for i in indices]\n\n    images = Parallel(n_jobs=-1)(\n        delayed(load_dicom_slice)(filename, ess, hu) for filename in t_paths\n    )\n    volume = np.array(images)\n    # transpose volume from (z, x, y) to (x, y, z)\n    volume = np.transpose(volume, (1, 2, 0))\n    return volume\n\n\ndef clip_hu(volume, min_hu=-1000, max_hu=400):\n    \"\"\"Clip bones (≥400) and air (≤-1000) HU\"\"\"\n    volume[volume < min_hu] = min_hu\n    volume[volume > max_hu] = max_hu\n    return volume\n\n\ndef normalize(volume, hu=True):\n    \"\"\"Normalize the volume\"\"\"\n    if hu:\n        volume = clip_hu(volume)\n    volume = (volume - np.min(volume)) / (np.max(volume) - np.min(volume))\n    # using uint8 to save memory\n    volume = (volume * 255).astype(np.uint8)\n    return volume\n\n\ndef resize_volume(img, desired_width=128, desired_height=128, desired_depth=64):\n    \"\"\"Resize the volume\"\"\"\n    # Compute zoom factors\n    width_factor = desired_width / img.shape[0]\n    height_factor = desired_height / img.shape[1]\n    depth_factor = desired_depth / img.shape[-1]\n\n    # Resize the volume using spline interpolated zoom (SIZ)\n    img = zoom(img, (width_factor, height_factor, depth_factor), order=1)\n    return img\n\n\ndef process_scan(path, ess=False, hu=True):\n    \"\"\"Read and resize a volume\"\"\"\n    # Read scan\n    volume = load_dicom_volume(path, ess, hu)\n    # Normalize\n    volume = normalize(volume, hu)\n    if not ess:\n        # Resize width, height, and depth\n        volume = resize_volume(\n            volume, config['img_size'], config['img_size'], config['depth']\n        )\n    return volume\n\n\ndef save_3d_voxels(dicom_path, output_path, ess=False, hu=True):\n    \"\"\"Save volume as .npy\"\"\"\n    image = process_scan(dicom_path, ess=ess, hu=hu)\n    np.savez_compressed(f'{output_path}{dicom_path.split(\"/\")[-1]}.npz', data=image)\n\n    del image","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.896424Z","iopub.execute_input":"2022-11-29T00:31:15.896929Z","iopub.status.idle":"2022-11-29T00:31:15.916762Z","shell.execute_reply.started":"2022-11-29T00:31:15.896897Z","shell.execute_reply":"2022-11-29T00:31:15.915646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Separating Data with Bounding Box Annotations from the Training Data\n\nAs you can see, the study cases covered by the training dataset collide with the study cases with bounding box annotations provided. Since we want to develop a weekly-supervised model for object detection, we will use the latter study cases for testing only.","metadata":{}},{"cell_type":"code","source":"common = [\n    t\n    for t in set(train_bbox['StudyInstanceUID'].values)\n    if t in set(train['StudyInstanceUID'].values)\n]\nprint(f'There are {len(common)} study cases that are in both train and train_bbox.')","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.918271Z","iopub.execute_input":"2022-11-29T00:31:15.918604Z","iopub.status.idle":"2022-11-29T00:31:15.969488Z","shell.execute_reply.started":"2022-11-29T00:31:15.918576Z","shell.execute_reply":"2022-11-29T00:31:15.968672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# use study cases with bbox annotations for testing\ntest_patients = common\ntrain_patients = train.StudyInstanceUID.to_list()\n# remove test_patients' elements from train_patients\ntrain_patients = [patient for patient in train_patients if patient not in test_patients]","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.97075Z","iopub.execute_input":"2022-11-29T00:31:15.971412Z","iopub.status.idle":"2022-11-29T00:31:15.985154Z","shell.execute_reply.started":"2022-11-29T00:31:15.971379Z","shell.execute_reply":"2022-11-29T00:31:15.983913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualize HU Clipping","metadata":{}},{"cell_type":"code","source":"rc('animation', html='jshtml')\n\n\ndef create_double_animation(array, array_hu, case):\n    \"\"\"Create animation of two volumes\"\"\"\n    # transpose volume from (x, y, z) to (z, x, y)\n    array = np.transpose(array, (2, 0, 1))\n    array_hu = np.transpose(array_hu, (2, 0, 1))\n\n    fig = plt.figure(figsize=(12, 6))\n    fig.tight_layout()\n    plt.axis('off')\n    plt.suptitle(f'Patient ID: {case}', fontsize=16, fontweight='bold')\n\n    ax1 = fig.add_subplot(1,2,1)\n    ax1.set_title('HU Clipping Off')\n    ax2 = fig.add_subplot(1,2,2)\n    ax2.set_title('HU Clipping On')\n\n    images = []\n    for image, image_hu in zip(array, array_hu):\n        im1 = ax1.imshow(image, animated=True, cmap='bone')\n        im2 = ax2.imshow(image_hu, animated=True, cmap='bone')\n        images.append([im1, im2])\n\n    ani = animation.ArtistAnimation(\n        fig, images, interval=5000 // len(array), blit=False, repeat_delay=1000\n    )\n    plt.close()\n    return ani","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:15.986627Z","iopub.execute_input":"2022-11-29T00:31:15.987045Z","iopub.status.idle":"2022-11-29T00:31:15.998139Z","shell.execute_reply.started":"2022-11-29T00:31:15.986974Z","shell.execute_reply":"2022-11-29T00:31:15.997064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_patient = train_patients[0]\nsample_path = f'{TRAIN_IMG_PATH}{sample_patient}'\nimage = process_scan(sample_path, ess=False, hu=False)\nimage_hu = process_scan(sample_path, ess=False, hu=True)\n\ncreate_double_animation(image, image_hu, sample_patient)","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:16.002182Z","iopub.execute_input":"2022-11-29T00:31:16.00257Z","iopub.status.idle":"2022-11-29T00:31:38.834463Z","shell.execute_reply.started":"2022-11-29T00:31:16.002538Z","shell.execute_reply":"2022-11-29T00:31:38.833417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Saving 3D Numpy Volumes","metadata":{}},{"cell_type":"code","source":"for i in tqdm(range(len(train_patients))):\n    case = TRAIN_IMG_PATH + train_patients[i]\n    save_3d_voxels(case, TRAIN_OUTPUT_PATH, ess=config['ess'], hu=config['hu'])\n\ngc.collect()\n\nfor i in tqdm(range(len(test_patients))):\n    case = TRAIN_IMG_PATH + test_patients[i]\n    save_3d_voxels(case, TEST_OUTPUT_PATH, ess=config['ess'], hu=config['hu'])\n    \ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-11-29T00:31:38.836056Z","iopub.execute_input":"2022-11-29T00:31:38.836584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization","metadata":{}},{"cell_type":"code","source":"rc('animation', html='jshtml')\n\n\ndef create_animation(array, case):\n    \"\"\"Create animation of a volume\"\"\"\n    # transpose volume from (x, y, z) to (z, x, y)\n    array = np.transpose(array, (2, 0, 1))\n\n    fig = plt.figure(figsize=(6, 6))\n    images = []\n    for image in array:\n        image_plot = plt.imshow(image, animated=True, cmap='bone')\n        plt.axis('off')\n        images.append([image_plot])\n\n    plt.title(f'Patient ID: {case}', fontsize=16)\n    ani = animation.ArtistAnimation(\n        fig, images, interval=5000 // len(array), blit=False, repeat_delay=1000\n    )\n    plt.close()\n    return ani","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Animation\n\nHere we will create an animation of the whole 3D volume corresponding to a single study case.","metadata":{}},{"cell_type":"code","source":"train_volumes = os.listdir(TRAIN_OUTPUT_PATH)\nvolume_path = TRAIN_OUTPUT_PATH + train_volumes[np.random.randint(len(train_volumes))]\nvolume = np.load(volume_path)\nvolume = volume['data']\n\nprint(f'Volume shape: {volume.shape}')\ncase_id = volume_path.split('/')[-1][:-4]\ncreate_animation(volume, case_id)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Zip Output","metadata":{}},{"cell_type":"code","source":"!tar -czf output.tar.gz --remove-files ./*","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}