{"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":"# RSNA-2023-1st-Place-Best-Model-Infer (Cleaned)\n\nGreetings,\n\nThis is a notebook that simplifies our final inference ensemble pipeline.\n\nIn this notebook, for each stage, we only used a single model -- No ensemble here anymore.\n\nAlso the code has been cleaned up for demonstrating how to predicting on a single patient.\n\n---\n\n\nOur brief summary of winning solution: https://www.kaggle.com/competitions/rsna-2023-abdominal-trauma-detection/discussion/447449\n        \nNotebook for training 3d semantic segmentation: https://www.kaggle.com/code/haqishen/rsna-2023-1st-place-solution-train-3d-seg\n\n\nIf you find these notebooks helpful please upvote. Thanks!","metadata":{}},{"cell_type":"code","source":"!cp -r '/kaggle/input/contrails-libraries/pretrainedmodels-0.7.4/' './'\n!cp -r '/kaggle/input/contrails-libraries/efficientnet_pytorch-0.7.1/' './'\n\n!pip -q install /kaggle/input/dicomsdl--0-109-2/dicomsdl-0.109.2-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl\n!pip -q install '/kaggle/input/contrails-libraries/segmentation_models_pytorch-0.3.3-py3-none-any.whl' --no-deps\n!pip -q install /kaggle/input/contrails-model-def1/einops-0.6.1-py3-none-any.whl","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:13:36.889447Z","iopub.execute_input":"2023-10-19T08:13:36.889774Z","iopub.status.idle":"2023-10-19T08:13:58.16416Z","shell.execute_reply.started":"2023-10-19T08:13:36.889741Z","shell.execute_reply":"2023-10-19T08:13:58.163147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.append('./pretrainedmodels-0.7.4/pretrainedmodels-0.7.4/')\nsys.path.append('./efficientnet_pytorch-0.7.1/efficientnet_pytorch-0.7.1/')\nsys.path.append(\"/kaggle/input/rsna-abd-models-classes/\")\n","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:13:58.165717Z","iopub.execute_input":"2023-10-19T08:13:58.166014Z","iopub.status.idle":"2023-10-19T08:13:58.170749Z","shell.execute_reply.started":"2023-10-19T08:13:58.165989Z","shell.execute_reply":"2023-10-19T08:13:58.169922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport gc\nimport copy\nimport time\nimport numpy as np\nimport pandas as pd\nfrom glob import glob\nfrom tqdm import tqdm\n\nimport cv2\nfrom PIL import Image\nimport pydicom\nimport albumentations as A\nfrom albumentations.pytorch import ToTensorV2\nimport matplotlib.pyplot as plt\n\nimport torch\nfrom torch import nn\nimport torch.nn.functional as F\n\nimport timm\nimport segmentation_models_pytorch as smp\nfrom models import *\n\nimport dicomsdl\ndef __dataset__to_numpy_image(self, index=0):\n    info = self.getPixelDataInfo()\n    dtype = info['dtype']\n    if info['SamplesPerPixel'] != 1:\n        raise RuntimeError('SamplesPerPixel != 1')\n    else:\n        shape = [info['Rows'], info['Cols']]\n    outarr = np.empty(shape, dtype=dtype)\n    self.copyFrameData(index, outarr)\n    return outarr\ndicomsdl._dicomsdl.DataSet.to_numpy_image = __dataset__to_numpy_image   \n\n\ntorch.cuda.set_device('cuda:0')","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:13:58.171679Z","iopub.execute_input":"2023-10-19T08:13:58.17188Z","iopub.status.idle":"2023-10-19T08:14:05.285802Z","shell.execute_reply.started":"2023-10-19T08:13:58.171863Z","shell.execute_reply":"2023-10-19T08:14:05.284765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocess Util","metadata":{}},{"cell_type":"code","source":"def glob_sorted(path):\n    return sorted(glob(path), key=lambda x: int(x.split('/')[-1].split('.')[0]))\n\ndef get_rescaled_image(dcm, img):\n    resI, resS = dcm.RescaleIntercept, dcm.RescaleSlope\n    img = resS * img + resI\n    return img\n\ndef get_windowed_image(img, WL=50, WW=400):\n    upper, lower = WL+WW//2, WL-WW//2\n    X = np.clip(img.copy(), lower, upper)\n    X = X - np.min(X)\n    X = X / np.max(X)\n    X = (X*255.0).astype('uint8')\n    \n    return X\n\ndef standardize_pixel_array(dcm, pixel_array):\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    #pixel_array = dcm.pixel_array\n    \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\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\n\ndef load_volume(dcms):\n    volume = []\n    pos_zs = []\n    \n    for dcm_path in dcms:\n        pydcm = pydicom.dcmread(dcm_path)\n        \n        pos_z = pydcm[(0x20, 0x32)].value[-1]\n        pos_zs.append(pos_z)\n        \n        dcm = dicomsdl.open(dcm_path)\n        \n        orig_image = dcm.to_numpy_image()\n        image = get_rescaled_image(dcm, orig_image)\n        image = get_windowed_image(image)\n        \n        if np.min(image)<0:\n            image = image + np.abs(np.min(image))\n        \n        image = image / image.max()\n        image = (image * 255).astype(np.uint8)\n        volume.append(image)\n    \n    return np.stack(volume)\n\n\ndef process_volume(volume):\n    volume = np.stack([cv2.resize(x, (128, 128)) for x in volume])\n    \n    volumes = []\n    cuts = [(x, x+32) for x in np.arange(0, volume.shape[0], 32)[:-1]]\n    \n    if cuts:\n        for cut in cuts:\n            volumes.append(volume[cut[0]:cut[1]])\n        volumes = np.stack(volumes)\n    else:\n        volumes = np.zeros((1, 32, 128, 128), dtype=np.uint8)\n        volumes[0, :len(volume)] = volume\n    \n    if cuts:\n        last_volume = np.zeros((1, 32, 128, 128), dtype=np.uint8)\n        last_volume[0, :volume[cuts[-1][1]:].shape[0]] =  volume[cuts[-1][1]:]\n        volumes = np.concatenate([volumes, last_volume])\n    \n    volumes = torch.as_tensor(volumes).float()\n    \n    return volumes\n\n\ndef get_volume_data(grd, step=96, stride=1, stride_cutoff=200):\n    volumes = []\n    \n    if len(grd)>stride_cutoff:\n        grd = grd[::stride]\n\n    take_last = False\n    if not str(len(grd)/step).endswith('.0'):\n        take_last = True\n\n    started = False\n    for i in range(len(grd)//step):\n        rows = grd[i*step:(i+1)*step]\n\n        if len(rows)!=step:\n            rows = pd.DataFrame([rows.iloc[int(x*len(rows))] for x in np.arange(0, 1, 1/step)])\n\n        volumes.append(rows)\n\n        started = True\n\n    if not started:\n        rows = grd\n        rows = pd.DataFrame([rows.iloc[int(x*len(rows))] for x in np.arange(0, 1, 1/step)])\n        volumes.append(rows)\n\n    if take_last:\n        rows = grd[-step:]\n        if len(rows)==step:\n            volumes.append(rows)\n\n    return volumes","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:05.287755Z","iopub.execute_input":"2023-10-19T08:14:05.288049Z","iopub.status.idle":"2023-10-19T08:14:05.302886Z","shell.execute_reply.started":"2023-10-19T08:14:05.288026Z","shell.execute_reply":"2023-10-19T08:14:05.302039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"IMAGE_FOLDER = '/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/'\n\npatient = '26501'  # We only predict a single patient in this notebook\n\ntest_augs = A.Compose([\n    A.Resize(384, 384),\n    ToTensorV2()\n])\n\npatient","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:05.303928Z","iopub.execute_input":"2023-10-19T08:14:05.3042Z","iopub.status.idle":"2023-10-19T08:14:05.319871Z","shell.execute_reply.started":"2023-10-19T08:14:05.304154Z","shell.execute_reply":"2023-10-19T08:14:05.319191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"code","source":"path = f'/kaggle/input/rsna-abd-models/try3_seg_resnet18d_v3/zip/0.pth'\nst = torch.load(path, map_location='cpu')\nmodel_3dseg = convert_3d(SegmentationModel())\nmodel_3dseg.load_state_dict(st)\nmodel_3dseg.eval()\nmodel_3dseg.cuda()\n    \n    \npath = f\"/kaggle/input/coatmed384ourdataseed6969/3.pth\"\nst = torch.load(path, map_location='cpu')\nmodel_organs = Model4(num_classes=10, seg_classes=4, arch='medium', mask_head=False)\nmodel_organs.load_state_dict(st)\nmodel_organs.cuda()\nmodel_organs.eval()\n\n\npath = f\"/kaggle/input/coatsmall384extravast4funet/3.pth\"\nst = torch.load(path, map_location='cpu')\nmodel_extrav = Model4(num_classes=2, seg_classes=4, arch='small', mask_head=False)\nmodel_extrav.load_state_dict(st)\nmodel_extrav.cuda()\nmodel_extrav.eval()\n\nprint('''We load 3 models here:\n1: 3d semantic segmentation model for segment organs\n2: 2.5d classification model for classify organs\n3: 2.5d classification model for classify extravasation\n''')","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:05.321076Z","iopub.execute_input":"2023-10-19T08:14:05.32131Z","iopub.status.idle":"2023-10-19T08:14:15.940016Z","shell.execute_reply.started":"2023-10-19T08:14:05.32129Z","shell.execute_reply":"2023-10-19T08:14:15.939104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predict","metadata":{}},{"cell_type":"code","source":"PATIENT_TO_PREDICTION = {}\nPATIENT_TO_PREDICTION2 = {}\n\nfinal_outputs = []\nfinal_outputs2 = []\n\nstudies = os.listdir(f'{IMAGE_FOLDER}/{patient}')\nfor study in studies:\n\n    files = glob_sorted(f\"{IMAGE_FOLDER}/{patient}/{study}/*\")\n\n    volume = load_volume(files)\n    file_to_volume = {file: vol for file, vol in zip(files, volume)}\n\n    volumes = process_volume(volume)\n    volumes_seg = predict_segmentation(volumes, [model_3dseg])\n    volume_seg = np.concatenate(volumes_seg.transpose(0, 2, 1, 3, 4))[:len(volume)]\n    \n    vis_seg_0 = volumes[volumes.shape[0]//2, 16].numpy().astype(np.uint8)\n    vis_seg_1 = volumes_seg[volumes_seg.shape[0]//2, :, 16]\n    vis_seg_1[0] += vis_seg_1[3]\n    vis_seg_1[1] += vis_seg_1[4]\n    vis_seg_1 = (vis_seg_1[:3].transpose(1,2,0).clip(0, 1) * 255).astype(np.uint8)\n\n#     print(volumes.shape, volumes_seg.shape)  # torch.Size([7, 32, 128, 128]) (7, 5, 32, 128, 128)\n\n    msk = volume_seg.max(0).max(0)\n    ys, xs = np.where(msk)\n    y1, y2, x1, x2 = np.min(ys) / 128, np.max(ys) / 128, np.min(xs) / 128, np.max(xs) / 128\n\n    files = pd.DataFrame({\"file\": files})\n    files_volumes = get_volume_data(files, step=96, stride=2, stride_cutoff=400)\n\n    first = True\n\n    del volumes, volumes_seg, volume_seg, volume\n    gc.collect()\n\n    for file_volume in files_volumes:\n        volume = np.stack([file_to_volume[file] for file in file_volume.file])\n\n        if first:\n            h, w = volume.shape[1:]\n            y1, y2, x1, x2 = int(y1*h), int(y2*h), int(x1*w), int(x2*w)\n        volume2 = volume\n        \n        #### CROPPED #####\n        volume = volume[:, y1:y2, x1:x2]\n\n        vols = []\n        NC = 3\n        for i in range(len(volume)//NC):\n            vols.append(volume[i*NC:(i+1)*NC])\n        vol = np.stack(vols, 0).transpose(0, 2, 3, 1)\n\n        volume_ = []\n        for image in vol:\n            image = image.astype(np.float32) / 255\n            transformed = test_augs(image=image)\n            image = transformed['image']\n            volume_.append(image)\n        volume = torch.stack(volume_).float()\n        volume = volume.cuda()\n\n        #### UNCROPPED #####\n        vols = []\n        NC = 3\n        for i in range(len(volume2)//NC):\n            vols.append(volume2[i*NC:(i+1)*NC])\n        vol = np.stack(vols, 0).transpose(0, 2, 3, 1)\n\n        volume_ = []\n        for image in vol:\n            image = image.astype(np.float32) / 255\n            transformed = test_augs(image=image)\n            image = transformed['image']\n            volume_.append(image)\n\n        volume2 = torch.stack(volume_).float()\n        volume2 = volume2.cuda()\n\n        outputs = []\n        outputs2 = []\n\n        with torch.no_grad():\n            with torch.cuda.amp.autocast(enabled=True):\n\n                outs = model_organs(volume.unsqueeze(0))\n                outs = outs.float().sigmoid()\n                outputs.append(outs)\n\n                outs = model_extrav(volume2.unsqueeze(0))[:, :, [1, 0]]\n                outs = outs.float().sigmoid()\n                outputs2.append(outs)\n\n        torch.cuda.empty_cache()\n\n        outputs = torch.stack(outputs)[:, 0].mean(0)\n        outputs2 = torch.stack(outputs2)[:, 0].mean(0)\n\n        final_outputs.append(outputs.detach().cpu().numpy())\n        final_outputs2.append(outputs2.detach().cpu().numpy())\n\n        first = False\n\n        torch.cuda.empty_cache()\n\n\nlast_final_outputs = final_outputs.copy()\nlast_final_outputs2 = final_outputs2.copy()\n\nfinal_outputs = np.concatenate(final_outputs)\nfinal_outputs2 = np.concatenate(final_outputs2)\n\nfinal_predictions = final_outputs.max(0)\nfinal_predictions2 = final_outputs2.max(0)\n\nPATIENT_TO_PREDICTION[patient] = final_predictions\nPATIENT_TO_PREDICTION2[patient] = final_predictions2\n","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:15.941586Z","iopub.execute_input":"2023-10-19T08:14:15.941848Z","iopub.status.idle":"2023-10-19T08:14:29.587434Z","shell.execute_reply.started":"2023-10-19T08:14:15.941825Z","shell.execute_reply":"2023-10-19T08:14:29.586736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axarr = plt.subplots(1, 2, figsize=(10, 4))\n\naxarr[0].imshow(vis_seg_0)\naxarr[0].axis('off') \n\naxarr[1].imshow(vis_seg_1)\naxarr[1].axis('off') \n\nplt.tight_layout() \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:29.588352Z","iopub.execute_input":"2023-10-19T08:14:29.588609Z","iopub.status.idle":"2023-10-19T08:14:29.864257Z","shell.execute_reply.started":"2023-10-19T08:14:29.588589Z","shell.execute_reply":"2023-10-19T08:14:29.863449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Postprocess","metadata":{}},{"cell_type":"code","source":"bowel_w = 2\nextrav_w = 6\nlow_w = 2\nhigh_w = 4\n\nFINAL_SUB = {'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\nfor patient in PATIENT_TO_PREDICTION:\n    prediction = PATIENT_TO_PREDICTION[patient].copy()\n    prediction2 = PATIENT_TO_PREDICTION2[patient].copy()\n    \n    prediction[9] = (prediction[9] * 0.666) + (prediction2[1]*0.334)\n    \n    FINAL_SUB['patient_id'].append(patient)\n    \n    FINAL_SUB['bowel_healthy'].append(1 - prediction[9])\n    FINAL_SUB['bowel_injury'].append(prediction[9] * bowel_w)\n    FINAL_SUB['extravasation_healthy'].append(1 - prediction2[0])\n    FINAL_SUB['extravasation_injury'].append(0.06355258976803305 + (prediction2[0] * extrav_w))\n    \n    FINAL_SUB['liver_healthy'].append(1 - prediction[0])\n    FINAL_SUB['liver_low'].append(prediction[3]*low_w)\n    FINAL_SUB['liver_high'].append(prediction[4]*high_w)\n\n    FINAL_SUB['spleen_healthy'].append(1 - prediction[1])\n    FINAL_SUB['spleen_low'].append(prediction[5]*low_w)\n    FINAL_SUB['spleen_high'].append(prediction[6]*high_w)\n\n    FINAL_SUB['kidney_healthy'].append(1 - prediction[2])\n    FINAL_SUB['kidney_low'].append((prediction[7])*low_w)\n    FINAL_SUB['kidney_high'].append(prediction[8]*high_w)\n","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:29.865584Z","iopub.execute_input":"2023-10-19T08:14:29.865953Z","iopub.status.idle":"2023-10-19T08:14:29.876668Z","shell.execute_reply.started":"2023-10-19T08:14:29.86592Z","shell.execute_reply":"2023-10-19T08:14:29.875596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nsubmission = pd.DataFrame(FINAL_SUB)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:29.87861Z","iopub.execute_input":"2023-10-19T08:14:29.878888Z","iopub.status.idle":"2023-10-19T08:14:29.905381Z","shell.execute_reply.started":"2023-10-19T08:14:29.878866Z","shell.execute_reply":"2023-10-19T08:14:29.904634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_id = submission['patient_id'].iloc[0]\ndata_to_plot = submission.drop(columns=['patient_id'])\nenglish_column_names = {\n    \"bowel_healthy\": \"Bowel Healthy\",\n    \"bowel_injury\": \"Bowel Injury\",\n    \"extravasation_healthy\": \"Extravasation Healthy\",\n    \"extravasation_injury\": \"Extravasation Injury\",\n    \"kidney_healthy\": \"Kidney Healthy\",\n    \"kidney_low\": \"Kidney Low\",\n    \"kidney_high\": \"Kidney High\",\n    \"liver_healthy\": \"Liver Healthy\",\n    \"liver_low\": \"Liver Low\",\n    \"liver_high\": \"Liver High\",\n    \"spleen_healthy\": \"Spleen Healthy\",\n    \"spleen_low\": \"Spleen Low\",\n    \"spleen_high\": \"Spleen High\"\n}\n\ndata_to_plot.columns = [english_column_names[col] for col in data_to_plot.columns]\n\nplt.figure(figsize=(15, 7))\ndata_to_plot.T.plot(kind='bar', legend=False, ax=plt.gca())\nplt.title(f\"Probability of Diseases for Patient {patient_id}\")\nplt.ylabel(\"Probability\")\nplt.xticks(rotation=45, ha='right')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:22:14.238435Z","iopub.execute_input":"2023-10-19T08:22:14.23876Z","iopub.status.idle":"2023-10-19T08:22:14.619042Z","shell.execute_reply.started":"2023-10-19T08:22:14.238734Z","shell.execute_reply":"2023-10-19T08:22:14.61821Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -rf /kaggle/working/*","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('./submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:14:30.944712Z","iopub.execute_input":"2023-10-19T08:14:30.945086Z","iopub.status.idle":"2023-10-19T08:14:30.956847Z","shell.execute_reply.started":"2023-10-19T08:14:30.945055Z","shell.execute_reply":"2023-10-19T08:14:30.956081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}