{"cells":[{"metadata":{"papermill":{"duration":0.007492,"end_time":"2020-09-15T11:13:40.187079","exception":false,"start_time":"2020-09-15T11:13:40.179587","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# OSIC Pulmonary Fibrosis Progression\n## Predicting patients' lung function decline based on chest CT scans\n\nAuthor: Hunter S. DiCicco (`hsd`, `igokaskowitz`)\n\nIn this notebook I'll be examining the OSIC P. Fibrosis dataset and hopefully moving on to some EDA soon.\n\n### TODO:\n\n* [DONE] Get an idea of how many scans each patient has available\n* [DONE] Visualize some sequential images from random patients\n* Verify that ALL scans examine the chest, subset if not\n* Batch-convert all of the scans into images (?)\n* Decide on an approach for the image regression\n\n### Imports and Setup"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":1.407138,"end_time":"2020-09-15T11:15:12.400393","exception":false,"start_time":"2020-09-15T11:15:10.993255","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import os\nimport glob\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torchvision\n\nimport pydicom\n\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","papermill":{"duration":0.109874,"end_time":"2020-09-15T11:15:12.609921","exception":false,"start_time":"2020-09-15T11:15:12.500047","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"DEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint('Using device:', DEVICE)\nDATA_PATH = '/kaggle/input/osic-pulmonary-fibrosis-progression/'","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.119036,"end_time":"2020-09-15T11:15:12.836553","exception":false,"start_time":"2020-09-15T11:15:12.717517","status":"completed"},"tags":[]},"cell_type":"markdown","source":"### Examining some of the training data\n\nHow many scans does each patient have?"},{"metadata":{"papermill":{"duration":3.845818,"end_time":"2020-09-15T11:15:16.794439","exception":false,"start_time":"2020-09-15T11:15:12.948621","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"scan_nums_per_patient = {os.path.split(dirname)[-1]:len(filenames)\n                         for dirname, _, filenames\n                         in os.walk(DATA_PATH + 'train/')}\nscan_nums_per_patient = pd.Series(tuple(scan_nums_per_patient.values()),\n                                  index=scan_nums_per_patient.keys())\n_, ax = plt.subplots(figsize=(10, 8))\nscan_nums_per_patient.plot(kind='hist', bins=30, ax=ax);","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.102791,"end_time":"2020-09-15T11:15:17.001332","exception":false,"start_time":"2020-09-15T11:15:16.898541","status":"completed"},"tags":[]},"cell_type":"markdown","source":"This shows that it might be a good idea to filter for patients with greater than 100 scans or so. Intuition says that we want the most granular chest scans we can get.\n\nJust to get a feel for things, let's pick a patient at random and look at a random one of their DICOM scans:"},{"metadata":{"papermill":{"duration":0.16574,"end_time":"2020-09-15T11:15:17.271116","exception":false,"start_time":"2020-09-15T11:15:17.105376","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# each patient has their own directory containing their\n# entire history of DICOM scan sets\n\n# list of all patients available in the training set\npatients = glob.glob(DATA_PATH + 'train/*')\n\n# pick a patient at random from the above list\n# with >0 shots\nsample_patient = np.random.choice(patients[1:])\n\n# now we need a list of all the paths to their scans\npatient_scans = glob.glob(sample_patient + '/*')\n\n# pick one of this patient's scans at random\nsample_scan = np.random.choice(patient_scans)\nprint(f\"Patient: {os.path.split(sample_patient)[-1]}\")\nprint(f\"Scan: {os.path.split(sample_scan)[-1]}\")\n\n# `pydicom.dcmread` will read a `.dcm` file\n# using the correct image encoding\n# it returns a `DataSet` iterable containing\n# one image slice and all of the patient's information\ndataset = pydicom.dcmread(sample_scan)\n\n# the `dir` keyword lets us explore the attributes,\n# class methods and instance methods of a particular object:\n#dir(dataset)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.114507,"end_time":"2020-09-15T11:15:17.495579","exception":false,"start_time":"2020-09-15T11:15:17.381072","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"type(dataset)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These files are *not just* images from CT scans. They are in a special format called DICOM, which is designed to immutably link a patient's medical information to the results of diagnostic procedures like CT scans. That said, our data contain lots more technical information about the scan itself that we wouldn't have seen from the surface.\n\nHere is a link to the [DICOM Standard Browser](https://dicom.innolitics.com/ciods). We're interested in the tags under the CT Image section. That page helped me understand a lot of what I'm showing here.\n\n`pydicom` has a robust way of accessing the DICOM taxonomy of patient data. Its `FileDataset` objects emulate Python dictionaries, and so can be reassembled as such using a dictionary comprehension.\n\nFrom there we could hand that dictionary over to the `pd.DataFrame` constructor, and eventually we'd have a big `DataFrame` with relevant technical information about the CT scans for each patient.\n\nBelow we'll show the string representation of the top level `DataElement`s in the `FileDataset` we chose:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And now we can go forth and reconstruct it as a dictionary that maps the DICOM keyword to its value.\nWe omit `PixelData` from the dictionary comprehension, because we already have a way to access the pixel array for every scan."},{"metadata":{"trusted":true},"cell_type":"code","source":"{data_elem.keyword: data_elem.value for data_elem in dataset.values() if not data_elem.keyword=='PixelData'}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So far I'm interested in the following features:\n* `Manufacturer` and `ManufacturerModelName`\n    * Do different CT scan machines have 'quirks?' It makes sense to think that these machines are engineered to give reproduceable scans, but there still may be useful differences\n* `ImagePositionPatient` and `ImageOrientationPatient`\n    * My hunch says that `PatientPosition` to `ImageOrientationPatient` is one to one, but the image position relative to the frame of reference is good for knowing generally what area of the chest we're looking at.\n* `SliceLocation`\n    * Same as above, I think that this is an indexed version of the `ImagePositionPatient` metric.\n* `KVP` (stands for Kilovoltage peak, the highest voltage applied to the X-ray tube during scanning) and `XRayTubeCurrent`\n    * Does electrical power have an effect on other qualities of a scan? Thinking the other way around, perhaps it would also be helpful to know if something anomalous happened to the power supply if we see something odd in a particular scan.\n\nFinally, let's show the image associated with our `FileDataset`!"},{"metadata":{"papermill":{"duration":0.284769,"end_time":"2020-09-15T11:15:17.882386","exception":false,"start_time":"2020-09-15T11:15:17.597617","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"plt.imshow(dataset.pixel_array, cmap=plt.cm.bone);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Very cool. I'd also like to get a succinct view of many slices of this patient's CT scan. "},{"metadata":{"papermill":{"duration":0.368323,"end_time":"2020-09-15T11:15:18.365675","exception":false,"start_time":"2020-09-15T11:15:17.997352","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# sort our sample patient's scan by the number in the filename\npatient_scans_by_number = sorted(patient_scans,\n                                 key=lambda x: \n                                 int(os.path.splitext(os.path.split(x)[-1])[0]))\n\n# choose nine equally-spaced sequential images from this patient's set of scans\nidx = np.round(np.linspace(0, len(patient_scans_by_number) - 1, 9)).astype(int)\n# subset the scan images with the indices we made above\nsample_scan2 = [patient_scans_by_number[i] for i in idx]\n\n# read each of the nine images with `pydicom`\nsample_sets = [pydicom.dcmread(image) for image in sample_scan2]\n# extract the pixel array from each of the nine images\nimages = [image.pixel_array for image in sample_sets]\n\n# tricky thing here:\n# so far i've noticed that most of these images are encoded as floats\n# but there are some edge cases encoded as `uint16`s, which are too big\n# for `torch` to construct a `torch.Tensor` out of.\n# if we happen upon one of those cases, we need to convert the image\n# array to `uint8` so that `torch` doesn't throw a `TypeError`\n\n# we use `torchvision.transforms.ToTensor` because if that edge case arises,\n# it will handle the logic required to scale a [0,255] RGB image into \n# the [0,1] RGB format desired by `plt`\ntry:\n    images_tensor = [torchvision.transforms.ToTensor()(image) for image in images]\nexcept TypeError:\n    images_tensor = [torchvision.transforms.ToTensor()(image.astype('uint8')) for image in images]\n\n# use `make_grid` to arrange a 3x3 visualization of the entire chest scan\n# add some padding so that the grid isn't so crowded\ngrid = torchvision.utils.make_grid(images_tensor, nrow=3, padding=100)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.681224,"end_time":"2020-09-15T11:15:19.155847","exception":false,"start_time":"2020-09-15T11:15:18.474623","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# create some axes with a nice, big `figsize` so we can look closely\n_, ax = plt.subplots(figsize=(12,10))\n\n# permute the grid tensor so it has the image number as the last channel\n# finally, show the image grid with a pretty CT scan color palette!\nax.imshow(grid.permute((1,2,0)), cmap=plt.cm.bone);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Awesome. While I was assembling that visualization I had some more questions about the validity of the data we're given, so I wrote some sanity check queries:"},{"metadata":{"papermill":{"duration":0.121609,"end_time":"2020-09-15T11:15:19.390868","exception":false,"start_time":"2020-09-15T11:15:19.269259","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# verifying that all of these are indeed part of a chest scan\nprint('All samples are chest scans:',\n      all([sample['BodyPartExamined'].value == 'Chest' for sample in sample_sets]))","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.121609,"end_time":"2020-09-15T11:15:19.390868","exception":false,"start_time":"2020-09-15T11:15:19.269259","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# verifying that all of these are indeed part of a chest scan\nprint('All samples used the LUNG kernel for reconstruction:',\n      all([sample['ConvolutionKernel'].value in ('LUNG', 'L') for sample in sample_sets]))","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.112386,"end_time":"2020-09-15T11:15:19.619278","exception":false,"start_time":"2020-09-15T11:15:19.506892","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}