{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Reconstructing 3D volumes from metadata\nIn this kernel we will use the `StudyInstanceUID` in the metadata to group together images from the same scan. We will then sort the images using `ImagePositionPatient_2` and create 3D volumes. The 2D DICOM images supplied to us are *axial* slices. Once we have a 3D volume we will be able to make *sagittal* and *coronal* slices\n\n![](https://www.radiologycafe.com/images/basics/ct-planes.png)\n\nCredit due to this discussion https://www.kaggle.com/c/rsna-intracranial-hemorrhage-detection/discussion/109953#latest-639253"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport pydicom\nimport matplotlib.pylab as plt\nfrom tqdm import tqdm_notebook\n%matplotlib inline\n\ndata_path = \"../input/rsna-intracranial-hemorrhage-detection\"\nmetadata_path = \"../input/rsna-ich-metadata\"\nos.listdir(metadata_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prepare the labels & metadata\nThe metadata was extracted beforehand using `pydicom`. This takes a while so I saved the results in these parquet files so they don't need to be generated each time."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(f'{data_path}/stage_1_train.csv').drop_duplicates()\ntrain_df['ImageID'] = train_df['ID'].str.slice(stop=12)\ntrain_df['Diagnosis'] = train_df['ID'].str.slice(start=13)\ntrain_labels = train_df.pivot(index=\"ImageID\", columns=\"Diagnosis\", values=\"Label\")\ntrain_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_metadata(image_dir):\n\n    labels = [\n        'BitsAllocated', 'BitsStored', 'Columns', 'HighBit', \n        'ImageOrientationPatient_0', 'ImageOrientationPatient_1', 'ImageOrientationPatient_2',\n        'ImageOrientationPatient_3', 'ImageOrientationPatient_4', 'ImageOrientationPatient_5',\n        'ImagePositionPatient_0', 'ImagePositionPatient_1', 'ImagePositionPatient_2',\n        'Modality', 'PatientID', 'PhotometricInterpretation', 'PixelRepresentation',\n        'PixelSpacing_0', 'PixelSpacing_1', 'RescaleIntercept', 'RescaleSlope', 'Rows', 'SOPInstanceUID',\n        'SamplesPerPixel', 'SeriesInstanceUID', 'StudyID', 'StudyInstanceUID', \n        'WindowCenter', 'WindowWidth', 'Image',\n    ]\n\n    data = {l: [] for l in labels}\n\n    for image in tqdm_notebook(os.listdir(image_dir)):\n        data[\"Image\"].append(image[:-4])\n\n        ds = pydicom.dcmread(os.path.join(image_dir, image))\n\n        for metadata in ds.dir():\n            if metadata != \"PixelData\":\n                metadata_values = getattr(ds, metadata)\n                if type(metadata_values) == pydicom.multival.MultiValue and metadata not in [\"WindowCenter\", \"WindowWidth\"]:\n                    for i, v in enumerate(metadata_values):\n                        data[f\"{metadata}_{i}\"].append(v)\n                else:\n                    if type(metadata_values) == pydicom.multival.MultiValue and metadata in [\"WindowCenter\", \"WindowWidth\"]:\n                        data[metadata].append(metadata_values[0])\n                    else:\n                        data[metadata].append(metadata_values)\n\n    return pd.DataFrame(data).set_index(\"Image\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Generate metadata dataframes\ntrain_metadata = get_metadata(os.path.join(data_path, \"stage_1_train_images\"))\ntest_metadata = get_metadata(os.path.join(data_path, \"stage_1_test_images\"))\n\ntrain_metadata.to_parquet(f'{data_path}/train_metadata.parquet.gzip', compression='gzip')\ntest_metadata.to_parquet(f'{data_path}/test_metadata.parquet.gzip', compression='gzip')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_metadata = pd.read_parquet(f'{metadata_path}/train_metadata.parquet.gzip')\ntest_metadata = pd.read_parquet(f'{metadata_path}/test_metadata.parquet.gzip')\n\ntrain_metadata[\"Dataset\"] = \"train\"\ntest_metadata[\"Dataset\"] = \"test\"\n\ntrain_metadata = train_metadata.join(train_labels)\n\nmetadata = pd.concat([train_metadata, test_metadata], sort=True)\nmetadata.sort_values(by=\"ImagePositionPatient_2\", inplace=True, ascending=False)\nmetadata.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"metadata[\"StudyInstanceUID\"].nunique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Group the unique studies"},{"metadata":{"trusted":true},"cell_type":"code","source":"studies = metadata.groupby(\"StudyInstanceUID\")\nstudies_list = list(studies)\n\nstudy_name, study_df = studies_list[0]\nstudy_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"studies.size().describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(studies.size());","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems like the unique studies can have anywhere between 20 and 60 images (i.e. axial slices). Perhaps they need to be resized to a constant z-dimension? Or if x & y is constant, use an architecture that can cope with this e.g. adaptive pooling layers?\n\nLet's check if any studies straddle train/test:"},{"metadata":{"trusted":true},"cell_type":"code","source":"studies.filter(lambda x: x[\"Dataset\"].nunique() > 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks like none - which is good, however we still need to remember that patients can have multiple scans and these can be across the train & stage 1 test."},{"metadata":{},"cell_type":"markdown","source":"# Create a 3D volume for a single study\nWe'll use the first study in the grouped `studies` which is comprised of 40 individial axial DICOM images\n\nThanks to this notebook for the windowing functions https://www.kaggle.com/wfwiggins203/eda-dicom-tags-windowing-head-cts"},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_img(dcm, width=None, level=None, norm=True):\n    pixels = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    \n    # Pad non-square images\n    if pixels.shape[0] != pixels.shape[1]:\n        (a,b) = pixels.shape\n        if a > b:\n            padding = ((0, 0), ((a-b) // 2, (a-b) // 2))\n        else:\n            padding = (((b-a) // 2, (b-a) // 2), (0, 0))\n        pixels = np.pad(pixels, padding, mode='constant', constant_values=0)\n            \n    if not width:\n        width = dcm.WindowWidth\n        if type(width) != pydicom.valuerep.DSfloat:\n            width = width[0]\n    if not level:\n        level = dcm.WindowCenter\n        if type(level) != pydicom.valuerep.DSfloat:\n            level = level[0]\n    lower = level - (width / 2)\n    upper = level + (width / 2)\n    img = np.clip(pixels, lower, upper)\n    \n    if norm:\n        return (img - lower) / (upper - lower)\n    else:\n        return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"volume, labels = [], []\nfor index, row in study_df.iterrows():\n    if row[\"Dataset\"] == \"train\":\n        dcm = pydicom.dcmread(os.path.join(data_path, \"stage_1_train_images\", index+\".dcm\"))\n    else:\n        dcm = pydicom.dcmread(os.path.join(data_path, \"stage_1_test_images\", index+\".dcm\"))\n        \n    img = window_img(dcm)\n    label = row[[\"any\", \"epidural\", \"intraparenchymal\", \"intraventricular\", \"subarachnoid\", \"subdural\"]]\n    volume.append(img)\n    labels.append(label)\n    \nvolume = np.array(volume)\nlabels = np.array(labels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"volume.shape, labels.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The provided DICOM images are axial slices. Let's use our new 3D volume to create sagittal and coronal slices\n* Red line - axial plane\n* Green line - sagittal plane\n* Blue line - coronal plane"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Axial\nplt.figure(figsize=(8, 8))\nplt.imshow(volume[20, :, :], cmap=plt.cm.bone)\nplt.vlines(300, 0, 512, colors='g')\nplt.hlines(300, 0, 512, colors='b');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Sagittal\nplt.figure(figsize=(8, 8))\nplt.imshow(volume[:, :, 300], aspect=9, cmap=plt.cm.bone)\nplt.vlines(300, 0, 40, colors='b')\nplt.hlines(20, 0, 512, colors='r');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Coronal\nplt.figure(figsize=(8, 8))\nplt.imshow(volume[:, 300, :], aspect=9, cmap=plt.cm.bone)\nplt.vlines(300, 0, 40, colors='g')\nplt.hlines(20, 0, 512, colors='r');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.1"}},"nbformat":4,"nbformat_minor":1}