{"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":"In this notebook, we'll **consolidate series metadata** from different files, perform some basic exploratory data analysis, and **transform 2D DICOM files into 3D NIfTI arrays**. This initial step is a prerequisite before progressing to further data processing and model feeding. But the latter stage of the project is still quite distant from where we are now :)\n\n**Important note**! The transformation process takes approximately 24 hours to complete. I will upload the results to Google Drive and share the link on the discussion forum shortly. Please refrain from running this code on Kaggle - the transformed 3D NIfTI arrays have almost the same size as the original DICOM files and you'll run out of disc space well before reaching the finish line.","metadata":{}},{"cell_type":"code","source":"pip install dicom2nifti","metadata":{"execution":{"iopub.status.busy":"2023-08-09T08:36:00.100372Z","iopub.execute_input":"2023-08-09T08:36:00.100838Z","iopub.status.idle":"2023-08-09T08:36:12.845306Z","shell.execute_reply.started":"2023-08-09T08:36:00.100806Z","shell.execute_reply":"2023-08-09T08:36:12.84372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import datetime\nimport os\nimport shutil\nimport sys\nimport numpy as np\nimport pandas as pd\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nimport pydicom\nimport nibabel as nib\nfrom nibabel import processing\nimport dicom2nifti\nfrom typing import Generator","metadata":{"execution":{"iopub.status.busy":"2023-08-09T08:36:12.848142Z","iopub.execute_input":"2023-08-09T08:36:12.848592Z","iopub.status.idle":"2023-08-09T08:36:12.869288Z","shell.execute_reply.started":"2023-08-09T08:36:12.848553Z","shell.execute_reply":"2023-08-09T08:36:12.867697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before we dive in, the function below sets all the parameters that are used across the functions within the scope of the project. It might seem a bit mundane at the start, but as the codebase grows to hundreds of lines across dozens of different functions, it becomes extremely useful in maintaining control. It also simplifies your life when transitioning between your local machine, Google Colab, and Kaggle. Currently, the function has only two parameters, but it will gradually expand as we continue to develop the project.","metadata":{}},{"cell_type":"code","source":"def set_parameters(root: str, working: str) -> dict:\n    \"\"\"\n    Set uniform parameters for the functions in the scope of the project.\n    :param root: Path to the project directory, where 'train_images' and 'test_images' directories are located,\n    :param working: Path to the working directory,\n    :return: Dictionary of uniform parameters.\n    \"\"\"\n\n    # Get dictionary of uniform parameters\n    parameters = {\n        'datasets': {\n            'train': {\n                'path': os.path.join(root, 'train_images'),\n                'tags': os.path.join(root, 'train_dicom_tags.parquet'),\n                'meta': os.path.join(root, 'train_series_meta.csv'),\n                'patient_labels': os.path.join(root, 'train.csv'),\n                'image_labels': os.path.join(root, 'image_level_labels.csv'),\n                'segmentation': os.path.join(root, 'segmentations/')},\n            'test': {\n                'path': os.path.join(root, 'test_images/'),\n                'tags': os.path.join(root, 'test_dicom_tags.parquet'),\n                'meta': os.path.join(root, 'test_series_meta.csv')},\n            'totalsegmentator': {\n                'path': os.path.join(root, 'Totalsegmentator_dataset/'),\n                'meta': os.path.join(root, 'Totalsegmentator_dataset/meta.csv')}},\n        'root': root,\n        'working': working}\n\n    return parameters","metadata":{"execution":{"iopub.status.busy":"2023-08-09T08:38:53.486441Z","iopub.execute_input":"2023-08-09T08:38:53.486918Z","iopub.status.idle":"2023-08-09T08:38:53.497653Z","shell.execute_reply.started":"2023-08-09T08:38:53.486884Z","shell.execute_reply":"2023-08-09T08:38:53.496273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"parameters = set_parameters(root='/kaggle/input/rsna-2023-abdominal-trauma-detection', working='/kaggle/working/')","metadata":{"execution":{"iopub.status.busy":"2023-08-09T08:39:26.713255Z","iopub.execute_input":"2023-08-09T08:39:26.713721Z","iopub.status.idle":"2023-08-09T08:39:26.721118Z","shell.execute_reply.started":"2023-08-09T08:39:26.713686Z","shell.execute_reply":"2023-08-09T08:39:26.719543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, let's consolidate all series metadata into a single file. We have two primary sources for this consolidation:\n1. 'train_series_meta.csv' with patient- and series IDs, Hounsfield Unit, and incomplete organ indicator - 4711 lines; \n2. 'train_dicom_tags.parquet' with image-level attributes extracted from DICOM files, encompassing several dozen attributes - 1.5 million lines. \n\nSince all image attributes (except for ImagePositionPatient) are specific to series, we'll combine them into a unified series metadata file. Working with 4.7 thousand instead of 1.5 million lines should make the life easier.","metadata":{}},{"cell_type":"code","source":"def get_metadata(dataset: str, **kwargs) -> pd.DataFrame:\n    \"\"\"\n    Add extra attributes to series metadata.\n\n    :param dataset:  Dataset (train, test).\n    :param kwargs: Dictionary of uniform parameters.\n    :return: pd.DataFrame.\n    \"\"\"\n\n    # Open previously saved metadata.csv, if any\n    path = os.path.join(kwargs['working'], 'metadata.csv')\n    if os.path.exists(path):\n        meta = pd.read_csv(path)\n        cols = ['ImageOrientationPatient', 'PixelSpacing']\n        for col in cols:\n            meta[col] = meta[col].apply(lambda x: tuple(round(float(value), 4) for value in x.strip('()').split(',')))\n\n    # If not - add extra attributes to series metadata\n    else:\n\n        # Read series metadata and image tags\n        cols = ['PatientID','SeriesID','AorticHU','IncompleteOrgan']\n        meta = pd.read_csv(kwargs['datasets'][dataset]['meta'], names=cols, skiprows=1)\n        tags = pd.read_parquet(kwargs['datasets'][dataset]['tags'])\n\n        # Drop missing series from tags\n        tags['SeriesID'] = tags['SeriesInstanceUID'].apply(lambda x: int(x.split('.')[-1]))\n        tags = tags[tags['SeriesID'].isin(meta['SeriesID'].unique())]\n\n        # # Add index of first and last image and image count per series\n        tmp = tags.groupby(by='SeriesID')['InstanceNumber'].agg(['min', 'max', 'count'])\n        cols = {'min': 'ImageFirst', 'max': 'ImageLast', 'count': 'ImageCount'}\n        tmp.rename(columns=cols, inplace=True)\n        meta = meta.merge(right=tmp, left_on='SeriesID', right_index=True, how='left')\n\n        # Convert image attributes to tuples of floats, where applicable\n        cols = ['ImageOrientationPatient', 'PixelSpacing']\n        for col in cols:\n            tags[col] = tags[col].apply(lambda x: tuple(round(float(value), 4) for value in x.strip('[]').split(',')))\n\n        # Add image attributes to metadata\n        cols = ['SliceThickness', 'KVP', 'PatientPosition', 'ImageOrientationPatient',\n                'SamplesPerPixel', 'PhotometricInterpretation', 'Rows', 'Columns', 'PixelSpacing',\n                'BitsAllocated', 'BitsStored', 'HighBit', 'PixelRepresentation',\n                'WindowCenter', 'WindowWidth', 'RescaleIntercept', 'RescaleSlope', 'RescaleType']\n        for col in cols:\n            meta[col] = meta['SeriesID'].apply(lambda x: tags[tags['SeriesID'] == x][col].unique()[0])\n\n        # Add 'Isotropic' (square pixel) and 'Square' (square image) attributes\n        meta['Isotropic'] = meta['PixelSpacing'].apply(lambda x: 1 if x[0] == x[1] else 0)\n        meta['Square'] = meta.apply(lambda x: 1 if x['Rows'] == x['Columns'] else 0, axis=1)\n\n        # Convert 'PatientPosition' {HFS: 1, FFS: 0} and 'RescaleType' {HU: 1, None: 0} to binary\n        meta['PatientPosition'] = meta['PatientPosition'].apply(lambda x: 1 if x == 'HFS' else 0)\n        meta['RescaleType'] = meta['RescaleType'].apply(lambda x: 1 if x in ['HU', 'Houndsfield Unit'] else 0)\n\n        # Write to metadata.csv\n        meta.to_csv(os.path.join(kwargs['working'], 'metadata.csv'), index=False)\n\n    return meta","metadata":{"execution":{"iopub.status.busy":"2023-08-09T09:16:23.749064Z","iopub.execute_input":"2023-08-09T09:16:23.749608Z","iopub.status.idle":"2023-08-09T09:16:23.768778Z","shell.execute_reply.started":"2023-08-09T09:16:23.749572Z","shell.execute_reply":"2023-08-09T09:16:23.766723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta = get_metadata(dataset='train', **parameters)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T09:16:33.725996Z","iopub.execute_input":"2023-08-09T09:16:33.726432Z","iopub.status.idle":"2023-08-09T09:20:48.88581Z","shell.execute_reply.started":"2023-08-09T09:16:33.7264Z","shell.execute_reply":"2023-08-09T09:20:48.884715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now part for some EDA.Since I'm planning to use 3D images, I'm particulary interested in image and pixel size, dimensions, etc. I will not go into class distribution and label analysis at this stage, so it will not take much time, I promice :)","metadata":{}},{"cell_type":"code","source":"def eda(metadata) -> None:\n    \"\"\"\n    Perform exploratory data analysis.\n    :param metadata: metadata pd.Dataframe\n    :return: None\n    \"\"\"\n\n    # Set visual patterns\n    sns.set_theme(style=\"darkgrid\")\n\n    # Image and pixel dimensions\n    fig, ax = plt.subplots(2)\n    dist = pd.crosstab(meta['Square'], meta['Isotropic'])\n    cmap = sns.color_palette(\"light:b\", as_cmap=True)\n    sns.heatmap(dist, annot=True, fmt=\"d\", cmap=cmap, cbar=False, xticklabels=['Anisotropic', 'Isotropic'],\n                yticklabels=['Non-square', 'Square'], ax=ax[0], linewidths=1)\n    tmp = meta[meta['Isotropic'] == 1]['PixelSpacing'].apply(lambda x: x[0])\n    ax[1] = sns.histplot(data=tmp, color=cmap(0.95), ax=ax[1])\n    ax[0].set(xlabel=\"\", ylabel=\"\")\n    ax[1].set(xlabel='Pixel size, mm', ylabel='Number of series', xlim=(0.4, 1.1))\n    plt.suptitle('Distribution of series by image/pixel size', size=16)\n    plt.show()\n\n    # Slice thickness and image count per series\n    fig, ax = plt.subplots(nrows=2, ncols=1, sharex='all')\n    tmp = meta['SliceThickness'].round(1).value_counts().reset_index()\n    sns.barplot(data=tmp, x='index', y='SliceThickness', color=cmap(0.95), ax=ax[0])\n    sns.boxplot(x=meta['SliceThickness'].round(1), y=meta['ImageCount'], color=cmap(0.95), whis=[0,100], ax=ax[1])\n    ax[0].set(xlabel=None, ylabel='Number of series', ylim=(0,2000))\n    ax[1].set(xlabel='Slice thickness, mm', ylabel='Images per series', ylim=(0,1750))\n    plt.suptitle(t='Slice thickness and image count per series', size=16)\n    plt.show()\n\n    return None","metadata":{"execution":{"iopub.status.busy":"2023-08-09T09:21:15.80866Z","iopub.execute_input":"2023-08-09T09:21:15.809052Z","iopub.status.idle":"2023-08-09T09:21:15.822089Z","shell.execute_reply.started":"2023-08-09T09:21:15.809021Z","shell.execute_reply":"2023-08-09T09:21:15.820663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Just a few short observations:\n1. Most of the series consist of square images (512x512 pixels) with isotropic pixel spacing (pixels with equal x- and y- axis dimensions). 166 series have anisotropic pixel spasing and (or) non-square images. We might want to handle that at preprocessing stage.\n2. Pixel sizes vary a lot across series - from 0.4 to 1.2 mm. I'm not sure if this would be a problem, will leave just as a side note for now.\n3. Slice thickness and numer of images per series have even larger variation: from 0.5 to 5 mm and from 41 to 1727 images per series. Looks like if we want to feed 3D images to the model (which is exactly the plan I have in mind) - we need to resample the series at the preprocessing stage.","metadata":{}},{"cell_type":"code","source":"eda(metadata=meta)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T09:21:22.222228Z","iopub.execute_input":"2023-08-09T09:21:22.222634Z","iopub.status.idle":"2023-08-09T09:21:23.724068Z","shell.execute_reply.started":"2023-08-09T09:21:22.222603Z","shell.execute_reply":"2023-08-09T09:21:23.722301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now the last part of the puzzle for today: a function to stack images in 3D arrays and convert from DICOM to NIfTI format. The next stop in a few days would be to preprocess 3D images (resampling, resizing, cropping, etc.), combine them with a Totalsegmentation dataset, and train the model to detect abdominal organs on the individual images. This should help to move from patient-level (which seem to be pretty weak) to image-level labels. Let me know if you want me to share the results, if there would be any :)","metadata":{}},{"cell_type":"code","source":"def dcm2nii(dataset: str, **kwargs) -> None:\n    \"\"\"\n    Transform DICOM images to NIfTI format and save to dataset directory. Delete DICOM files.\n    :param dataset:  Dataset (train, test).\n    :param kwargs: Dictionary of uniform parameters.\n    :return: None\n    \"\"\"\n\n    # Iterate over subjects and series\n    start = datetime.datetime.now()\n    counter = 0\n    count = len(pd.read_csv(kwargs['datasets'][dataset]['meta']))\n    dataset = kwargs['datasets'][dataset]['path']\n    subjects = os.listdir(dataset)\n    for subject in subjects:\n        images = os.listdir(os.path.join(dataset, subject))\n        for image in images:\n\n            # Assign values to Modality and SeriesNumber attributes\n            files = os.listdir(os.path.join(dataset, subject, image))\n            for file in files:\n                path = os.path.join(dataset, subject, image, file)\n                ct = pydicom.dcmread(path)\n                ct.Modality = 'CT'\n                ct.SeriesNumber = image\n                ct.save_as(path)\n\n            # Convert to NIfTI format and delete image directory\n            path = os.path.join(dataset, subject, image)\n            dicom2nifti.convert_dir.convert_directory(dicom_directory=path, output_folder=dataset)\n            shutil.rmtree(path)\n\n            # Print progress report\n            counter += 1\n            progress = '[' + '.' * int(counter / 100) + ' ' * int((count - counter) / 100) + ']'\n            sys.stdout.write('\\rProcessed {}/{} images: {} elapsed time: {}'.format(\n                counter, count, progress, datetime.datetime.now() - start))\n            sys.stdout.flush()\n\n        # Delete subject directory\n        shutil.rmtree(os.path.join(dataset, subject))\n\n    return None","metadata":{},"execution_count":null,"outputs":[]}]}