{"cells":[{"metadata":{},"cell_type":"markdown","source":"## setup"},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install -q ../input/monai030/monai-0.3.0-202010042353-py3-none-any.whl","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import sys\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport PIL.Image\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm as tqdm\nfrom sklearn.metrics import cohen_kappa_score, confusion_matrix\nfrom sklearn.model_selection import train_test_split, StratifiedKFold\nimport random\n\nimport torch\nfrom torch.utils.data import TensorDataset, DataLoader,Dataset\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torchvision\nimport torchvision.transforms as transforms\nimport torch.optim as optim\nfrom torch.optim import lr_scheduler\nfrom torch.utils.data.sampler import SubsetRandomSampler, RandomSampler, SequentialSampler\nfrom torch.optim.lr_scheduler import StepLR, ReduceLROnPlateau, CosineAnnealingLR\nfrom sklearn.metrics import roc_auc_score\nimport albumentations\n\nimport pydicom\nimport os, os.path as osp\n\nimport monai\nfrom monai.transforms import LoadNifti, Randomizable, apply_transform\nfrom monai.transforms import AddChannel, Compose, RandRotate90, Resize, ScaleIntensity, ToTensor\nfrom monai.utils import get_seed\n\nfrom scipy.ndimage.interpolation import zoom\nfrom tqdm import tqdm\nfrom glob import glob\n\ndevice = torch.device('cuda')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_size = 160\nout_dim = 9\nbatch_size = 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!cp ../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv submission.csv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train= pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/train.csv')\nsub =  pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv')\ntest = pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/test.csv')\n\nDEBUG = (test.shape[0]==146853)\nDEBUG","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_study = test.drop_duplicates('StudyInstanceUID')[['StudyInstanceUID','SeriesInstanceUID']]\n\n# save time during commit\nif DEBUG:\n    test_study = test_study.head(25)\ntest_study.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_dicom_array(f):\n    dicom_files = glob(osp.join(f, '*.dcm'))\n    dicoms = [pydicom.dcmread(d) for d in dicom_files]\n    M = float(dicoms[0].RescaleSlope)\n    B = float(dicoms[0].RescaleIntercept)\n    # Assume all images are axial\n    z_pos = [float(d.ImagePositionPatient[-1]) for d in dicoms]\n    dicoms = np.asarray([d.pixel_array for d in dicoms])\n    dicoms = dicoms[np.argsort(z_pos)]\n    dicoms = dicoms * M\n    dicoms = dicoms + B\n    return dicoms, np.asarray(dicom_files)[np.argsort(z_pos)]\n\ndef window(img, WL=50, WW=350):\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    return X\n\ndef read_dicom(dcm_path, image_size=256):\n    image, files = load_dicom_array(dcm_path)\n    # Windows from https://pubs.rsna.org/doi/pdf/10.1148/rg.245045008\n    image_lung = np.expand_dims(window(image, WL=-600, WW=1500), axis=3)\n    image_mediastinal = np.expand_dims(window(image, WL=40, WW=400), axis=3)\n    image_pe_specific = np.expand_dims(window(image, WL=100, WW=700), axis=3)\n    image = np.concatenate([image_mediastinal, image_pe_specific, image_lung], axis=3)\n    rat = image_size / np.max(image.shape[1:])\n    image = zoom(image, [1.,rat,rat,1.], prefilter=False, order=1)\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class RSNADataset3D(Dataset, Randomizable):\n    def __init__(self, csv, mode, transform=None):\n        self.csv = csv.reset_index()\n        self.mode = mode\n        self.transform = transform\n    def __len__(self):\n        return self.csv.shape[0]\n    def randomize(self) -> None:\n        MAX_SEED = np.iinfo(np.uint32).max + 1\n        self._seed = self.R.randint(MAX_SEED, dtype=\"uint32\")    \n    def __getitem__(self, index):\n        self.randomize()\n        row = self.csv.iloc[index]\n        try:\n            img = read_dicom(os.path.join('../input/rsna-str-pulmonary-embolism-detection/test', row.StudyInstanceUID, row.SeriesInstanceUID))\n        except:\n            img = np.zeros((144, 256, 256, 3),dtype=np.uint8)\n            \n        # (144, 256, 256, 3)  Z, H, W, ch\n        img = img[:,:,:,::-1].transpose(3,1,2,0) # -> ch, H, W, Z\n        if self.transform is not None:\n            if isinstance(self.transform, Randomizable):\n                self.transform.set_random_state(seed=self._seed)\n            img = apply_transform(self.transform, img)\n        if self.mode == 'test':\n            return img   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"val_transforms = Compose([ScaleIntensity(), Resize((image_size, image_size, image_size)), ToTensor()])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if DEBUG:\n    dataset_show = RSNADataset3D(test_study.head(5), 'test', transform=val_transforms)\n    from pylab import rcParams\n    rcParams['figure.figsize'] = 12,5\n    for i in range(3):\n        f, axarr = plt.subplots(1,3)\n        img = dataset_show[i]\n        print(img.shape)\n        for j in range(3):            \n            axarr[j].imshow(img.numpy().transpose(1,2,3,0).mean(axis=j))        \n            axarr[j].set_title(i)\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## load model trained locally"},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls -lrt ../input/monai3d-160-3ch-1e-5-20ep-aug/monai3d_160_3ch_1e-5_20ep_aug_best_fold0.pth","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_files = [f'../input/monai3d-160-3ch-1e-5-20ep-aug/monai3d_160_3ch_1e-5_20ep_aug_best_fold{i}.pth' for i in range(1)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_model(model_file):\n    model = monai.networks.nets.densenet.densenet121(spatial_dims=3, in_channels=3, out_channels=out_dim).to(device)\n\n    try:  # single GPU model_file\n        model.load_state_dict(torch.load(model_file), strict=True)\n    except:  # multi GPU model_file\n        state_dict = torch.load(model_file)\n        state_dict = {k[7:] if k.startswith('module.') else k: state_dict[k] for k in state_dict.keys()}\n        model.load_state_dict(state_dict, strict=True)\n\n    model.eval()    \n    print()\n    return model\n\nmodels = [load_model(model) for model in model_files]\nlen(models)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset_test = RSNADataset3D(test_study, 'test', transform=val_transforms)\ntest_loader = DataLoader(dataset_test, batch_size=8, num_workers=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"LOGITS = []\nwith torch.no_grad():\n    for data in tqdm(test_loader):\n        data = data.to(device)\n        logits = torch.zeros((data.shape[0], out_dim)).to(device)\n        for model in models:\n            l = model(data)\n            logits += l\n        logits /= len(models)\n        LOGITS.append(logits.detach().cpu())\nPROBS = torch.sigmoid(torch.cat(LOGITS)).numpy().squeeze()    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"PROBS.max(), PROBS.min()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## weighted mean prediction per slice location for single images\n\nhttps://www.kaggle.com/osciiart/baseline-with-no-image/"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get dicom paths\ndf_test = test\ndf_test['path'] = (\"../input/rsna-str-pulmonary-embolism-detection/test/\" \n                   + df_test['StudyInstanceUID'].values + \"/\"\n                   + df_test['SeriesInstanceUID'].values + \"/\"\n                   + df_test['SOPInstanceUID'].values + \".dcm\"\n                  )\nprint(df_test['path'][0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# extract exam (study) level data\ncol_index = 'SOPInstanceUID'\ncol_groupby = 'StudyInstanceUID'\ndf_test_study = df_test[df_test[col_groupby].duplicated()==False].reset_index(drop=True)\ndf_tmp = df_test.groupby(col_groupby)[col_index].agg(len).reset_index()\ndf_tmp.columns = [col_groupby, 'num_images']\ndf_test_study = pd.merge(df_test_study, df_tmp, on=col_groupby, how='left')\ndf_test = pd.merge(df_test, df_test_study[[col_groupby, 'num_images']], on=col_groupby, how='left')\nprint(df_test.shape)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get series index of image\ndef task(i):\n    if (i+1)%10000==0:\n        print(\"{}/{} {:.1f}\".format(i+1, len(df_test), time.time()-starttime))\n    path = df_test['path'][i]\n    tmp_dcm = pydicom.dcmread(path)\n    return tmp_dcm.ImagePositionPatient[-1]\n\nimport time\nimport multiprocessing\nfrom concurrent.futures import ProcessPoolExecutor\n\nstarttime = time.time()\nexecutor = ProcessPoolExecutor(max_workers=multiprocessing.cpu_count())\n# futures = [executor.submit(task, i) for i in range(10000)]\nfutures = [executor.submit(task, i) for i in range(len(df_test))]\nresult_list = []\nfor i in range(len(futures)):\n    result_list.append(futures[i].result())\ndf_test['z_pos'] = result_list\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate slice location\ndf_tmp = []\nfor i in range(len(df_test_study)):\n    if (i+1)%100==0: print(\"{}/{}\".format(i+1, len(df_test_study)))\n    study = df_test_study[col_groupby][i]\n    df_study = df_test[df_test[col_groupby]==study].sort_values('z_pos').reset_index(drop=True)\n    df_study['series_index'] = np.arange(len(df_study))\n    df_tmp.append(df_study[[col_index, 'series_index']])\ndf_tmp = pd.concat(df_tmp)\n\ndf_test = pd.merge(df_test, df_tmp, on=col_index, how='left')\n# df_test = pd.merge(df_test, df_test_study[[col_groupby, 'num_images']], on=col_groupby, how='left')\ndf_test['slice_location'] = df_test['series_index'] / (df_test['num_images'] - 1)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get weighted mean prediction per slice location\nq_weighted_means = np.array([0.00326324, 0.05970682, 0.32645303, 0.67452216, 0.71344817, 0.4734337, 0.0740926, 0.00369781])\n\ndf_test = df_test.copy()\nbins = 8\ndf_test['bins'] = bins-1\nfor i in range(bins):\n    df_test['bins'][(df_test['slice_location']>=(i/bins)) & (df_test['slice_location']<((i+1)/bins))] = i\ndf_test['q_weighted_means'] = df_test['bins'].apply(lambda x: q_weighted_means[x])\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## merge the study-level predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub =  pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv')\nsub = pd.merge(sub[['id']], df_test[['SOPInstanceUID','q_weighted_means']].rename(columns={'SOPInstanceUID':'id','q_weighted_means':'label'}), on='id', how='left')\nsub.fillna(0.5,inplace=True)\nsub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"target_cols = [\n        'negative_exam_for_pe', # exam level\n        'rv_lv_ratio_gte_1', # exam level\n        'rv_lv_ratio_lt_1', # exam level\n        'leftsided_pe', # exam level\n        'chronic_pe', # exam level\n        'rightsided_pe', # exam level\n        'acute_and_chronic_pe', # exam level\n        'central_pe', # exam level\n        'indeterminate' # exam level\n    ]\nsub.set_index('id',inplace=True)\nfor i,target in enumerate(target_cols):\n    sub.loc[[x+'_'+target for x in test_study.StudyInstanceUID.values],'label']= PROBS[:,i]\nsub.reset_index(inplace=True)\nsub.to_csv('submission.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}