{"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":"## Introduction\n\nIn this notebook, we are going to use our previously trained CNN model (refer Baseline 1 - 2D CNN Image-level Predictions [Training](https://www.kaggle.com/code/pankajpansari/baseline-1-2d-cnn-image-level-prediction-train) notebook) to make prediction on unseen CT scans. We are also going to use this notebook to make our submission for the competition. We make use of the [fastai](https://docs.fast.ai) library.\n\nThe advantage of separating the training and inference parts of the pipeline is that we don't need to spend time and compute retraining our model in the cloud. We can simply import the saved model and make predictions. Moreover, this separation enables more flexible development. We can focus on improving training or speeding up inference in an indepenedent way.\n\nOur model makes predictions on individual images. There are sequence of images for each patient. We need to take the image-level predictions and aggregrate them to make prediction for each patient. Instead of taking mean over all images in a scan, we'll take 95% quantile.","metadata":{}},{"cell_type":"markdown","source":"## Code\n\n#### Import libaries\n\nWe're going to import the necessary libraries.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os, random\nfrom fastai.vision.all import *\nfrom fastai.medical.imaging import *\nimport shutil\nimport pydicom\nimport cv2\nimport glob\nimport time\nfrom rsna_2023_atd_metric import score\nimport tqdm\n\nfrom PIL import Image\n\nrandom.seed(1441)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:05.388887Z","iopub.execute_input":"2023-10-05T12:37:05.389779Z","iopub.status.idle":"2023-10-05T12:37:18.87781Z","shell.execute_reply.started":"2023-10-05T12:37:05.389742Z","shell.execute_reply":"2023-10-05T12:37:18.876814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### DICOM Preprocessing\n\nThe following resolves a particular issue in the dataset with regards to appearances of some scan images. Here, we also do image rescaling for better clarity.","metadata":{}},{"cell_type":"code","source":"def standardize_pixel_array(dcm):\n    \"\"\"\n    Source : https://www.kaggle.com/competitions/rsna-2023-abdominal-trauma-detection/discussion/427217\n    \"\"\"\n    # Correct DICOM pixel_array if PixelRepresentation == 1.\n    #dcm = pydicom.dcmread(fn)\n    pixel_array = dcm.pixel_array\n    if dcm.PixelRepresentation == 1:\n        bit_shift = dcm.BitsAllocated - dcm.BitsStored\n        dtype = pixel_array.dtype \n        pixel_array = (pixel_array << bit_shift).astype(dtype) >>  bit_shift\n#         pixel_array = pydicom.pixel_data_handlers.util.apply_modality_lut(new_array, dcm)\n\n    intercept = float(dcm.RescaleIntercept)\n    slope = float(dcm.RescaleSlope)\n    center = int(dcm.WindowCenter)\n    width = int(dcm.WindowWidth)\n    low = center - width / 2\n    high = center + width / 2    \n    \n    pixel_array = (pixel_array * slope) + intercept\n    pixel_array = np.clip(pixel_array, low, high)\n\n    return pixel_array","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:18.879552Z","iopub.execute_input":"2023-10-05T12:37:18.880753Z","iopub.status.idle":"2023-10-05T12:37:18.887353Z","shell.execute_reply.started":"2023-10-05T12:37:18.88072Z","shell.execute_reply":"2023-10-05T12:37:18.886437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Import Custom Code and Model\n\nFor us to import the saved model, fast.ai requires that we keep custom code such as custom model, custom loss, and custom metrics in the same place here. ","metadata":{}},{"cell_type":"code","source":"class MultiHeadModel(Module):\n    \n    def __init__(self, body):\n    \n        self.body = body\n        nf = num_features_model(nn.Sequential(*self.body.children()))\n\n        self.bowel = create_head(nf, 1)\n        self.extravasation = create_head(nf, 1)\n        self.kidney = create_head(nf, 3)\n        self.liver = create_head(nf, 3)\n        self.spleen = create_head(nf, 3)\n        \n    def forward(self, x):\n        \n        y = self.body(x)\n        bowel = self.bowel(y)\n        extravasation = self.extravasation(y)\n        kidney = self.kidney(y)\n        liver = self.liver(y)\n        spleen = self.spleen(y)\n        return [bowel, extravasation, kidney, liver, spleen]","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:18.889096Z","iopub.execute_input":"2023-10-05T12:37:18.889856Z","iopub.status.idle":"2023-10-05T12:37:18.900392Z","shell.execute_reply.started":"2023-10-05T12:37:18.889825Z","shell.execute_reply":"2023-10-05T12:37:18.899386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CombinationLoss(Module):\n    \"Cross entropy loss on multiple targets\"\n    def __init__(self, func = F.cross_entropy, weights = [2, 6, 3, 3, 3]):\n        self.func = func\n        self.w = weights\n        \n    def forward(self, xs, *ys, reduction = 'mean'):\n        loss = 0\n    \n        for i, w, x, y in zip(range(len(xs)), self.w, xs, ys):\n            if i < 2:\n                loss += w*F.binary_cross_entropy_with_logits(x, y.unsqueeze(1).float(), reduction = reduction)\n            else:\n                #import pdb;pdb.set_trace()\n                loss += w*F.cross_entropy(x, y, reduction = reduction)\n        return loss","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:18.903175Z","iopub.execute_input":"2023-10-05T12:37:18.903827Z","iopub.status.idle":"2023-10-05T12:37:18.912363Z","shell.execute_reply.started":"2023-10-05T12:37:18.903707Z","shell.execute_reply":"2023-10-05T12:37:18.911367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import recall_score\n\nclass RecallPartial(Metric):\n    \"Stores predictions and targets on CPU in accumulate to perform final calculations with `func`.\"\n    def __init__(self, a=0, **kwargs):\n        self.func = partial(recall_score, average='macro', zero_division=0)\n        self.a = a\n\n    def reset(self): self.targs,self.preds = [],[]\n\n    def accumulate(self, learn):\n        pred = learn.pred[self.a].argmax(-1)\n        targ = learn.y[self.a]\n        pred,targ = to_detach(pred),to_detach(targ)\n        pred,targ = flatten_check(pred,targ)\n        self.preds.append(pred)\n        self.targs.append(targ)\n\n    @property\n    def value(self):\n        if len(self.preds) == 0: return\n        preds,targs = torch.cat(self.preds),torch.cat(self.targs)\n        return self.func(targs, preds)\n\n    @property\n    def name(self): return 'recall_' + str(self.a+1)\n    \nclass RecallCombine(Metric):\n    def accumulate(self, learn):\n        scores = [learn.metrics[i].value for i in range(3)]\n        self.combine = np.average(scores, weights=[2,1,1])\n\n    @property\n    def value(self):\n        return self.combine","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:18.91356Z","iopub.execute_input":"2023-10-05T12:37:18.914507Z","iopub.status.idle":"2023-10-05T12:37:18.927442Z","shell.execute_reply.started":"2023-10-05T12:37:18.914477Z","shell.execute_reply":"2023-10-05T12:37:18.926482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We load the saved model. Note that we want our model to load on the GPU, because we're going to use GPU for inference. To ensure this, we need to set the _cpu_ option to __False__, else the model gets loaded on the CPU by default.","metadata":{}},{"cell_type":"code","source":"learn = load_learner('/kaggle/input/rsna2023-atd-2d-cnn-image-level-model-2/model_2.pt', cpu = False)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:18.928638Z","iopub.execute_input":"2023-10-05T12:37:18.929613Z","iopub.status.idle":"2023-10-05T12:37:25.915687Z","shell.execute_reply.started":"2023-10-05T12:37:18.929581Z","shell.execute_reply":"2023-10-05T12:37:25.914742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### DICOM to PNG\n\nThe original train/test CT scan images are in DICOM format, but we're going to use the PNG format for making predictions.\n\nWe're going to use a shortcut. Instead of making predictions on all test DICOM images, we are going to use every 10th (stride) image in each scan series. Using all images will take up a lot of time for inference, mainly due to conversion of DICOM to PNG.","metadata":{}},{"cell_type":"code","source":"TEST_PATH = '/kaggle/input/rsna-2023-abdominal-trauma-detection/test_images/'\nSAVE_FOLDER = 'temp_folder/'\nSIZE = 128\nSTRIDE = 10\n\nif not os.path.exists(SAVE_FOLDER):\n    os.makedirs(SAVE_FOLDER)\n\nprint('Number of test patients:', len(os.listdir(TEST_PATH)))","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:25.916922Z","iopub.execute_input":"2023-10-05T12:37:25.917536Z","iopub.status.idle":"2023-10-05T12:37:25.927356Z","shell.execute_reply.started":"2023-10-05T12:37:25.917499Z","shell.execute_reply":"2023-10-05T12:37:25.926314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following function processes the scan images for a given patients. For each series, it takes every STRIDE-th image, reads and resizes it, and saves it as PNG.","metadata":{}},{"cell_type":"code","source":"def convert_dicom_to_png(patient, size = 128):\n    \n    for study in (sorted(os.listdir(TEST_PATH + patient))):\n        imgs = {}\n        for f in sorted(glob.glob(TEST_PATH + f\"{patient}/{study}/*.dcm\"))[::STRIDE]:\n            \n            dicom = pydicom.dcmread(f)\n            pos_z = dicom[(0x20, 0x32)].value[-1]\n            img = standardize_pixel_array(dicom)\n            \n            img = (img - img.min())/(img.max() - img.min() + 1e-6)\n            imgs[pos_z] = img\n                \n        for i, k in enumerate(sorted(imgs.keys())):\n            \n            img = imgs[k]\n            \n            img = cv2.resize(img, (size, size))\n            cv2.imwrite(SAVE_FOLDER + f\"{patient}_{study}_{i}.png\", (img * 255).astype(np.uint8))\n        \n#_ = Parallel(n_jobs = 2)(\n#    delayed(convert_dicom_to_png)(patient, size=SIZE)\n#    for patient in tqdm(os.listdir(TEST_PATH))\n#    )\n    ","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:25.928727Z","iopub.execute_input":"2023-10-05T12:37:25.929698Z","iopub.status.idle":"2023-10-05T12:37:25.939377Z","shell.execute_reply.started":"2023-10-05T12:37:25.929665Z","shell.execute_reply":"2023-10-05T12:37:25.93831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def merge_arr(a, b):\n    return np.concatenate((a, b.numpy()), axis = 0)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:25.94056Z","iopub.execute_input":"2023-10-05T12:37:25.940902Z","iopub.status.idle":"2023-10-05T12:37:25.950911Z","shell.execute_reply.started":"2023-10-05T12:37:25.940871Z","shell.execute_reply":"2023-10-05T12:37:25.950045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Obtain Predictions\n\nIn the following, we loop over the patients. For each patient, we convert the DICOM images to rescaled PNG images and save them in a buffer folder. Then, we make a dataloader using these images and make predictions. We keep the batch size relatively high to make good use of GPU. We save the predictions in a list and then delete the images since we don't need them further.","metadata":{}},{"cell_type":"code","source":"patients = os.listdir(TEST_PATH)\n\nbowel_preds, extrav_preds = np.array([]).reshape(0), np.array([]).reshape(0)\nkidney_preds, liver_preds, spleen_preds = np.array([]).reshape(0, 3), np.array([]).reshape(0, 3), np.array([]).reshape(0, 3)\nfnames_list = []\n\nstart = time.time()\nsigm = torch.nn.Sigmoid()\nsoftm = torch.nn.Softmax(dim = 1)\n\nfor idx, patient in enumerate(patients):\n    \n    convert_dicom_to_png(patient, SIZE)\n    files = get_image_files(SAVE_FOLDER)\n    test_dl = learn.dls.test_dl(files, with_labels = False, device = 'cuda', bs = 128)\n\n    preds = learn.get_preds(dl = test_dl)[0]\n        \n    bowel_preds = merge_arr(bowel_preds, sigm(preds[0]).squeeze(-1))\n    extrav_preds = merge_arr(extrav_preds, sigm(preds[1]).squeeze(-1))\n    kidney_preds = merge_arr(kidney_preds, softm(preds[2]))\n    liver_preds = merge_arr(liver_preds, softm(preds[3]))\n    spleen_preds = merge_arr(spleen_preds, softm(preds[4]))\n       \n    fnames_list.append(files)\n    \n    for file in files:\n        os.remove(file)\n    \n    if (idx + 1) % 5 == 0:\n        end = time.time()\n        print(f'{idx + 1} patients processed.')\n        print(f'Time elapsed: {end - start} ')\n        print(f'Avg time per patient: {(end - start)/(idx + 1)}')","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:25.954415Z","iopub.execute_input":"2023-10-05T12:37:25.955109Z","iopub.status.idle":"2023-10-05T12:37:32.653932Z","shell.execute_reply.started":"2023-10-05T12:37:25.955077Z","shell.execute_reply":"2023-10-05T12:37:32.652891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Post-processing output\n\nWe're going to take the predictions made on each scan image, and aggregate them to derive predictions for each patient.","metadata":{}},{"cell_type":"code","source":"from itertools import chain\nfnames_list = list(chain.from_iterable(fnames_list))","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.655555Z","iopub.execute_input":"2023-10-05T12:37:32.655819Z","iopub.status.idle":"2023-10-05T12:37:32.662875Z","shell.execute_reply.started":"2023-10-05T12:37:32.655796Z","shell.execute_reply":"2023-10-05T12:37:32.661896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We assign the appropriate probabilities to various conditions using the model predictions.","metadata":{}},{"cell_type":"code","source":"test_files_probs = pd.DataFrame()\n\ntest_files_probs['fname'] = pd.Series(fnames_list, dtype = 'string')\n\ntest_files_probs['bowel_injury'] = pd.Series(bowel_preds)\ntest_files_probs['extravasation_injury'] = pd.Series(extrav_preds)\ntest_files_probs['kidney_low'] = pd.Series(kidney_preds[:, 1])\ntest_files_probs['kidney_high'] = pd.Series(kidney_preds[:, 2])\ntest_files_probs['liver_low'] = pd.Series(liver_preds[:, 1])\ntest_files_probs['liver_high'] = pd.Series(liver_preds[:, 2])\ntest_files_probs['spleen_low'] = pd.Series(spleen_preds[:, 1])\ntest_files_probs['spleen_high'] = pd.Series(spleen_preds[:, 2])\n\n#test_files_probs","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.664376Z","iopub.execute_input":"2023-10-05T12:37:32.66519Z","iopub.status.idle":"2023-10-05T12:37:32.689032Z","shell.execute_reply.started":"2023-10-05T12:37:32.665156Z","shell.execute_reply":"2023-10-05T12:37:32.688042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_files_probs.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.690397Z","iopub.execute_input":"2023-10-05T12:37:32.690763Z","iopub.status.idle":"2023-10-05T12:37:32.714493Z","shell.execute_reply.started":"2023-10-05T12:37:32.690714Z","shell.execute_reply":"2023-10-05T12:37:32.713508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We add a new column _'patient_id'_ derived from file name. This will help us group all the image predictions based on *patient_id*, and then we can use appropriate statistics for aggregation.","metadata":{}},{"cell_type":"code","source":"patient_id_list = []\nfor idx, fname in enumerate(test_files_probs['fname']):\n    patient_id_list.append(fname.split('/')[1].split('_')[0])\n    \ntest_files_probs['patient_id'] = pd.Series(patient_id_list, dtype = 'string')","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.715684Z","iopub.execute_input":"2023-10-05T12:37:32.716002Z","iopub.status.idle":"2023-10-05T12:37:32.722816Z","shell.execute_reply.started":"2023-10-05T12:37:32.715969Z","shell.execute_reply":"2023-10-05T12:37:32.721685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We first obtain all the unique *patient_ids*. Then for each patient, we derive 95% quantile statistics for disease values from all the entries which correspond to images for that patient.\n\nWe take the value for the healthy tags as 1 - sum of value of disease tags. We also make use of the scaling factor given in the evaluation metric to be able to improve score on the test set. We put a floor of eps so that values are not negative.\n\nNote that these values are not probabilities; they don't add up to 1. The scoring function normalizes them to add up to 1.","metadata":{}},{"cell_type":"code","source":"patients = set(test_files_probs.patient_id)\n\ncol_names = ['patient_id', 'bowel_healthy', 'bowel_injury',\n            'extravasation_healthy', 'extravasation_injury',\n            'kidney_healthy', 'kidney_low', 'kidney_high',\n            'liver_healthy', 'liver_low', 'liver_high',\n            'spleen_healthy', 'spleen_low', 'spleen_high']\n\npatient_probs = pd.DataFrame([], columns = col_names)\n    \n\nfor idx, pat in enumerate(patients):\n    p = test_files_probs[test_files_probs.patient_id == pat]\n    \n    bi = p.bowel_injury.quantile(q = 0.95)\n    ei = p.extravasation_injury.quantile(q = 0.95)\n    kl = p.kidney_low.quantile(q = 0.95)\n    kh = p.kidney_high.quantile(q = 0.95)\n    ll = p.liver_low.quantile(q = 0.95)\n    lh = p.liver_high.quantile(q = 0.95)\n    sl = p.spleen_low.quantile(q = 0.95)\n    sh = p.spleen_high.quantile(q = 0.95)\n    \n    eps = 1e-5\n    patient_probs.loc[idx] = [str(pat), max(eps, 1 - bi), 2*bi, max(eps, 1 - ei), 6*ei, max(eps, 1 - kl - kh), 2*kl, 4*kh, max(eps, 1 - ll - lh), 2*ll, 4*lh, max(eps, 1 - sl - sh), 2*sl, 4*sh]","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.724741Z","iopub.execute_input":"2023-10-05T12:37:32.725813Z","iopub.status.idle":"2023-10-05T12:37:32.7767Z","shell.execute_reply.started":"2023-10-05T12:37:32.725777Z","shell.execute_reply":"2023-10-05T12:37:32.775758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_probs.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.778206Z","iopub.execute_input":"2023-10-05T12:37:32.778533Z","iopub.status.idle":"2023-10-05T12:37:32.792595Z","shell.execute_reply.started":"2023-10-05T12:37:32.778504Z","shell.execute_reply":"2023-10-05T12:37:32.791645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We make sure that our submission is aligned with respect to the patient_ids.","metadata":{}},{"cell_type":"code","source":"sample_submission = pd.read_csv('/kaggle/input/rsna-2023-abdominal-trauma-detection/sample_submission.csv')\nsample_submission.patient_id = sample_submission.patient_id.astype(str)\n\n\npatient_probs = patient_probs.set_index('patient_id')\npatient_probs = patient_probs.reindex(index = sample_submission['patient_id'])\npatient_probs = patient_probs.reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.793943Z","iopub.execute_input":"2023-10-05T12:37:32.794729Z","iopub.status.idle":"2023-10-05T12:37:32.812819Z","shell.execute_reply.started":"2023-10-05T12:37:32.794693Z","shell.execute_reply":"2023-10-05T12:37:32.811833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we write the predictions for all patients to _submission.csv_ file and submit our notebook.","metadata":{}},{"cell_type":"code","source":"patient_probs.to_csv('submission.csv', header = True, index = False)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T12:37:32.81413Z","iopub.execute_input":"2023-10-05T12:37:32.814947Z","iopub.status.idle":"2023-10-05T12:37:32.822987Z","shell.execute_reply.started":"2023-10-05T12:37:32.814915Z","shell.execute_reply":"2023-10-05T12:37:32.821971Z"},"trusted":true},"execution_count":null,"outputs":[]}]}