{"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":"This notebook explores the DICOM metadata of the provided dataset.","metadata":{}},{"cell_type":"code","source":"!pip install dosma -q","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sn\nimport os\n\nfrom dataclasses import dataclass\nfrom pathlib import Path\nfrom tqdm.notebook import tqdm\n\nimport pydicom\nimport dosma\nfrom tqdm.notebook import tqdm\nfrom skimage.util import montage\n\nimport torch\nimport torch.nn as nn\nimport torchvision.transforms as TF\n\nfrom typing import Callable, Optional, Sequence, Tuple, Union\n\ntqdm.pandas()\n\n# /kaggle/input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.17625/12.dcm\n# (/kaggle/input/rsna-2022-cervical-spine-fracture-detection)/(train_images)/(1.2.826.0.1.3680043.17625)/(12.dcm)\n# root/<train/test>/<case>/<slice>\n# root/train.csv\n# root/test.csv\n# root/train_bounding_boxes.csv\n# root/sample_submission.csv\n\n# specific weights for loss\n# https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/340392\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-01T07:00:40.823317Z","iopub.execute_input":"2022-08-01T07:00:40.823856Z","iopub.status.idle":"2022-08-01T07:00:43.150252Z","shell.execute_reply.started":"2022-08-01T07:00:40.823814Z","shell.execute_reply":"2022-08-01T07:00:43.148967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define paths\n@dataclass\nclass folders:\n    root   = Path(\"/kaggle/input/rsna-2022-cervical-spine-fracture-detection\")\n    train  = Path(\"/kaggle/input/rsna-2022-cervical-spine-fracture-detection/train_images\")\n    test   = Path(\"/kaggle/input/rsna-2022-cervical-spine-fracture-detection/test_images\")\n    output = Path(\"/kaggle/working/output\")\n    temp   = Path(\"/kaggle/working/temp\")\n\n# CT windows\n@dataclass\nclass windows:\n    bone = (400,1000)\n    lung = (-600, 800)\n    soft = (100,200)\n\n# Orientations\n@dataclass\nclass orientation:\n    axial = (\"SI\",\"AP\",\"RL\")\n    sagittal = (\"RL\",\"SI\",\"AP\")\n    coronal = (\"AP\",\"SI\",\"RL\")","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:22:33.783133Z","iopub.execute_input":"2022-08-01T07:22:33.78364Z","iopub.status.idle":"2022-08-01T07:22:33.793815Z","shell.execute_reply.started":"2022-08-01T07:22:33.7836Z","shell.execute_reply":"2022-08-01T07:22:33.792581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for p in [folders.output, folders.temp]:\n    p.mkdir(0o755, exist_ok=True)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:56:37.165942Z","iopub.execute_input":"2022-08-01T06:56:37.166346Z","iopub.status.idle":"2022-08-01T06:56:37.18153Z","shell.execute_reply.started":"2022-08-01T06:56:37.166312Z","shell.execute_reply":"2022-08-01T06:56:37.180467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:07:22.795773Z","iopub.execute_input":"2022-08-01T07:07:22.796226Z","iopub.status.idle":"2022-08-01T07:07:29.068026Z","shell.execute_reply.started":"2022-08-01T07:07:22.796185Z","shell.execute_reply":"2022-08-01T07:07:29.066087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Build dicom index, should take ~ 2 min, from a previous run we know there are 711601 dcm-files\ncached = folders.output/\"full_dcm_index.csv\"\nif(os.path.exists(cached)):\n    dcm_index = pd.read_csv(cached, low_memory=False)\n    print(\"Loaded dcm_index from file.\")\nelse:\n    dcm_index = pd.DataFrame(tqdm(folders.train.glob(\"*/*.dcm\"), total=711601))\n    dcm_index.to_csv(cached, index=False)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:10:45.077239Z","iopub.execute_input":"2022-08-01T07:10:45.077672Z","iopub.status.idle":"2022-08-01T07:10:46.413199Z","shell.execute_reply.started":"2022-08-01T07:10:45.07764Z","shell.execute_reply":"2022-08-01T07:10:46.411819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Identify individual studies\nif(\"series\" not in dcm_index.columns):\n    dcm_index[\"series\"] = dcm_index[0].progress_apply(lambda x: x.parent.name)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:10:32.730901Z","iopub.execute_input":"2022-08-01T07:10:32.731749Z","iopub.status.idle":"2022-08-01T07:10:32.73912Z","shell.execute_reply.started":"2022-08-01T07:10:32.731698Z","shell.execute_reply":"2022-08-01T07:10:32.737679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# how many slices are there per study?\ndist = dcm_index.groupby([\"series\"]).count()\ndist.columns = [\"Number of Scans\"]\nprint(\"{} studies with on avg. {} slices (sd: {}, range: {} - {})\".format(*dist.describe().loc[[\"count\",\"mean\",\"std\",\"min\",\"max\"]].T.apply(lambda x: int(np.round(x))).T))\nprint(f\"95% of the studies contain between {np.percentile(dist, 2.5):0.0f} and {np.percentile(dist, 97.5):0.0f} slices.\")\np = sn.displot(dist, legend=False)\np.set(xlabel = \"Slices per study\", ylabel = \"Studies\", title = \"Histogram\")","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:58:56.06121Z","iopub.execute_input":"2022-08-01T06:58:56.06168Z","iopub.status.idle":"2022-08-01T06:58:56.75208Z","shell.execute_reply.started":"2022-08-01T06:58:56.061642Z","shell.execute_reply":"2022-08-01T06:58:56.751197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# select the elements we are interested in \nof_interest = ['BitsAllocated', 'BitsStored', 'Columns', 'ContentDate', 'ContentTime', 'HighBit', \n               'ImageOrientationPatient', 'ImagePositionPatient', 'InstanceNumber',  'PatientID',\n               'PatientName',  'PhotometricInterpretation',  'PixelRepresentation',  'PixelSpacing',\n               'RescaleIntercept',  'RescaleSlope', 'Rows',  'SOPInstanceUID',  'SamplesPerPixel',\n               'SeriesInstanceUID',  'SliceThickness',  'StudyInstanceUID',  'WindowCenter', 'WindowWidth']","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:58:58.858168Z","iopub.execute_input":"2022-08-01T06:58:58.858922Z","iopub.status.idle":"2022-08-01T06:58:58.86493Z","shell.execute_reply.started":"2022-08-01T06:58:58.858883Z","shell.execute_reply":"2022-08-01T06:58:58.863945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# collect metadata from the first slice of every study ETA ~25s\nfirst_slices = dcm_index.groupby(\"series\").first()\n\ndcm_metadata = []\nfor i,x in tqdm(first_slices.iterrows(), total=len(first_slices)):\n    current = pydicom.read_file(x[0], stop_before_pixels=True)\n    dcm_metadata.append([x[0]] + [current.get(x) for x in of_interest])\n\ndcm_metadata = pd.DataFrame(dcm_metadata, columns=([\"path\"]+of_interest))","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:59:00.409963Z","iopub.execute_input":"2022-08-01T06:59:00.410776Z","iopub.status.idle":"2022-08-01T06:59:27.004454Z","shell.execute_reply.started":"2022-08-01T06:59:00.410719Z","shell.execute_reply":"2022-08-01T06:59:27.003188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# what columns give a signal\n# look how many unique values are in each column, and what their distribution is\nprint(f'{\"column\": <32} {\"n_unique\": <5}  value(s)')\nfor c in dcm_metadata.columns:\n    try:\n        n_unique = dcm_metadata[c].nunique()\n        if(n_unique == 1):\n            print(f\"{c: <25} {n_unique: >15}  {dcm_metadata[c][0]}\")\n        elif(1 < n_unique < 10):\n            print(f\"{c: <25} {n_unique: >15} \", dcm_metadata[c].value_counts().to_dict())\n        else:\n            print(f\"{c: <25} {n_unique: >15} \", \"many uniques\")\n    except Exception as e:\n        print(f\"{c: <25} {str(e): >46}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:59:27.930675Z","iopub.execute_input":"2022-08-01T06:59:27.931148Z","iopub.status.idle":"2022-08-01T06:59:27.985668Z","shell.execute_reply.started":"2022-08-01T06:59:27.931111Z","shell.execute_reply":"2022-08-01T06:59:27.984378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### The columns BitStored, Columns, HighBit, PixelRepresentation have a few outliers\n\n# BitStored: Number of bits stored for each pixel sample. Each sample shall have the same number of bits stored. See PS3.5 for further explanation.\n# Source: https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280101\n\n# HighBit: Most significant bit for pixel sample data. Each sample shall have the same high bit. High Bit (0028,0102) shall be one less than Bits Stored (0028,0101). See PS3.5 for further explanation.\n# Source: https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280102\n\n# PixelRepresentation:  Data representation of the pixel samples. Each sample shall have the same pixel representation.\n# Enumerated Values: 0000H unsigned integer. 0001H 2's complement\n# Source: https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280103\n\n### Rescale Intercept and Slope\n# different values, so that a scan-specific rescaling is necessary to obtain HU values\n\n### Slice thickness, many different values\n# plot histogram\n\n### Rows and Cols\n# few outliers, check if there are sag/cor reconstructed CTs in the data\n\n# BitsAllocated, ContentDate, SamplesPerPixel, PhotometricInterpretation don't need to be stored as they are the same for every row\n# path, PatientID, PatientName, SOPInstanceUID, SeriesInstanceUID, StudyInstanceUID are identifiers, no patient seems to have been included twice.\n\n### => Further investigation of: ContentTime, ","metadata":{"execution":{"iopub.status.busy":"2022-07-30T07:23:25.640661Z","iopub.execute_input":"2022-07-30T07:23:25.642588Z","iopub.status.idle":"2022-07-30T07:23:25.650536Z","shell.execute_reply.started":"2022-07-30T07:23:25.642515Z","shell.execute_reply":"2022-07-30T07:23:25.648652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Slice Thickness\nsn.displot(dcm_metadata[\"SliceThickness\"])\ndisplay(dcm_metadata[\"SliceThickness\"].describe())","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:59:32.555009Z","iopub.execute_input":"2022-08-01T06:59:32.555442Z","iopub.status.idle":"2022-08-01T06:59:32.806139Z","shell.execute_reply.started":"2022-08-01T06:59:32.555408Z","shell.execute_reply":"2022-08-01T06:59:32.804775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Rows and Cols\nmatrix_outliers = dcm_metadata[(dcm_metadata[\"Rows\"] != 512) | (dcm_metadata[\"Columns\"] != 512)]\ndisplay(matrix_outliers)\n\n# TODO: Look at CTs of outliers 860, 908, 1951","metadata":{"execution":{"iopub.status.busy":"2022-08-01T06:59:35.365033Z","iopub.execute_input":"2022-08-01T06:59:35.3657Z","iopub.status.idle":"2022-08-01T06:59:35.405614Z","shell.execute_reply.started":"2022-08-01T06:59:35.365664Z","shell.execute_reply":"2022-08-01T06:59:35.404408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load and view CT scans","metadata":{}},{"cell_type":"code","source":"# CT ops\ndef rescale_to_HU(scan: Union[torch.Tensor, np.ndarray],\n                  intercept: Union[int,float], \n                  slope: float\n                 ) -> Union[torch.Tensor,np.ndarray]:\n    \"\"\"Rescale tensor or ndarray to HU\"\"\"\n    return scan*slope + intercept\n\ndef read_scan(path: Path, \n              reformat: Optional[Tuple[str]] = ('SI', 'AP', 'RL'),\n              rescale_hu: bool = True,\n              transform: Optional[Union[nn.Module, Callable]] = None\n             ) -> dosma.MedicalVolume:\n    \"\"\"Read DICOM file from path. Returns dosma.MedicalVolume\"\"\"\n\n    if(os.path.exists(path) and os.path.isdir(path)):\n        # read\n        try:\n            scan = dosma.read(path, \"dicom\", group_by=\"PatientName\", unpack=True)\n        except Exception as e:\n            print(f\"Could not load {path.as_posix()}, Error: {str(e)}\")\n        \n        metadata = scan.headers()[0][0][0]\n        \n        # reformat\n        if(reformat is not None):\n            scan.reformat(reformat)\n        \n        if(transform is not None):\n            scan = transform(scan)\n\n        # rescale\n        if(rescale_hu):\n            scan = rescale_to_HU(scan, \n                                 intercept=metadata.RescaleIntercept,\n                                 slope=metadata.RescaleSlope)\n        return scan\n\n    else:\n        raise NotADirectoryError(f\"path should point to an existing directory, received {str(path)}\")\n\n\ndef window_ct(scan: Union[torch.Tensor, np.ndarray], \n              wlevel: int, \n              wwidth: int, \n              rescale_to: Optional[Tuple[Union[int,float],Union[int,float]]]\n              ) -> Union[torch.Tensor, np.ndarray]:\n    \"\"\"Window CT scan\"\"\"\n\n    new_min = wlevel - (wwidth//2)\n    new_max = wlevel + (wwidth//2)\n    \n    clip_fn = torch.clamp if isinstance(scan, torch.Tensor) else np.clip\n    scan = clip_fn(scan, new_min, new_max)\n\n    if(rescale_to is not None):\n        min_fn, max_fn = (torch.min, torch.max) if isinstance(scan, torch.Tensor) else (np.min,np.max)\n        scan = (scan - new_min) / (new_max - new_min)\n        scan *= max_fn(rescale_to) - min_fn(rescale_to)\n        scan += min_fn(rescale_to)\n\n    return scan","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:00:46.963842Z","iopub.execute_input":"2022-08-01T07:00:46.965235Z","iopub.status.idle":"2022-08-01T07:00:46.984109Z","shell.execute_reply.started":"2022-08-01T07:00:46.965187Z","shell.execute_reply":"2022-08-01T07:00:46.982626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(25,20))\ntest = read_scan(first_slices[0][0].parent)\nwindowed1 = window_ct(test.reformat(orientation.coronal).volume, *windows.bone, rescale_to=(0.,1.))\nwindowed2 = window_ct(test.reformat(orientation.sagittal).volume, *windows.bone, rescale_to=(0.,1.))\n\nplt.imshow(montage(np.vstack([windowed1[125:301:25,:,:], windowed2[125:301:25,:,:]]),\n                   fill=windowed1.min(),\n                   grid_shape=(4,4)),\n           cmap=\"bone\")\nplt.show();plt.clf()","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:29:49.112017Z","iopub.execute_input":"2022-08-01T07:29:49.112437Z","iopub.status.idle":"2022-08-01T07:29:54.376837Z","shell.execute_reply.started":"2022-08-01T07:29:49.112406Z","shell.execute_reply":"2022-08-01T07:29:54.375521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at matrix size outliers\nfor i,x in matrix_outliers.iterrows():\n    test = read_scan(x.path.parent)\n    print(test)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T07:05:05.392379Z","iopub.execute_input":"2022-08-01T07:05:05.392831Z","iopub.status.idle":"2022-08-01T07:05:35.397446Z","shell.execute_reply.started":"2022-08-01T07:05:05.392795Z","shell.execute_reply":"2022-08-01T07:05:35.396558Z"},"trusted":true},"execution_count":null,"outputs":[]}]}