{"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 2022 Cervical Spine Fracture Detection : Ensembling best solutions\n\n\nOriginal Notebooks : \n\n- [[infer] PyTorch-EffNetV2 single-model LB:0.49](https://www.kaggle.com/code/vslaykovsky/infer-pytorch-effnetv2-single-model-lb-0-49/data?scriptVersionId=104434908)\n\n- [RNSA - 3D model [Infer] [PyTorch]](https://www.kaggle.com/code/samuelcortinhas/rnsa-3d-model-infer-pytorch/data)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-info\" style=\"text-align:center; font-size:28px;\">\n    Solution 1 : [infer] PyTorch-EffNetV2 single-model LB:0.49\n</div>","metadata":{}},{"cell_type":"markdown","source":"### 1. Imports, constants, dependencies","metadata":{}},{"cell_type":"code","source":"try:\n    import pylibjpeg\nexcept:\n    # Offline dependencies:\n    !mkdir -p /root/.cache/torch/hub/checkpoints/\n    !cp ../input/rsna-2022-whl/efficientnet_v2_s-dd5fe13b.pth  /root/.cache/torch/hub/checkpoints/\n\n    !pip install /kaggle/input/rsna-2022-whl/{pydicom-2.3.0-py3-none-any.whl,pylibjpeg-1.4.0-py3-none-any.whl,python_gdcm-3.0.15-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl}\n    !pip install /kaggle/input/rsna-2022-whl/{torch-1.12.1-cp37-cp37m-manylinux1_x86_64.whl,torchvision-0.13.1-cp37-cp37m-manylinux1_x86_64.whl}","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:20:12.510071Z","iopub.execute_input":"2022-10-27T20:20:12.511558Z","iopub.status.idle":"2022-10-27T20:20:12.529787Z","shell.execute_reply.started":"2022-10-27T20:20:12.511511Z","shell.execute_reply":"2022-10-27T20:20:12.528694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\nimport glob\nimport os\nimport re\n\nimport cv2\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pydicom as dicom\nimport torch\nimport torchvision as tv\nfrom sklearn.model_selection import GroupKFold\nfrom torch.cuda.amp import GradScaler, autocast\nfrom torchvision.models.feature_extraction import create_feature_extractor\nfrom tqdm.notebook import tqdm\n\nimport wandb\n\npd.set_option('display.max_rows', 1000)\npd.set_option('display.max_columns', 1000)\nplt.rcParams['figure.figsize'] = (20, 5)\n\n\n# Effnet\nWEIGHTS = tv.models.efficientnet.EfficientNet_V2_S_Weights.DEFAULT\nRSNA_2022_PATH = '../input/rsna-2022-cervical-spine-fracture-detection'\nTRAIN_IMAGES_PATH = f'{RSNA_2022_PATH}/train_images'\nTEST_IMAGES_PATH = f'{RSNA_2022_PATH}/test_images'\nEFFNET_CHECKPOINTS_PATH = '../input/rsna-2022-base-effnetv2'\n\n# MODEL_NAMES = [f'effnetv2']\n\n# This notebook supports ensembles and single model predictions. Uncomment to switch to ensemble prediction:\nMODEL_NAMES = [f'effnetv2-f{i}' for i in range(5)]\n\n# Common\nFRAC_COLS = [f'C{i}_effnet_frac' for i in range(1, 8)]\nVERT_COLS = [f'C{i}_effnet_vert' for i in range(1, 8)]\n\ntry:\n    from kaggle_secrets import UserSecretsClient\n    IS_KAGGLE = True\nexcept:\n    IS_KAGGLE = False\n\n\n# Switch to offline for submission\nos.environ[\"WANDB_MODE\"] = \"offline\"\n\nif os.environ[\"WANDB_MODE\"] == \"online\":\n    if IS_KAGGLE:\n        os.environ['WANDB_API_KEY'] = UserSecretsClient().get_secret(\"WANDB_API_KEY\")\n\nif not IS_KAGGLE:\n    print('Running locally')\n    RSNA_2022_PATH = '/mnt/rsna2022'\n    TRAIN_IMAGES_PATH = '/mnt/rsna2022/train_images'\n    TEST_IMAGES_PATH = '/mnt/rsna2022/test_images'\n    METADATA_PATH = '/home/vslaykovsky/Downloads/'\n    EFFNET_CHECKPOINTS_PATH = 'frac_checkpoints'\n    os.environ['WANDB_API_KEY'] = 'yourkeyhere'\n\n\nDEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'\nif DEVICE == 'cuda':\n    BATCH_SIZE = 32\nelse:\n    BATCH_SIZE = 2","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-10-27T20:20:12.536387Z","iopub.execute_input":"2022-10-27T20:20:12.540225Z","iopub.status.idle":"2022-10-27T20:20:12.559154Z","shell.execute_reply.started":"2022-10-27T20:20:12.54018Z","shell.execute_reply":"2022-10-27T20:20:12.557738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2. Loading train/eval/test dataframes","metadata":{}},{"cell_type":"code","source":"def load_df_test():\n    df_test = pd.read_csv(f'{RSNA_2022_PATH}/test.csv')\n\n    if df_test.iloc[0].row_id == '1.2.826.0.1.3680043.10197_C1':\n        # test_images and test.csv are inconsistent in the dev dataset, fixing labels for the dev run.\n        df_test = pd.DataFrame({\n            \"row_id\": ['1.2.826.0.1.3680043.22327_C1', '1.2.826.0.1.3680043.25399_C1', '1.2.826.0.1.3680043.5876_C1'],\n            \"StudyInstanceUID\": ['1.2.826.0.1.3680043.22327', '1.2.826.0.1.3680043.25399', '1.2.826.0.1.3680043.5876'],\n            \"prediction_type\": [\"C1\", \"C1\", \"patient_overall\"]}\n        )\n    return df_test\n\ndf_test = load_df_test()\ndf_test","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:20:12.566512Z","iopub.execute_input":"2022-10-27T20:20:12.569329Z","iopub.status.idle":"2022-10-27T20:20:12.598235Z","shell.execute_reply.started":"2022-10-27T20:20:12.569278Z","shell.execute_reply":"2022-10-27T20:20:12.597095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_slices = glob.glob(f'{TEST_IMAGES_PATH}/*/*')\ntest_slices = [re.findall(f'{TEST_IMAGES_PATH}/(.*)/(.*).dcm', s)[0] for s in test_slices]\ndf_test_slices = pd.DataFrame(data=test_slices, columns=['StudyInstanceUID', 'Slice']).astype({'Slice': int}).sort_values(['StudyInstanceUID', 'Slice']).reset_index(drop=True)\ndf_test_slices","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:20:12.603689Z","iopub.execute_input":"2022-10-27T20:20:12.606247Z","iopub.status.idle":"2022-10-27T20:20:12.649333Z","shell.execute_reply.started":"2022-10-27T20:20:12.606204Z","shell.execute_reply":"2022-10-27T20:20:12.648218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3. Dataset class","metadata":{}},{"cell_type":"code","source":"def load_dicom(path):\n    \"\"\"\n    This supports loading both regular and compressed JPEG images. \n    See the first sell with `pip install` commands for the necessary dependencies\n    \"\"\"\n    img=dicom.dcmread(path)\n    img.PhotometricInterpretation = 'YBR_FULL'\n    data = img.pixel_array    \n    data = data - np.min(data)\n    if np.max(data) != 0:\n        data = data / np.max(data)\n    data=(data * 255).astype(np.uint8)\n    return cv2.cvtColor(data, cv2.COLOR_GRAY2RGB), img\n\n\nim, meta = load_dicom(f'{TRAIN_IMAGES_PATH}/1.2.826.0.1.3680043.10001/1.dcm')\nplt.figure()\nplt.imshow(im)\nplt.title('regular image')\n\nim, meta = load_dicom(f'{TRAIN_IMAGES_PATH}/1.2.826.0.1.3680043.10014/1.dcm')\nplt.figure()\nplt.imshow(im)\nplt.title('jpeg')","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:20:12.780406Z","iopub.execute_input":"2022-10-27T20:20:12.781703Z","iopub.status.idle":"2022-10-27T20:20:13.405387Z","shell.execute_reply.started":"2022-10-27T20:20:12.781652Z","shell.execute_reply":"2022-10-27T20:20:13.404418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class EffnetDataSet(torch.utils.data.Dataset):    \n    def __init__(self, df, path, transforms=None):\n        super().__init__()\n        self.df = df\n        self.path = path\n        self.transforms = transforms\n        \n    def __getitem__(self, i):\n        path = os.path.join(self.path, self.df.iloc[i].StudyInstanceUID, f'{self.df.iloc[i].Slice}.dcm')        \n        \n        try:\n            img = load_dicom(path)[0]         \n            img = np.transpose(img, (2, 0, 1))  # Pytorch uses (batch, channel, height, width) order. Converting (height, width, channel) -> (channel, height, width)\n            if self.transforms is not None:\n                img = self.transforms(torch.as_tensor(img))\n        except Exception as ex:\n            print(ex)\n            return None\n        \n        if 'C1_fracture' in self.df:\n            frac_targets = torch.as_tensor(self.df.iloc[i][['C1_fracture', 'C2_fracture', 'C3_fracture', 'C4_fracture', 'C5_fracture', 'C6_fracture', 'C7_fracture']].astype('float32').values)\n            vert_targets = torch.as_tensor(self.df.iloc[i][['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7']].astype('float32').values)\n            frac_targets = frac_targets * vert_targets   # we only enable targets that are visible on the current slice\n            return img, frac_targets, vert_targets\n        return img        \n    \n    def __len__(self):\n        return len(self.df)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:20:13.407539Z","iopub.execute_input":"2022-10-27T20:20:13.408552Z","iopub.status.idle":"2022-10-27T20:20:13.419645Z","shell.execute_reply.started":"2022-10-27T20:20:13.408506Z","shell.execute_reply":"2022-10-27T20:20:13.418348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Only X values returned by the test dataset\nds_test = EffnetDataSet(df_test_slices, TEST_IMAGES_PATH, WEIGHTS.transforms())\nX = ds_test[42]\nX.shape","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:13.422584Z","iopub.execute_input":"2022-10-27T20:20:13.423033Z","iopub.status.idle":"2022-10-27T20:20:13.451773Z","shell.execute_reply.started":"2022-10-27T20:20:13.422991Z","shell.execute_reply":"2022-10-27T20:20:13.450647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class EffnetModel(torch.nn.Module):\n    def __init__(self):\n        super().__init__()\n        effnet = tv.models.efficientnet_v2_s()\n        self.model = create_feature_extractor(effnet, ['flatten'])\n        self.nn_fracture = torch.nn.Sequential(\n            torch.nn.Linear(1280, 7),\n        )\n        self.nn_vertebrae = torch.nn.Sequential(\n            torch.nn.Linear(1280, 7),\n        )\n\n    def forward(self, x):\n        # returns logits\n        x = self.model(x)['flatten']\n        return self.nn_fracture(x), self.nn_vertebrae(x)\n\n    def predict(self, x):\n        frac, vert = self.forward(x)\n        return torch.sigmoid(frac), torch.sigmoid(vert)\n\nmodel = EffnetModel()\nmodel.predict(torch.randn(1, 3, 512, 512))\ndel model","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:13.454893Z","iopub.execute_input":"2022-10-27T20:20:13.455294Z","iopub.status.idle":"2022-10-27T20:20:15.267397Z","shell.execute_reply.started":"2022-10-27T20:20:13.455255Z","shell.execute_reply":"2022-10-27T20:20:15.266377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_model(model, name, path='.'):\n    data = torch.load(os.path.join(path, f'{name}.tph'), map_location=DEVICE)\n    model.load_state_dict(data)\n    return model","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:15.268882Z","iopub.execute_input":"2022-10-27T20:20:15.269943Z","iopub.status.idle":"2022-10-27T20:20:15.276265Z","shell.execute_reply.started":"2022-10-27T20:20:15.2699Z","shell.execute_reply":"2022-10-27T20:20:15.275065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"effnet_models = [load_model(EffnetModel(), name, EFFNET_CHECKPOINTS_PATH).to(DEVICE) for name in MODEL_NAMES]","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:15.278013Z","iopub.execute_input":"2022-10-27T20:20:15.278455Z","iopub.status.idle":"2022-10-27T20:20:20.12112Z","shell.execute_reply.started":"2022-10-27T20:20:15.278416Z","shell.execute_reply":"2022-10-27T20:20:20.120094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 7. Submission","metadata":{}},{"cell_type":"code","source":"from typing import List\n\n\ndef predict_effnet(models: List[EffnetModel], ds, max_batches=1e9):\n    dl_test = torch.utils.data.DataLoader(ds, batch_size=BATCH_SIZE, shuffle=False, num_workers=os.cpu_count())\n    for m in models:\n        m.eval()\n\n    with torch.no_grad():\n        predictions = []\n        for idx, X in enumerate(tqdm(dl_test, miniters=10)):\n            pred = torch.zeros(len(X), 14).to(DEVICE)\n            for m in models:\n                y1, y2 = m.predict(X.to(DEVICE))\n                pred += torch.concat([y1, y2], dim=1) / len(models)\n            predictions.append(pred)\n            if idx >= max_batches:\n                break\n        return torch.concat(predictions).cpu().numpy()\n\n# Quick test\npredict_effnet([EffnetModel().to(DEVICE)], ds_test, max_batches=2).shape","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:20.124292Z","iopub.execute_input":"2022-10-27T20:20:20.124611Z","iopub.status.idle":"2022-10-27T20:20:23.171149Z","shell.execute_reply.started":"2022-10-27T20:20:20.124571Z","shell.execute_reply":"2022-10-27T20:20:23.169868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"effnet_pred = predict_effnet(effnet_models, ds_test)\n\ndf_effnet_pred = pd.DataFrame(\n    data=effnet_pred, columns=[f'C{i}_effnet_frac' for i in range(1, 8)] + [f'C{i}_effnet_vert' for i in range(1, 8)]\n)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:20:23.173043Z","iopub.execute_input":"2022-10-27T20:20:23.173847Z","iopub.status.idle":"2022-10-27T20:21:00.201847Z","shell.execute_reply.started":"2022-10-27T20:20:23.173799Z","shell.execute_reply":"2022-10-27T20:21:00.200649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test_pred = pd.concat([df_test_slices, df_effnet_pred], axis=1).sort_values(['StudyInstanceUID', 'Slice'])\ndf_test_pred","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:00.203818Z","iopub.execute_input":"2022-10-27T20:21:00.204508Z","iopub.status.idle":"2022-10-27T20:21:00.236218Z","shell.execute_reply.started":"2022-10-27T20:21:00.204458Z","shell.execute_reply":"2022-10-27T20:21:00.235167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_sample_patient(df_pred):\n    patient = np.random.choice(df_pred.StudyInstanceUID)\n    df = df_pred.query('StudyInstanceUID == @patient').reset_index()\n\n    df[[f'C{i}_effnet_frac' for i in range(1, 8)]].plot(\n        title=f'Patient {patient}, fracture prediction',\n        ax=(plt.subplot(1, 2, 1)))\n\n    df[[f'C{i}_effnet_vert' for i in range(1, 8)]].plot(\n        title=f'Patient {patient}, vertebrae prediction',\n        ax=plt.subplot(1, 2, 2)\n    )\n\nplot_sample_patient(df_test_pred)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:00.240571Z","iopub.execute_input":"2022-10-27T20:21:00.240897Z","iopub.status.idle":"2022-10-27T20:21:01.125467Z","shell.execute_reply.started":"2022-10-27T20:21:00.240868Z","shell.execute_reply":"2022-10-27T20:21:01.12449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def patient_prediction(df):\n    c1c7 = np.average(df[FRAC_COLS].values, axis=0, weights=df[VERT_COLS].values)\n    pred_patient_overall = 1 - np.prod(1 - c1c7)\n    return pd.Series(data=np.concatenate([[pred_patient_overall], c1c7]), index=['patient_overall'] + [f'C{i}' for i in range(1, 8)])\n\ndf_patient_pred = df_test_pred.groupby('StudyInstanceUID').apply(lambda df: patient_prediction(df))\ndf_patient_pred","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:01.129583Z","iopub.execute_input":"2022-10-27T20:21:01.132202Z","iopub.status.idle":"2022-10-27T20:21:01.166198Z","shell.execute_reply.started":"2022-10-27T20:21:01.132161Z","shell.execute_reply":"2022-10-27T20:21:01.164579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sub = df_test.copy()\ndf_sub = df_sub.set_index('StudyInstanceUID').join(df_patient_pred)\ndf_sub['fractured'] = df_sub.apply(lambda r: r[r.prediction_type], axis=1)\n#df_sub[['row_id', 'fractured']].to_csv('submission1.csv', index=False)\ndf_sub","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:01.170383Z","iopub.execute_input":"2022-10-27T20:21:01.173036Z","iopub.status.idle":"2022-10-27T20:21:01.201844Z","shell.execute_reply.started":"2022-10-27T20:21:01.172994Z","shell.execute_reply":"2022-10-27T20:21:01.200863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-info\" style=\"text-align:center; font-size:28px;\">\n    Solution 2 : RNSA - 3D model [Infer] [PyTorch]\n</div>","metadata":{}},{"cell_type":"code","source":"!pip install -qU ../input/for-pydicom/python_gdcm-3.0.14-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl ../input/for-pydicom/pylibjpeg-1.4.0-py3-none-any.whl --find-links frozen_packages --no-index","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:01.206127Z","iopub.execute_input":"2022-10-27T20:21:01.20887Z","iopub.status.idle":"2022-10-27T20:21:11.311507Z","shell.execute_reply.started":"2022-10-27T20:21:01.208825Z","shell.execute_reply":"2022-10-27T20:21:11.310241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -q kaggle_vol3d_classify -f ../input/cervical-spine-fracture-detection-npz-3d-volumes/frozen_packages --no-index","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:11.313548Z","iopub.execute_input":"2022-10-27T20:21:11.314362Z","iopub.status.idle":"2022-10-27T20:21:21.778184Z","shell.execute_reply.started":"2022-10-27T20:21:11.314313Z","shell.execute_reply":"2022-10-27T20:21:21.776801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport matplotlib.patches as patches\nimport seaborn as sns\nsns.set(style='darkgrid', font_scale=1.6)\nimport cv2\nimport os\nfrom os import listdir\nimport re\nimport gc\nimport random\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom tqdm.auto import tqdm\nfrom pprint import pprint\nfrom time import time\nimport itertools\nfrom skimage import measure\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport nibabel as nib\nfrom glob import glob\nimport warnings\n#warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n#warnings.filterwarnings(\"ignore\", category=UserWarning)\n#warnings.filterwarnings(\"ignore\", category=FutureWarning)\nimport zipfile\nfrom scipy import ndimage\nfrom sklearn.model_selection import train_test_split\nfrom joblib import Parallel, delayed\nfrom PIL import Image\nfrom dipy.denoise.nlmeans import nlmeans\nfrom dipy.denoise.noise_estimate import estimate_sigma\nfrom kaggle_volclassif.utils import interpolate_volume\nfrom skimage import exposure\n\n# Pytorch\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.optim.lr_scheduler as lr_scheduler\nfrom torch.utils.data import Dataset, DataLoader\nimport torchvision.transforms as transforms\nimport torch.nn.functional as F","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-10-27T20:21:21.7803Z","iopub.execute_input":"2022-10-27T20:21:21.780739Z","iopub.status.idle":"2022-10-27T20:21:21.801383Z","shell.execute_reply.started":"2022-10-27T20:21:21.7807Z","shell.execute_reply":"2022-10-27T20:21:21.800274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Reproducibility","metadata":{}},{"cell_type":"code","source":"# Set random seeds\ndef set_seed(seed=0):\n    np.random.seed(seed)\n    random.seed(seed)\n    torch.manual_seed(seed)\nset_seed()","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:21.803069Z","iopub.execute_input":"2022-10-27T20:21:21.803949Z","iopub.status.idle":"2022-10-27T20:21:21.813971Z","shell.execute_reply.started":"2022-10-27T20:21:21.803902Z","shell.execute_reply":"2022-10-27T20:21:21.812913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Config","metadata":{}},{"cell_type":"code","source":"# Hyperparameters\nBATCH_SIZE = 1\n\n# Config device\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\ndevice","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:21.817219Z","iopub.execute_input":"2022-10-27T20:21:21.817531Z","iopub.status.idle":"2022-10-27T20:21:21.828571Z","shell.execute_reply.started":"2022-10-27T20:21:21.817504Z","shell.execute_reply":"2022-10-27T20:21:21.827608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load tables\n","metadata":{}},{"cell_type":"code","source":"# Load metadata\ntrain_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntrain_bbox = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv\")\ntest_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\nss = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv\")\n\n# Print dataframe shapes\nprint('train shape:', train_df.shape)\nprint('train bbox shape:', train_bbox.shape)\nprint('test shape:', test_df.shape)\nprint('ss shape:', ss.shape)\nprint('')\n\n# Show first few entries\ntrain_df.head(3)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:21.830289Z","iopub.execute_input":"2022-10-27T20:21:21.830979Z","iopub.status.idle":"2022-10-27T20:21:21.872751Z","shell.execute_reply.started":"2022-10-27T20:21:21.83094Z","shell.execute_reply":"2022-10-27T20:21:21.87174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Debug","metadata":{}},{"cell_type":"code","source":"debug = False\nif len(ss)==3:\n    debug = True\n    \n    # Fix mismatch with test_images folder\n    test_df = pd.DataFrame(columns = ['row_id','StudyInstanceUID','prediction_type'])\n    for i in ['1.2.826.0.1.3680043.22327','1.2.826.0.1.3680043.25399','1.2.826.0.1.3680043.5876']:\n        for j in ['C1','C2','C3','C4','C5','C6','C7','patient_overall']:\n            test_df = test_df.append({'row_id':i+'_'+j,'StudyInstanceUID':i,'prediction_type':j},ignore_index=True)\n    \n    # Sample submission\n    ss = pd.DataFrame(test_df['row_id'])\n    ss['fractured'] = 0.5\n    \n    display(test_df.head(3))","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:21.874395Z","iopub.execute_input":"2022-10-27T20:21:21.87509Z","iopub.status.idle":"2022-10-27T20:21:21.937558Z","shell.execute_reply.started":"2022-10-27T20:21:21.87505Z","shell.execute_reply":"2022-10-27T20:21:21.936628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load volumes","metadata":{}},{"cell_type":"code","source":"# Convert dicom images to 3d tensor\ndef convert_volume(dir_path, out_dir = \"test_volumes\", size = (224, 224, 224)):\n    ls_imgs = glob(os.path.join(dir_path, \"*.dcm\"))\n    ls_imgs = sorted(ls_imgs, key=lambda p: int(os.path.splitext(os.path.basename(p))[0]))\n\n    imgs = []\n    for p_img in ls_imgs:\n        dicom = pydicom.dcmread(p_img)\n        img = apply_voi_lut(dicom.pixel_array, dicom)\n        img = cv2.resize(img, size[:2], interpolation=cv2.INTER_LINEAR)\n        imgs.append(img.tolist())\n    vol = torch.tensor(imgs, dtype=torch.float32)\n\n    vol = (vol - vol.min()) / float(vol.max() - vol.min())\n    vol = interpolate_volume(vol, size).numpy()\n    \n    # https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_adapt_hist_eq_3d.html\n    vol = exposure.equalize_adapthist(vol, kernel_size=np.array([64, 64, 64]), clip_limit=0.01)\n    # vol = exposure.equalize_hist(vol)\n    vol = np.clip(vol * 255, 0, 255).astype(np.uint8)\n    \n    path_pt = os.path.join(out_dir, f\"{os.path.basename(dir_path)}.pt\")\n    torch.save(torch.tensor(vol), path_pt)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:21.939279Z","iopub.execute_input":"2022-10-27T20:21:21.939693Z","iopub.status.idle":"2022-10-27T20:21:21.950434Z","shell.execute_reply.started":"2022-10-27T20:21:21.939654Z","shell.execute_reply":"2022-10-27T20:21:21.949019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make directory\nos.mkdir('/kaggle/working/test_volumes')\n\n# Get paths\nls_dirs = [p for p in glob(os.path.join(\"../input/rsna-2022-cervical-spine-fracture-detection\", \"test_images\", \"*\")) if os.path.isdir(p)]\nprint(f\"volumes: {len(ls_dirs)}\")\n\n# Convert volumes\n_= Parallel(n_jobs=3)(delayed(convert_volume)(p_dir, out_dir='/kaggle/working/test_volumes') for p_dir in tqdm(ls_dirs))","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:21.952393Z","iopub.execute_input":"2022-10-27T20:21:21.953675Z","iopub.status.idle":"2022-10-27T20:21:22.162846Z","shell.execute_reply.started":"2022-10-27T20:21:21.953624Z","shell.execute_reply":"2022-10-27T20:21:22.159536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Torch dataset","metadata":{}},{"cell_type":"code","source":"# Dataset for test set only\nclass RSNADataset(Dataset):\n    # Initialise\n    def __init__(self, subset='test', df_table=test_df):\n        super().__init__()\n        \n        self.subset = subset\n        self.df_table = df_table\n        \n        # Image paths\n        self.volume_dir = '/kaggle/working/test_volumes/'\n        \n    # Get item in position given by index\n    def __getitem__(self, index):\n        \n        # load 3d volume\n        patient = self.df_table.loc[index,'StudyInstanceUID']\n        path = os.path.join(self.volume_dir, f'{patient}.pt')\n        vol = torch.load(path).to(torch.float32)\n        \n        return (vol.unsqueeze(0), patient)\n\n    # Length of dataset\n    def __len__(self):\n        return len(self.df_table)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.164101Z","iopub.status.idle":"2022-10-27T20:21:22.164808Z","shell.execute_reply.started":"2022-10-27T20:21:22.164482Z","shell.execute_reply":"2022-10-27T20:21:22.16451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Test dataset\ntest_table = pd.DataFrame(pd.unique(test_df['StudyInstanceUID']),columns=['StudyInstanceUID'])\ntest_dataset = RSNADataset(subset='test', df_table = test_table)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.167163Z","iopub.status.idle":"2022-10-27T20:21:22.167692Z","shell.execute_reply.started":"2022-10-27T20:21:22.167408Z","shell.execute_reply":"2022-10-27T20:21:22.167437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Torch dataloader","metadata":{}},{"cell_type":"code","source":"# Dataloader\ntest_loader = DataLoader(dataset=test_dataset, batch_size=BATCH_SIZE, shuffle=False)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.169494Z","iopub.status.idle":"2022-10-27T20:21:22.170024Z","shell.execute_reply.started":"2022-10-27T20:21:22.169756Z","shell.execute_reply":"2022-10-27T20:21:22.169779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3. Model","metadata":{}},{"cell_type":"code","source":"# 3D convolutional neural network\nclass Conv3DNet(nn.Module):\n    def __init__(self):\n        super().__init__()\n        \n        # Layers\n        self.conv1 = nn.Conv3d(in_channels=1, out_channels=16, kernel_size=7, stride=1, padding=0)\n        self.pool = nn.MaxPool3d(kernel_size=2, stride=2, padding=0)\n        self.norm1 = nn.BatchNorm3d(num_features=16)\n        self.conv2 = nn.Conv3d(in_channels=16, out_channels=32, kernel_size=3, stride=1, padding=0)\n        self.norm2 = nn.BatchNorm3d(num_features=32)\n        self.conv3 = nn.Conv3d(in_channels=32, out_channels=64, kernel_size=3, stride=1, padding=0)\n        self.norm3 = nn.BatchNorm3d(num_features=64)\n        self.avg = nn.AdaptiveAvgPool3d((7, 1, 1))\n        self.flat = nn.Flatten()\n        self.relu = nn.ReLU()\n        self.lin1 = nn.Linear(in_features=448, out_features=128)\n        self.lin2 = nn.Linear(in_features=128, out_features=8)\n        \n    def forward(self, x):\n        # Conv block 1\n        out = self.conv1(x)\n        out = self.relu(out)\n        out = self.pool(out)\n        out = self.norm1(out)\n        \n        # Conv block 2\n        out = self.conv2(out)\n        out = self.relu(out)\n        out = self.pool(out)\n        out = self.norm2(out)\n        \n        # Conv block 3\n        out = self.conv3(out)\n        out = self.relu(out)\n        out = self.pool(out)\n        out = self.norm3(out)\n        \n        # Average & flatten\n        out = self.avg(out)\n        out = self.flat(out)\n        \n        # Fully connected layer\n        out = self.lin1(out)\n        out = self.relu(out)\n        \n        # Output layer (no sigmoid needed)\n        out = self.lin2(out)\n        \n        return out\n\nmodel = Conv3DNet().to(device)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.171825Z","iopub.status.idle":"2022-10-27T20:21:22.172311Z","shell.execute_reply.started":"2022-10-27T20:21:22.172052Z","shell.execute_reply":"2022-10-27T20:21:22.172076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load model","metadata":{}},{"cell_type":"code","source":"# Load checkpoint\nPATH='../input/rsna-trained-3d-model-weights-pytorch/Conv3DNet.pt'\nif torch.cuda.is_available():\n    checkpoint = torch.load(PATH)\nelse:\n    checkpoint = torch.load(PATH, map_location=torch.device('cpu'))\n\n# Load states\nmodel.load_state_dict(checkpoint['model_state_dict'])\nepoch = checkpoint['epoch']\nloss = checkpoint['loss']\nval_loss = checkpoint['val_loss']\n\n# Evaluation mode\nmodel.eval()\nmodel.to(device)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-27T20:21:22.174025Z","iopub.status.idle":"2022-10-27T20:21:22.17451Z","shell.execute_reply.started":"2022-10-27T20:21:22.174258Z","shell.execute_reply":"2022-10-27T20:21:22.174283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Print final loss and epoch\nprint('Final epoch:', epoch)\nprint('Final loss:', loss)\nprint('Final valid loss:', val_loss)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.176221Z","iopub.status.idle":"2022-10-27T20:21:22.176723Z","shell.execute_reply.started":"2022-10-27T20:21:22.176449Z","shell.execute_reply":"2022-10-27T20:21:22.176472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Inference on test set","metadata":{}},{"cell_type":"code","source":"test_df['fractured']=0.5\nwith torch.no_grad():\n    # Loop over batches\n    for i, (imgs, patient) in enumerate(test_loader):\n        print(f'Iteration {i+1}/{len(test_loader)}')\n        # Send to device\n        imgs = imgs.to(device)\n        \n        # Make predictions\n        preds = model(imgs)\n        \n        # Apply sigmoid\n        sig = nn.Sigmoid()\n        preds = sig(preds)\n        preds = preds.to('cpu')\n        \n        # Save preds\n        test_df.loc[test_df['StudyInstanceUID']==patient[0],'fractured'] = preds.numpy().squeeze()\n        \nprint('Inference complete!')","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.178139Z","iopub.status.idle":"2022-10-27T20:21:22.179083Z","shell.execute_reply.started":"2022-10-27T20:21:22.178822Z","shell.execute_reply":"2022-10-27T20:21:22.178847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Submission","metadata":{}},{"cell_type":"code","source":"submission = test_df[['row_id','fractured']]\n#submission.to_csv('submission2.csv', index=False)\nsubmission.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.187705Z","iopub.status.idle":"2022-10-27T20:21:22.188233Z","shell.execute_reply.started":"2022-10-27T20:21:22.18797Z","shell.execute_reply":"2022-10-27T20:21:22.187997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Adapted from  https://github.com/parthsuresh/3dvhog\n\n#MIT License\n\n#Copyright (c) 2019 Parth Suresh\n\n#Permission is hereby granted, free of charge, to any person obtaining a copy\n#of this software and associated documentation files (the \"Software\"), to deal\n#in the Software without restriction, including without limitation the rights\n#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n#copies of the Software, and to permit persons to whom the Software is\n#furnished to do so, subject to the following conditions:\n\n#The above copyright notice and this permission notice shall be included in all\n#copies or substantial portions of the Software.\n\n#THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\n#SOFTWARE.\n\nimport numpy as np\nimport math\nfrom scipy.ndimage import convolve\nfrom tqdm import tqdm\n\n\ndef hog3d(vox_volume, cell_size, block_size, theta_histogram_bins, phi_histogram_bins, step_size=None):\n    \"\"\"\n    Inputs\n    vox_volume : a \t[x x y x z] numpy array defining voxels with values in the range 0-1\n    cell_size : size of a 3d cell (int)\n    block_size : size of a 3d block defined in cells\n    theta_histogram_bins : number of bins to break the angles in the xy plane - 180 degrees\n    phi_histogram_bins : number of bins to break the angles in the xz plane - 360 degrees\n    step_size : OPTIONAL integer defining the number of cells the blocks should overlap by.\n\t\"\"\"\n\n    if step_size is None:\n        step_size = block_size\n\n    c = cell_size\n    b = block_size\n\n    sx, sy, sz = vox_volume.shape\n\n    num_x_cells = math.floor(sx / cell_size)\n    num_y_cells = math.floor(sy / cell_size)\n    num_z_cells = math.floor(sz / cell_size)\n\n    # Get cell positions\n    x_cell_positions = np.array(list(range(0, (num_x_cells * cell_size), cell_size)))\n    y_cell_positions = np.array(list(range(0, (num_y_cells * cell_size), cell_size)))\n    z_cell_positions = np.array(list(range(0, (num_z_cells * cell_size), cell_size)))\n\n    # Get block positions\n    x_block_positions = (x_cell_positions[0: num_x_cells: block_size])\n    y_block_positions = (y_cell_positions[0: num_y_cells: block_size])\n    z_block_positions = (z_cell_positions[0: num_z_cells: block_size])\n\n    # Check if last block in each dimension has enough voxels to be a full block. If not, discard it.\n    if x_block_positions[-1] > ((sx + 1) - (cell_size * block_size)):\n        x_block_positions = x_block_positions[:-2]\n    if y_block_positions[-1] > ((sy + 1) - (cell_size * block_size)):\n        y_block_positions = y_block_positions[:-2]\n    if z_block_positions[-1] > ((sz + 1) - (cell_size * block_size)):\n        z_block_positions = z_block_positions[:-2]\n\n    # Number of blocks\n    num_x_blocks = len(x_block_positions)\n    num_y_blocks = len(y_block_positions)\n    num_z_blocks = len(z_block_positions)\n\n    # Create 3D gradient vectors\n    # X filter and vector\n    x_filter = np.zeros((3, 3, 3))\n    x_filter[0, 1, 1], x_filter[2, 1, 1] = 1, -1\n    x_vector = convolve(vox_volume, x_filter, mode='constant', cval=0)\n\n    # Y filter and vector\n    y_filter = np.zeros((3, 3, 3))\n    y_filter[1, 0, 0], y_filter[1, 2, 0] = 1, -1\n    y_vector = convolve(vox_volume, y_filter, mode='constant', cval=0)\n\n    # Z filter and vector\n    z_filter = np.zeros((3, 3, 3))\n    z_filter[1, 1, 0], z_filter[1, 1, 2] = 1, -1\n    z_vector = convolve(vox_volume, z_filter, mode='constant', cval=0)\n\n    magnitudes = np.zeros([sx, sy, sz])\n    for i in range(sx):\n        for j in range(sy):\n            for k in range(sz):\n                magnitudes[i, j, k] = (x_vector[i, j, k] ** 2 + y_vector[i, j, k] ** 2 + z_vector[i, j, k] ** 2) ** (\n                    0.5)\n\n    # Voxel Weights\n    kernel_size = 3\n    voxel_filter = np.full((kernel_size, kernel_size, kernel_size), 1 / (kernel_size * kernel_size * kernel_size))\n    weights = convolve(vox_volume, voxel_filter, mode='constant', cval=0)\n    weights = weights + 1\n\n    # Gradient vector\n    grad_vector = np.zeros((sx, sy, sz, 3))\n    for i in range(sx):\n        for j in range(sy):\n            for k in range(sz):\n                grad_vector[i, j, k, 0] = x_vector[i, j, k]\n                grad_vector[i, j, k, 1] = y_vector[i, j, k]\n                grad_vector[i, j, k, 2] = z_vector[i, j, k]\n\n    theta = np.zeros((sx, sy, sz))\n    phi = np.zeros((sx, sy, sz))\n    for i in range(sx):\n        for j in range(sy):\n            for k in range(sz):\n                theta[i, j, k] = math.acos(grad_vector[i, j, k, 2])\n                phi[i, j, k] = math.atan2(grad_vector[i, j, k, 1], grad_vector[i, j, k, 0])\n                phi[i, j, k] += math.pi\n\n    # Binning\n    b_size_voxels = int(c * b)\n    t_hist_bins = math.pi / theta_histogram_bins\n    p_hist_bins = (2 * math.pi) / phi_histogram_bins\n\n    block_inds = np.zeros((num_x_blocks * num_y_blocks * num_z_blocks, 3))\n    i = 0\n    for z_block in range(num_z_blocks):\n        for y_block in range(num_y_blocks):\n            for x_block in range(num_x_blocks):\n                block_inds[i] = np.array(\n                    [x_block_positions[x_block], y_block_positions[y_block], z_block_positions[z_block]])\n                i += 1\n\n    num_blocks = len(block_inds)\n    error_count = 0\n    features = []\n    for i in range(num_blocks):\n        full_empty = vox_volume[int(block_inds[i, 0]):int(block_inds[i, 0] + b_size_voxels),\n                     int(block_inds[i, 1]):int(block_inds[i, 1] + b_size_voxels),\n                     int(block_inds[i, 2]):int(block_inds[i, 2] + b_size_voxels)]\n\n        if np.sum(full_empty) != 0 and np.sum(full_empty) != full_empty.size:\n            feature = np.zeros((b, b, b, theta_histogram_bins, phi_histogram_bins))\n            t_weights = weights[int(block_inds[i, 0]):int(block_inds[i, 0] + b_size_voxels),\n                        int(block_inds[i, 1]):int(block_inds[i, 1] + b_size_voxels),\n                        int(block_inds[i, 2]):int(block_inds[i, 2] + b_size_voxels)]\n            t_magnitudes = magnitudes[int(block_inds[i, 0]):int(block_inds[i, 0] + b_size_voxels),\n                           int(block_inds[i, 1]):int(block_inds[i, 1] + b_size_voxels),\n                           int(block_inds[i, 2]):int(block_inds[i, 2] + b_size_voxels)]\n            t_theta = theta[int(block_inds[i, 0]):int(block_inds[i, 0] + b_size_voxels),\n                      int(block_inds[i, 1]):int(block_inds[i, 1] + b_size_voxels),\n                      int(block_inds[i, 2]):int(block_inds[i, 2] + b_size_voxels)]\n            t_phi = phi[int(block_inds[i, 0]):int(block_inds[i, 0] + b_size_voxels),\n                    int(block_inds[i, 1]):int(block_inds[i, 1] + b_size_voxels),\n                    int(block_inds[i, 2]):int(block_inds[i, 2] + b_size_voxels)]\n\n            for l in range(b_size_voxels):\n                for m in range(b_size_voxels):\n                    for n in range(b_size_voxels):\n                        cell_pos_x = math.ceil(l / c) - 1\n                        cell_pos_y = math.ceil(m / c) - 1\n                        cell_pos_z = math.ceil(n / c) - 1\n\n                        hist_pos_theta = math.ceil(t_theta[l, m, n] / t_hist_bins) - 1\n                        hist_pos_phi = math.ceil(t_phi[l, m, n] / p_hist_bins) - 1\n\n                        if phi_histogram_bins >= hist_pos_phi > 0 and theta_histogram_bins >= hist_pos_theta > 0:\n                            feature[cell_pos_x, cell_pos_y, cell_pos_z, hist_pos_theta, hist_pos_phi] += (\n                                    t_magnitudes[l, m, n] * t_weights[l, m, n])\n                        else:\n                            error_count += 1\n\n            feature = np.reshape(feature, ((b * b * b), theta_histogram_bins, phi_histogram_bins))\n            l2 = np.linalg.norm(feature)\n            if l2 != 0:\n                norm_feature = feature / l2\n            else:\n                norm_feature = feature\n            norm_feature = np.reshape(norm_feature, ((b * b * b), (theta_histogram_bins * phi_histogram_bins)))\n\n            features.append(norm_feature)\n\n    features = np.array(features)\n\n    return features","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.190257Z","iopub.status.idle":"2022-10-27T20:21:22.190773Z","shell.execute_reply.started":"2022-10-27T20:21:22.190491Z","shell.execute_reply":"2022-10-27T20:21:22.190515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -qU ../input/for-pydicom/python_gdcm-3.0.14-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl ../input/for-pydicom/pylibjpeg-1.4.0-py3-none-any.whl --find-links frozen_packages --no-index","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.192552Z","iopub.status.idle":"2022-10-27T20:21:22.193486Z","shell.execute_reply.started":"2022-10-27T20:21:22.193235Z","shell.execute_reply":"2022-10-27T20:21:22.19326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -q kaggle_vol3d_classify -f ../input/cervical-spine-fracture-detection-npz-3d-volumes/frozen_packages --no-index","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.195112Z","iopub.status.idle":"2022-10-27T20:21:22.195607Z","shell.execute_reply.started":"2022-10-27T20:21:22.195343Z","shell.execute_reply":"2022-10-27T20:21:22.195367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport matplotlib.patches as patches\nimport seaborn as sns\nsns.set(style='darkgrid', font_scale=1.6)\nimport cv2\nimport os\nfrom os import listdir\nimport re\nimport gc\nimport random\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom tqdm.auto import tqdm\nfrom pprint import pprint\nfrom time import time\nimport itertools\nfrom skimage import measure\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport nibabel as nib\nfrom glob import glob\nimport warnings\n#warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n#warnings.filterwarnings(\"ignore\", category=UserWarning)\n#warnings.filterwarnings(\"ignore\", category=FutureWarning)\nimport zipfile\nfrom scipy import ndimage\nfrom sklearn.model_selection import train_test_split\nfrom joblib import Parallel, delayed\nfrom PIL import Image\nfrom dipy.denoise.nlmeans import nlmeans\nfrom dipy.denoise.noise_estimate import estimate_sigma\nfrom kaggle_volclassif.utils import interpolate_volume\nfrom skimage import exposure\n\n# Pytorch\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.optim.lr_scheduler as lr_scheduler\nfrom torch.utils.data import Dataset, DataLoader\nimport torchvision.transforms as transforms\nimport torch.nn.functional as F\nimport kornia\nimport kornia.augmentation as augmentation","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.197348Z","iopub.status.idle":"2022-10-27T20:21:22.197844Z","shell.execute_reply.started":"2022-10-27T20:21:22.197579Z","shell.execute_reply":"2022-10-27T20:21:22.197618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load metadata\ntrain_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntrain_bbox = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv\")\ntest_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\nss = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv\")\n\n# Print dataframe shapes\nprint('train shape:', train_df.shape)\nprint('train bbox shape:', train_bbox.shape)\nprint('test shape:', test_df.shape)\nprint('ss shape:', ss.shape)\nprint('')\n\n# Show first few entries\ntrain_df.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.199555Z","iopub.status.idle":"2022-10-27T20:21:22.200064Z","shell.execute_reply.started":"2022-10-27T20:21:22.19981Z","shell.execute_reply":"2022-10-27T20:21:22.199833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/344862\nbad_scans = ['1.2.826.0.1.3680043.20574','1.2.826.0.1.3680043.29952']\n\nfor uid in bad_scans:\n    train_df.drop(train_df[train_df['StudyInstanceUID']==uid].index, axis=0, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.201807Z","iopub.status.idle":"2022-10-27T20:21:22.202276Z","shell.execute_reply.started":"2022-10-27T20:21:22.202031Z","shell.execute_reply":"2022-10-27T20:21:22.202055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"debug = False\nif len(ss)==3:\n    debug = True\n    \n    # Fix mismatch with test_images folder\n    test_df = pd.DataFrame(columns = ['row_id','StudyInstanceUID','prediction_type'])\n    for i in ['1.2.826.0.1.3680043.22327','1.2.826.0.1.3680043.25399','1.2.826.0.1.3680043.5876']:\n        for j in ['C1','C2','C3','C4','C5','C6','C7','patient_overall']:\n            test_df = test_df.append({'row_id':i+'_'+j,'StudyInstanceUID':i,'prediction_type':j},ignore_index=True)\n    \n    # Sample submission\n    ss = pd.DataFrame(test_df['row_id'])\n    ss['fractured'] = 0.5\n    \n    display(test_df.head(3))","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.20382Z","iopub.status.idle":"2022-10-27T20:21:22.20448Z","shell.execute_reply.started":"2022-10-27T20:21:22.204207Z","shell.execute_reply":"2022-10-27T20:21:22.204233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"AUGMENTATIONS = True\n# Data augmentations (https://kornia.readthedocs.io/en/latest/augmentation.module.html#geometric)\nif AUGMENTATIONS:\n    augs = transforms.Compose([\n        augmentation.RandomRotation3D((0,0,30), resample='bilinear', p=0.5, same_on_batch=False, keepdim=True),\n        #augmentation.RandomHorizontalFlip3D(same_on_batch=False, p=0.5, keepdim=True),\n        ])\nelse:\n    augs=None","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.206427Z","iopub.status.idle":"2022-10-27T20:21:22.206936Z","shell.execute_reply.started":"2022-10-27T20:21:22.206681Z","shell.execute_reply":"2022-10-27T20:21:22.206705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dataset for train/valid sets only\nclass RSNADataset(Dataset):\n    # Initialise\n    def __init__(self, subset='train', df_table = train_df, transform=None):\n        super().__init__()\n        \n        self.subset = subset\n        self.df_table = df_table.reset_index(drop=True)\n        self.transform = transform\n        self.targets = ['C1','C2','C3','C4','C5','C6','C7','patient_overall']\n        \n        # Identify files in each of the two datasets\n        fh_paths = glob(os.path.join('../input/rsna-3d-train-tensors-first-half/train_volumes', \"*.pt\"))\n        sh_paths = glob(os.path.join('../input/rsna-3d-train-tensors-second-half/train_volumes', \"*.pt\"))\n        \n        fh_list = []\n        sh_list = []\n        for i in fh_paths:\n            fh_list.append(i.split('/')[-1][:-3])\n        \n        for i in sh_paths:\n            sh_list.append(i.split('/')[-1][:-3])\n        \n        self.df_table_fh = self.df_table[self.df_table['StudyInstanceUID'].isin(fh_list)]\n        self.df_table_sh = self.df_table[self.df_table['StudyInstanceUID'].isin(sh_list)]\n        \n        # Image paths\n        self.volume_dir1 = '../input/rsna-3d-train-tensors-first-half/train_volumes'  # <=1000 patient\n        self.volume_dir2 = '../input/rsna-3d-train-tensors-second-half/train_volumes' # >1000 patient\n\n        # Populate labels\n        self.labels = self.df_table[self.targets].values\n        \n    # Get item in position given by index\n    def __getitem__(self, index):\n        if index in self.df_table_fh.index:\n            patient = self.df_table_fh[self.df_table_fh.index==index]['StudyInstanceUID'].iloc[0]\n            path = os.path.join(self.volume_dir1, f\"{patient}.pt\")\n            vol = torch.load(path).to(torch.float32)\n        else:\n            patient = self.df_table_sh[self.df_table_sh.index==index]['StudyInstanceUID'].iloc[0]\n            path = os.path.join(self.volume_dir2, f\"{patient}.pt\")\n            vol = torch.load(path).to(torch.float32)\n        \n        # Data augmentations\n        if self.transform:\n            vol = self.transform(vol)\n        \n        return vol.unsqueeze(0), self.labels[index]\n\n    # Length of dataset\n    def __len__(self):\n        return len(self.df_table['StudyInstanceUID'])","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.208806Z","iopub.status.idle":"2022-10-27T20:21:22.209291Z","shell.execute_reply.started":"2022-10-27T20:21:22.209035Z","shell.execute_reply":"2022-10-27T20:21:22.20906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_table, valid_table = train_test_split(train_df, train_size=0.85, test_size=0.15, random_state=0)\ntrain_dataset = RSNADataset(subset='train', df_table = train_table, transform=augs)\nvalid_dataset = RSNADataset(subset='valid', df_table = valid_table)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.211Z","iopub.status.idle":"2022-10-27T20:21:22.211499Z","shell.execute_reply.started":"2022-10-27T20:21:22.21124Z","shell.execute_reply":"2022-10-27T20:21:22.211264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.ndimage import zoom\nfrom scipy import ndimage, misc\n\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.213252Z","iopub.status.idle":"2022-10-27T20:21:22.214361Z","shell.execute_reply.started":"2022-10-27T20:21:22.21408Z","shell.execute_reply":"2022-10-27T20:21:22.214107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cell_size = 15\nblock_size = 2\ntheta_histogram_bins = 4\nphi_histogram_bins = 4\n\n#grad_vec = hog3d(x2, cell_size, block_size, theta_histogram_bins, phi_histogram_bins)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.216069Z","iopub.status.idle":"2022-10-27T20:21:22.216602Z","shell.execute_reply.started":"2022-10-27T20:21:22.216331Z","shell.execute_reply":"2022-10-27T20:21:22.216355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compressed_data = []\n#train_size = len(train_dataset)\ntrain_size = 200\nfor i in range(train_size):\n    x = train_dataset[i][0][0].cpu().detach().numpy()/255\n    x2 = np.array(zoom(x,0.4))\n    #print(x2.max(),x2.min())\n    x2 = np.abs(x2)/x2.max()\n    #print(x2.max(), x2.min())\n    grad_vec = hog3d(x2, cell_size, block_size, theta_histogram_bins, phi_histogram_bins)\n    compressed_data.append(grad_vec.flatten())","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.218392Z","iopub.status.idle":"2022-10-27T20:21:22.219028Z","shell.execute_reply.started":"2022-10-27T20:21:22.21877Z","shell.execute_reply":"2022-10-27T20:21:22.218797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compressed_valid = []\n#train_size = len(valid_dataset)\nvalid_size = 2\nfor i in range(valid_size):\n    x = valid_dataset[i][0][0].cpu().detach().numpy()/255\n    x2 = np.array(zoom(x,0.4))\n    #print(x2.max(),x2.min())\n    x2 = np.abs(x2)/x2.max()\n    #print(x2.max(), x2.min())\n    grad_vec = hog3d(x2, cell_size, block_size, theta_histogram_bins, phi_histogram_bins)\n    compressed_valid.append(grad_vec.flatten())","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.220528Z","iopub.status.idle":"2022-10-27T20:21:22.221575Z","shell.execute_reply.started":"2022-10-27T20:21:22.221302Z","shell.execute_reply":"2022-10-27T20:21:22.221332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nlda_classifiers = []\n\nfor j in range(8):\n    y_train = []\n    for i in range(train_size):\n        y_train.append(train_dataset[i][1][j])\n    \n    y_valid = []\n    for i in range(valid_size):\n        y_valid.append(valid_dataset[i][1][j])\n    lda_clf = LinearDiscriminantAnalysis()\n    lda_clf.fit(compressed_data, y_train)\n    lda_classifiers.append(lda_clf)\n    print(lda_clf.score(compressed_data, y_train),lda_clf.score(compressed_valid, y_valid))","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.223461Z","iopub.status.idle":"2022-10-27T20:21:22.223959Z","shell.execute_reply.started":"2022-10-27T20:21:22.223708Z","shell.execute_reply":"2022-10-27T20:21:22.223732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compressed_test = []\ntest_size = len(test_dataset)\nfor i in range(test_size):\n    x = test_dataset[i][0][0].cpu().detach().numpy()/255\n    x2 = np.array(zoom(x,0.4))\n    #print(x2.max(),x2.min())\n    x2 = np.abs(x2)/x2.max()\n    #print(x2.max(), x2.min())\n    grad_vec = hog3d(x2, cell_size, block_size, theta_histogram_bins, phi_histogram_bins)\n    compressed_test.append(grad_vec.flatten())","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.225381Z","iopub.status.idle":"2022-10-27T20:21:22.226429Z","shell.execute_reply.started":"2022-10-27T20:21:22.226157Z","shell.execute_reply":"2022-10-27T20:21:22.226187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lda_predicts = []\nfor i in range(test_size):\n    y_test = []\n        \n    for j in range(8):\n        y_test.append(lda_classifiers[j].predict([compressed_test[i]])[0])\n    \n    lda_predicts.append(y_test)\n\nlda_predicts = np.array(lda_predicts).flatten()\n#print(lda_predicts)","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.22805Z","iopub.status.idle":"2022-10-27T20:21:22.22854Z","shell.execute_reply.started":"2022-10-27T20:21:22.22829Z","shell.execute_reply":"2022-10-27T20:21:22.228313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ensembling Solution","metadata":{}},{"cell_type":"code","source":"#sub = df_sub[['row_id', 'fractured']].copy()\nsub = submission.copy()","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.23022Z","iopub.status.idle":"2022-10-27T20:21:22.230715Z","shell.execute_reply.started":"2022-10-27T20:21:22.230447Z","shell.execute_reply":"2022-10-27T20:21:22.230471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(df_sub.shape[0]):\n    idx = df_sub['row_id'][i]\n    val1 = df_sub['fractured'][df_sub['row_id'] == str(idx)].values[0]\n    val2 = submission['fractured'][submission['row_id'] == str(idx)].values[0]\n\n    sub['fractured'][sub['row_id'] == str(idx)] = val1*0.459 + val2*0.341\n\nsub['fractured']+=lda_predicts*0.200\nsub.to_csv(\"submission.csv\", index=False)\nsub","metadata":{"execution":{"iopub.status.busy":"2022-10-27T20:21:22.232104Z","iopub.status.idle":"2022-10-27T20:21:22.23291Z","shell.execute_reply.started":"2022-10-27T20:21:22.232638Z","shell.execute_reply":"2022-10-27T20:21:22.232663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-warning\" style=\"text-align:center; font-size:28px;\">\n    Thanks for reading 🤗\n</div>","metadata":{}}]}