{"cells":[{"metadata":{},"cell_type":"markdown","source":"# [2020-PE] Preprocessing train data\nThis kernel was used to preprocess train data for the RSNA STR Pulmonary Embolism Detection competition. The code is largely from this [preprocessing and EDA kernel](https://www.kaggle.com/allunia/pulmonary-dicom-preprocessing) by [@allunia](https://www.kaggle.com/allunia) where the steps are explained in detail.\n\nAfter isolating the lung parts of each CT scan, all images are resized to 128x128 px and each exam is downsampled to 20 images. The final data is therefore one data cube of size 20x128x128 for each exam and can be found [here](https://www.kaggle.com/spacelx/2020pe-preprocessed-train-data) in npy format. Each exam is further accompanied by a ordered list of the original image names which were combined into the associated data cube.\n\nLastly, [this kernel](https://www.kaggle.com/spacelx/2020-pe-preprocessing-train-table) was used to compute a modified datatable for the train set to fit the format of the preprocessed data.\n\nhttps://www.kaggle.com/spacelx/2020pe-preprocessed-train-data\n\nhttps://www.kaggle.com/spacelx/2020-pe-preprocessing-train-table"},{"metadata":{"papermill":{"duration":0.005603,"end_time":"2020-09-19T15:12:29.831483","exception":false,"start_time":"2020-09-19T15:12:29.82588","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# Setup"},{"metadata":{"execution":{"iopub.execute_input":"2020-09-19T15:12:29.84615Z","iopub.status.busy":"2020-09-19T15:12:29.845482Z","iopub.status.idle":"2020-09-19T15:12:31.715847Z","shell.execute_reply":"2020-09-19T15:12:31.714727Z"},"papermill":{"duration":1.880109,"end_time":"2020-09-19T15:12:31.716005","exception":false,"start_time":"2020-09-19T15:12:29.835896","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\n# path management\nfrom pathlib import Path\n\n# DICOM file handling\nimport pydicom\n\n# image processing\nimport cv2\nfrom skimage import measure \nfrom skimage.morphology import disk, opening, closing\n\n# progress bars\nfrom tqdm import tqdm\n\n# zip files\nimport zipfile","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.execute_input":"2020-09-19T15:12:31.732601Z","iopub.status.busy":"2020-09-19T15:12:31.731318Z","iopub.status.idle":"2020-09-19T15:12:31.736686Z","shell.execute_reply":"2020-09-19T15:12:31.735951Z"},"papermill":{"duration":0.01553,"end_time":"2020-09-19T15:12:31.73682","exception":false,"start_time":"2020-09-19T15:12:31.72129","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"basepath = Path('../input/rsna-str-pulmonary-embolism-detection/')\nfor p in basepath.iterdir():\n    print(str(p))","execution_count":null,"outputs":[]},{"metadata":{"execution":{"iopub.execute_input":"2020-09-19T15:12:31.754837Z","iopub.status.busy":"2020-09-19T15:12:31.754144Z","iopub.status.idle":"2020-09-19T15:12:35.733772Z","shell.execute_reply":"2020-09-19T15:12:35.733138Z"},"papermill":{"duration":3.991822,"end_time":"2020-09-19T15:12:35.73389","exception":false,"start_time":"2020-09-19T15:12:31.742068","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath / 'train.csv')\ntrain['dcmpath'] = str(basepath) + '/train' + '/' + train.StudyInstanceUID + '/' + train.SeriesInstanceUID","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.004595,"end_time":"2020-09-19T15:12:35.743584","exception":false,"start_time":"2020-09-19T15:12:35.738989","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# Helper functions for preprocessing"},{"metadata":{"execution":{"iopub.execute_input":"2020-09-19T15:12:35.775787Z","iopub.status.busy":"2020-09-19T15:12:35.775155Z","iopub.status.idle":"2020-09-19T15:12:35.77821Z","shell.execute_reply":"2020-09-19T15:12:35.777659Z"},"papermill":{"duration":0.029841,"end_time":"2020-09-19T15:12:35.778331","exception":false,"start_time":"2020-09-19T15:12:35.74849","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"def open_files(dcmpath):\n    '''loads all scans in the given path and orders them'''\n    filelist = np.array(list(Path(dcmpath).iterdir()))\n    scans = np.array([pydicom.dcmread(str(file)) for file in filelist])\n    sortkey = np.argsort([float(x.ImagePositionPatient[2]) for x in scans])\n    return scans[sortkey], filelist[sortkey]\n\n\ndef transform_to_hu(scans):\n    '''transform all scans from raw data into Hounsfield units'''\n    # stack scans\n    imglist = []\n    for file in scans:\n        try:\n            tmp = file.pixel_array\n        except:\n            tmp = np.zeros((512,512))\n        imglist.append(tmp)\n    images = np.stack(imglist)\n    images = images.astype(np.int16)\n    \n    # threshold between air (0) and default mask value (-2000) using -1000 [raw values]\n    images[images <= -1000] = 0\n    \n    # convert to HU\n    for nnn in range(len(scans)):        \n        intercept = scans[nnn].RescaleIntercept\n        slope = scans[nnn].RescaleSlope        \n        if slope != 1:\n            images[nnn] = slope * images[nnn].astype(np.float64)\n            images[nnn] = images[nnn].astype(np.int16)            \n        images[nnn] += np.int16(intercept)    \n    return np.array(images, dtype=np.int16)\n\n\ndef segment_lung_mask(scans):\n    '''mask image and retain only lung sections'''\n    segmented = np.zeros(scans.shape)\n    for nnn in range(scans.shape[0]):\n\n        # segment into water-like (2) and air-like (1) parts\n        slice_binary = np.array((scans[nnn]>-320), dtype=np.int8) + 1\n        slice_label = measure.label(slice_binary)\n        \n        bad_labels = np.unique([\n            slice_label[0,:],\n            slice_label[-1,:],\n            slice_label[:,0],\n            slice_label[:,-1]\n        ])\n        for bbb in bad_labels:\n            slice_binary[slice_label == bbb] = 2\n\n        # invert, air-like is now 1\n        slice_binary -= 1\n        slice_binary = 1 - slice_binary\n        \n        segmented[nnn] = slice_binary.copy() * scans[nnn]\n    return segmented\n\n\ndef resize_scans(scans, NSCANS, NPX):\n    '''resize collections of scans to a common size'''\n    resized_scans = np.zeros((NSCANS, NPX, NPX))\n\n    split = np.linspace(0, scans.shape[0], num=NSCANS+1).astype(int)\n    for sss in range(NSCANS):\n        scan_selection = np.mean(scans[split[sss]:split[sss+1]], axis=0)\n        resized_scans[sss] = cv2.resize(scan_selection, (NPX, NPX))\n    return resized_scans\n\n\ndef load_scans(dcmpath, NSCANS, NPX):\n    '''load all files of a scan through the full pipeline'''\n    scans, filelist = open_files(dcmpath)\n    hu_scans = transform_to_hu(scans)\n    segmented_scans = segment_lung_mask(hu_scans)\n    resized_scans = resize_scans(segmented_scans, NSCANS, NPX)\n    return resized_scans, filelist\n\n\ndef preproc_scans(dcmpathlist, NSCANS, NPX, outdir):\n    '''preprocess list of scans through the full pipeline and save result in zipped output file'''\n    ziparchive = zipfile.ZipFile(Path(outdir), 'w', zipfile.ZIP_DEFLATED)\n    for dcmpath in tqdm(dcmpathlist):\n        # load and preprocess\n        scans, filelist = load_scans(dcmpath, NSCANS, NPX)\n        # get identifiers\n        series = Path(dcmpath).name\n        study = Path(dcmpath).parent.name\n        # save processed data\n        filename = (study + '_' + series + '_data.npy')\n        np.save(filename, scans)\n        ziparchive.write(filename)\n        Path(filename).unlink()\n        # save ordered filelist\n        filename = (study + '_' + series + '_list.npy')\n        np.save(filename, filelist)\n        ziparchive.write(filename)\n        Path(filename).unlink()\n    ziparchive.close()","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.00506,"end_time":"2020-09-19T15:12:35.788904","exception":false,"start_time":"2020-09-19T15:12:35.783844","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# Preprocess\nJust processing the first five exams here as an example"},{"metadata":{"execution":{"iopub.execute_input":"2020-09-19T15:12:35.808061Z","iopub.status.busy":"2020-09-19T15:12:35.806943Z","iopub.status.idle":"2020-09-20T00:04:30.73818Z","shell.execute_reply":"2020-09-20T00:04:30.738607Z"},"papermill":{"duration":31914.944532,"end_time":"2020-09-20T00:04:30.738791","exception":false,"start_time":"2020-09-19T15:12:35.794259","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"scanlist = np.unique(train.dcmpath.values)\npreproc_scans(scanlist[:5], 20, 128, 'proc_20_128_train.zip')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Plot example"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nsample_scans, filelist = load_scans(scanlist[150], 20, 128)\n\nfig, ax = plt.subplots(5, 4, figsize=(20,20))\nax = ax.flatten()\nfor m in range(20):\n    ax[m].imshow(sample_scans[m], cmap='Blues_r')","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}