{"cells":[{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# library import\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os, glob, pickle, gc, copy, sys, multiprocessing\nfrom joblib import Parallel, delayed\n\nimport warnings\nimport cv2, pydicom\nwarnings.filterwarnings('ignore')\npd.set_option('display.max_columns', 100)\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.parallel\nimport torch.backends.cudnn as cudnn\nimport torch.utils.data\nimport torchvision.transforms as transforms\nimport torchvision.datasets as datasets\nfrom torch.utils.data.dataset import Dataset\nfrom torch.utils.data import DataLoader\nimport torch.nn.functional as F\nfrom torch.optim.lr_scheduler import _LRScheduler\n\nsys.path.append('../input/timm-efficientnet/pytorch-image-models-master/')\nimport timm\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\ntorch.backends.cudnn.benchmark = True\n\n# params\ncol_index = 'SOPInstanceUID'\ncol_groupby = 'StudyInstanceUID'\nBATCH_SIZE = 64\ncol_targets = [\n    'negative_exam_for_pe',\n    'indeterminate',\n    'chronic_pe',\n    'acute_and_chronic_pe',\n    'central_pe',\n    'leftsided_pe',\n    'rightsided_pe',\n    'rv_lv_ratio_gte_1',\n    'rv_lv_ratio_lt_1',\n    'pe_present_on_image',\n]\ncol_targets_display = [\n    'Negative Exam for PE',\n    'Indeterminate',\n    'Chronic PE',\n    'Acute and Chronic PE',\n    'Central PE',\n    'Left-sided PE',\n    'Right-sided PE',\n    'RV/LV Ratio: > or = 1',\n    'RV/LV Ratio: < 1',\n    'PE Present on Image',\n]\nnum_features1_1 = 1280\n\n# data loading\ndf_valid = pd.read_csv(\"../input/rsna2020-demo/df_valid.csv\")\ndf_valid_exam = pd.read_csv(\"../input/rsna2020-demo/df_valid_exam.csv\")\ndf_valid_exam['3class'] = 1\ndf_valid_exam['3class'][df_valid_exam['negative_exam_for_pe']==1] = 0\ndf_valid_exam['3class'][df_valid_exam['indeterminate']==1] = 2\n\n# model loading\ndef get_images(idx):\n    exam = df_valid_exam[col_groupby][idx]\n    start_index = df_valid_exam['start_index'][idx]\n    end_index = start_index + df_valid_exam['num_series'][idx]\n    df_tmp = df_valid.iloc[start_index:end_index].reset_index(drop=True)\n    RescaleSlope = df_valid_exam['RescaleSlope'][idx]\n    RescaleIntercept = df_valid_exam['RescaleIntercept'][idx]\n    # load dicoms\n    images_exam = []\n    z_pos = []\n    for i in range(len(df_tmp)):\n        tmp_path = df_tmp['path'][i]\n#         print(tmp_path)\n        tmp_dcm = pydicom.dcmread(tmp_path)\n        tmp_npy = np.asarray(tmp_dcm.pixel_array)\n        images_exam.append(tmp_npy)\n\n    # process images\n    images_exam = np.array(images_exam)\n    images_exam_processed = (images_exam.astype(np.float32) * RescaleSlope + RescaleIntercept)/1000\n    images_exam_processed = images_exam_processed.reshape([-1, 1, 512, 512]).astype(np.float16)\n\n    return images_exam_processed\n\nclass nnWindow(nn.Module):\n    def __init__(self):\n        super(nnWindow, self).__init__()\n        wso = np.array(((40,80),(80,200),(40,400)))/1000\n        conv_ = nn.Conv2d(1,3, kernel_size=(1, 1))\n        conv_.weight.data.copy_(torch.tensor([[[[1./wso[0][1]]]],[[[1./wso[1][1]]]],[[[1./wso[2][1]]]]]))\n        conv_.bias.data.copy_(torch.tensor([0.5 - wso[0][0]/wso[0][1],\n                                            0.5 - wso[1][0]/wso[1][1],\n                                            0.5 -wso[2][0]/wso[2][1]]))\n        self.window = nn.Sequential(\n            conv_,\n            nn.Sigmoid(),\n            nn.InstanceNorm2d(3)\n        )\n    def forward(self, input1):\n        return self.window(input1)\n        \n        \nclass CNN_2D(nn.Module):\n    def __init__(self, num_classes=10, base_model='tf_efficientnet_b0_ns'):\n        super(CNN_2D, self).__init__()\n\n        self.num_classes = num_classes\n        self.mode = 'train'\n        self.window = nnWindow()\n#         self.base_model = pretrainedmodels.__dict__['resnet18'](num_classes=1000, pretrained='imagenet')\n        self.base_model = timm.create_model(base_model, pretrained=False, num_classes=10).to(device, non_blocking=True)\n        self.avgpool = nn.AdaptiveAvgPool2d(1)\n#         self.last_linear = nn.Linear(512, num_classes+1)\n        self.last_linear = nn.Linear(self.base_model.num_features, num_classes)\n\n    def forward(self, input1):\n        bs, ch, h, w = input1.size()\n        x = self.window(input1)\n        x = self.base_model.forward_features(x) #; print('layer conv1 ',x.size()) # [8, 64, 112, 112]\n        feature = self.avgpool(x).view(bs, -1)\n        y = self.last_linear(feature)\n\n        return y\n\n    def feature(self, input1):\n        bs, ch, h, w = input1.size()\n        x = self.window(input1)\n        x = self.base_model.forward_features(x) #; print('layer conv1 ',x.size()) # [8, 64, 112, 112]\n        feature = self.avgpool(x).view(bs, -1)\n        y = self.last_linear(feature)\n\n        return y, feature\n    \n\nclass SEModule(nn.Module):\n\n    def __init__(self, channels, reduction):\n        super(SEModule, self).__init__()\n        self.avg_pool = nn.AdaptiveAvgPool1d(1)\n        self.fc1 = nn.Conv1d(channels, channels // reduction, kernel_size=1,\n                             padding=0)\n        self.relu = nn.ReLU(inplace=True)\n        self.fc2 = nn.Conv1d(channels // reduction, channels, kernel_size=1,\n                             padding=0)\n        self.sigmoid = nn.Sigmoid()\n\n    def forward(self, x):\n        module_input = x\n        x = self.avg_pool(x)\n        x = self.fc1(x)\n        x = self.relu(x)\n        x = self.fc2(x)\n        x = self.sigmoid(x)\n        return module_input * x\n    \nclass CNN_1D(nn.Module):\n\n    def __init__(self, num_classes=400, input_ch=1280, verbose=False):\n\n        super(CNN_1D, self).__init__()\n        pool = 4\n        drop = 0.1\n        self.verbose = verbose\n        self.layer1 = nn.Sequential(\n                nn.Conv1d(input_ch//pool, 64, kernel_size=7, stride=1, padding=3, bias=False),\n                nn.BatchNorm1d(64),\n                nn.ReLU(inplace=True),\n                SEModule(64, 16),\n#                 nn.Dropout(drop),\n        )\n        self.fpool = nn.MaxPool1d(kernel_size=pool, stride=pool, padding=0)\n        self.maxpool = nn.MaxPool1d(kernel_size=3, stride=2, padding=1)\n#         self.upsample = nn.Upsample(scale_factor=2, mode='bilinear')\n        self.layer2 = nn.Sequential(\n                nn.Conv1d(64, 128, kernel_size=3, stride=1, padding=1, bias=False),\n                nn.BatchNorm1d(128),\n                nn.ReLU(inplace=True),\n                SEModule(128, 16),\n#                 nn.Dropout(drop),\n        )\n        self.layer3 = nn.Sequential(\n                nn.Conv1d(128, 256, kernel_size=3, stride=1, padding=1, bias=False),\n                nn.BatchNorm1d(256),\n                nn.ReLU(inplace=True),\n                SEModule(256, 16),\n#                 nn.Dropout(drop),\n        )\n        self.layer4 = nn.Sequential(\n                nn.Conv1d(256, 512, kernel_size=3, stride=1, padding=1, bias=False),\n                nn.BatchNorm1d(512),\n                nn.ReLU(inplace=True),\n                SEModule(512, 16),\n#                 nn.Dropout(drop),\n        )\n        self.avgpool = nn.AdaptiveAvgPool1d(1)\n        self.fc2 = nn.Conv1d(\n            input_ch//pool+64+128+256+512, \n            2, kernel_size=1)\n#         self.fc = nn.Linear(512, 9)\n        self.fc = nn.Sequential(\n                nn.Linear(512, 512),\n                nn.ReLU(inplace=True),\n                nn.Dropout(0.5),\n                nn.Linear(512, 512),\n                nn.ReLU(inplace=True),\n                nn.Dropout(0.5),\n                nn.Linear(512, 9),\n        )\n\n    def forward(self, x_input):\n        bs, ch, d = x_input.size()\n        x0 = torch.transpose(x_input, 1, 2)\n        x0 = self.fpool(x0)\n        x0 = torch.transpose(x0, 1, 2)\n        x1 = self.layer1(x0)\n        x1 = self.maxpool(x1)\n\n        x2 = self.layer2(x1)\n        x2 = self.maxpool(x2)\n        x3 = self.layer3(x2)\n        x3 = self.maxpool(x3)\n        x4 = self.layer4(x3)\n        \n#         tmp = F.adaptive_avg_pool1d(x1, d)\n#         print(tmp.shape)\n#         tmp = F.adaptive_avg_pool1d(x2, d)\n#         print(tmp.shape)\n        x5 = torch.cat([\n            x0,\n            F.adaptive_avg_pool1d(x1, d), \n            F.adaptive_avg_pool1d(x2, d), \n            F.adaptive_avg_pool1d(x3, d), \n            F.adaptive_avg_pool1d(x4, d), \n        ], axis=1)\n        y2 = self.fc2(x5)\n        \n        b, ch, d = x_input.size()\n#         x1 = self.fc(x)\n#         x1 = x1.view(b, -1, 1)\n            \n        y = self.avgpool(x4)\n        y = y.view(b, -1)\n        y = self.fc(y)\n        return y, y2\n    \ncnn_2d = CNN_2D().to(device, non_blocking=True)\ncnn_2d.load_state_dict(torch.load(\"../input/rsna2020-pretrained-weights/b0_stage1/weight_epoch_16_fold1.pth\"))\n# cnn_2d.load_state_dict(torch.load(\"../input/rsna2020-1/201022_6_CNN_b0_1loss_512_fp16_bs80_lr1e3 (2)/201022_6_CNN_b0_1loss_512_fp16_bs80_lr1e3/weight_epoch_16_fold1.pth\"))\ncnn_1d = CNN_1D().to(device, non_blocking=True)\ncnn_1d.load_state_dict(torch.load(\"../input/rsna2020-pretrained-weights/b0_stage2/1dcnn_weight_best_fold1.pth\"))\n# cnn_1d.load_state_dict(torch.load(\"../input/rsna2020-1/201025_3_2ndNNs_features_pool_flip_new/201025_3_2ndNNs_features_pool_flip_new/cnn_weight_best_fold1.pth\"))\nlastfunc = nn.Sigmoid().to(device, non_blocking=True)\ncnn_2d.eval()\ncnn_1d.eval()\n\n# function definition\ndef batch_padding(batch):\n    bs, ch, d = batch.shape\n    d_new = int(np.ceil(d/64)*64)\n#     d_new = int(np.ceil(1083/64)*64)\n    batch_new = torch.from_numpy(np.zeros([bs, ch, d_new], np.float32)).to(device, non_blocking=True)\n    batch_new[:, :, :d] = batch\n    return batch_new\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 get_demo(idx):\n    # get exam data\n    exam = df_valid_exam[col_groupby][idx]\n    print(\"index: {}\".format(idx))\n    print(\"{}: {}\".format(col_groupby, exam))\n    start_index = df_valid_exam['start_index'][idx]\n    end_index = start_index + df_valid_exam['num_series'][idx]\n    df_tmp = df_valid.iloc[start_index:end_index].reset_index(drop=True)\n    images = get_images(idx)\n    \n    \n    # model prediction\n    BATCH_SIZE = 64\n    batchs = torch.from_numpy(images).to(device, non_blocking=True)\n    num_batches = int(np.ceil(images.shape[0]/BATCH_SIZE))\n    features = []\n    output0s = []\n    for batch_index in range(num_batches):\n        with torch.no_grad():\n            with torch.cuda.amp.autocast():\n                batch = batchs[batch_index*BATCH_SIZE:(batch_index+1)*BATCH_SIZE].to(device, non_blocking=True)\n                output0, feature = cnn_2d.feature(batch)\n                output0 = lastfunc(output0)\n        features.append(feature)\n        output0s.append(output0)\n    features = torch.cat(features, axis=0) # bs=d, ch\n    features = torch.transpose(features, 0,1).reshape([1, num_features1_1, -1])\n    features = batch_padding(features)\n    output0s = torch.cat(output0s, axis=0).data.cpu().numpy()\n    with torch.no_grad():\n        output1, output2 = cnn_1d(features)\n        output2 = output2[:,-1:]\n        output1 = lastfunc(output1).data.cpu().numpy()[0]\n\n        output2 = lastfunc(output2)[:,:,:len(images)].data.cpu().numpy()[0,0]\n    \n    # print exam-level pred\n    print('Exam-Level Labels         Prediction    True')\n    for i, item in enumerate(col_targets[:-1]):\n        col_name = \"{}                          \".format(col_targets_display[i])[:25]\n        print(\"{} {:.6f}      {}\".format(col_name, output1[i], df_valid_exam[col_targets[i]][idx]))\n    \n    # show exam-level pred\n    for i in range(9):\n        plt.figure(figsize=(20,5))\n        plt.subplot(1,3,1)\n        plt.scatter(np.arange(len(df_tmp)), output0s[:,i], alpha=1, s=3, label=\"Pred {}\".format(col_targets[i]))\n        plt.scatter(output0s[:,i].argmax(), output0s[:,i].max(), alpha=1, s=50, label='Max')\n        plt.scatter(output0s[:,i].argmin(), output0s[:,i].min(), alpha=1, s=50, label='Min')\n        plt.ylim(-0.1,1.1)\n        plt.legend(fontsize=12)\n        plt.grid()\n        plt.ylabel('score', fontsize=12)\n        plt.xlabel('series index', fontsize=12)\n        plt.title(\"Label:                  {}\\nExam-lebel Pred: {:.6f}\\nTrue:                    {}\".format(\n            col_targets_display[i], output1[i], df_valid_exam[col_targets[i]][idx]), fontsize=18, loc='left',\n                 color=('tab:red' if df_valid_exam[col_targets[i]][idx]==1 else 'tab:blue')\n                 )\n        plt.subplot(1,3,2)\n        image = images[output0s[:,i].argmax(),0].astype(np.float32)\n        image = window(image*1000, 100, 700)\n        plt.imshow(image, cmap='gray')\n        plt.title(\"Max\\nSeries index:        {}\\nImage-level Pred: {:.6f}\".format(output0s[:,i].argmax(), output0s[:,i].max()), fontsize=18, loc='left')\n        plt.subplot(1,3,3)\n        image = images[output0s[:,i].argmin(),0].astype(np.float32)\n        image = window(image*1000, 100, 700)\n        plt.imshow(image, cmap='gray')\n        plt.title(\"Min\\nSeries index:        {}\\nImage-level Pred: {:.6f}\".format(output0s[:,i].argmin(), output0s[:,i].min()), fontsize=18, loc='left')\n        plt.show()\n    \n    # show image-level pred\n    plt.figure(figsize=(20,5))\n    plt.subplot(1,3,1)\n    plt.scatter(np.arange(len(df_tmp)), output2, alpha=1, s=3, label=\"Pred PE present\", color='#1f77b4')\n    plt.scatter(output2.argmax(), output2.max(), alpha=1, s=50, label='Max', color='#ff7f0e')\n    plt.scatter(output2.argmin(), output2.min(), alpha=1, s=50, label='Min', color='#2ca02c')\n    plt.scatter(np.arange(len(df_tmp)), df_tmp[col_targets[-1]], alpha=1, s=10, label=\"True PE present\", color='tab:red', zorder=-1)\n    plt.ylim(-0.1,1.1)\n    plt.legend(fontsize=12)\n    plt.grid()\n    plt.ylabel('score', fontsize=12)\n    plt.xlabel('series index', fontsize=12)\n    plt.title(\"Label: {}\\n\".format(col_targets_display[-1]), fontsize=18, loc='left')\n    plt.subplot(1,3,2)\n    image = images[output2.argmax(),0].astype(np.float32)\n    image = window(image*1000, 100, 700)\n    plt.imshow(image, cmap='gray')\n    plt.title(\"Max\\nSeries index:        {}\\nImage-level Pred: {:.6f}\".format(output2.argmax(), output2.max()), fontsize=18, loc='left')\n    plt.subplot(1,3,3)\n    image = images[output2.argmin(),0].astype(np.float32)\n    image = window(image*1000, 100, 700)\n    plt.imshow(image, cmap='gray')\n    plt.title(\"Min\\nSeries index:        {}\\nImage-level Pred: {:.6f}\".format(output2.argmin(), output2.min()), fontsize=18, loc='left')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Run the cell below.  \nIt shows a visualization of a model prediction of 1 CT scan selected at random."},{"metadata":{"trusted":true},"cell_type":"code","source":"idx = np.random.randint(len(df_valid_exam))# sampling at random\n# idx = np.random.choice(df_valid_exam[df_valid_exam['3class']==1].index.values) # sampling from PE-positive scan\n# idx = np.choice(df_valid_exam[df_valid_exam['3class']==0]) # sampling from PE-negative scan\nidx = 878 # please comment out this line to sample CT scan at random\nget_demo(idx)","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}