{"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":"## References:\n* https://www.kaggle.com/xuxu1234/efficientnet3d-for-mri","metadata":{}},{"cell_type":"code","source":"import os\nimport glob\nfrom tqdm import tqdm_notebook as tqdm\nimport random\nimport numpy as np\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.nn.functional as F\nfrom torchvision import transforms, utils\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-10-12T03:49:11.077949Z","iopub.execute_input":"2021-10-12T03:49:11.078287Z","iopub.status.idle":"2021-10-12T03:49:15.785176Z","shell.execute_reply.started":"2021-10-12T03:49:11.078206Z","shell.execute_reply":"2021-10-12T03:49:15.784363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def set_seed(seed):\n    random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed_all(seed)\n        torch.backends.cudnn.deterministic = True\n\n\nset_seed(42)","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:15.786628Z","iopub.execute_input":"2021-10-12T03:49:15.786882Z","iopub.status.idle":"2021-10-12T03:49:15.835291Z","shell.execute_reply.started":"2021-10-12T03:49:15.786847Z","shell.execute_reply":"2021-10-12T03:49:15.834513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification'\ntrain_data = pd.read_csv(os.path.join(path, 'train_labels.csv'))\nprint('Num of train samples:', len(train_data))\ntrain_data.head()\nimg_size = 128","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:15.836799Z","iopub.execute_input":"2021-10-12T03:49:15.837387Z","iopub.status.idle":"2021-10-12T03:49:15.854772Z","shell.execute_reply.started":"2021-10-12T03:49:15.837347Z","shell.execute_reply":"2021-10-12T03:49:15.853549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def dicom2array(path, voi_lut=True, fix_monochrome=True):\n    dicom = pydicom.read_file(path)\n    # VOI LUT (if available by DICOM device) is used to\n    # transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    data = cv2.resize(data, (img_size, img_size))\n    return data\n\ndef load_3d_dicom_images(scan_id, split = \"train\"):\n    \"\"\"\n    we will use some heuristics to choose the slices to avoid any numpy zero matrix (if possible)\n    \"\"\"\n    flair = sorted(glob.glob(f\"{path}/{split}/{scan_id}/FLAIR/*.dcm\"))\n    t1w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1w/*.dcm\"))\n    t1wce = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1wCE/*.dcm\"))\n    t2w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T2w/*.dcm\"))\n    \n    \n    flair_img = np.array([dicom2array(a) for a in flair[len(flair)//2 - 25:len(flair)//2 + 25]]).T\n    \n    if flair_img.shape[-1] < 50:\n        n_zero = 50 - flair_img.shape[-1]\n        flair_img = np.concatenate((flair_img, np.zeros((img_size, img_size, n_zero))), axis = -1)\n    #print(flair_img.shape)\n        \n    \n    \n    t1w_img = np.array([dicom2array(a) for a in t1w[len(t1w)//2 - 25:len(t1w)//2 + 25]]).T\n    if t1w_img.shape[-1] < 50:\n        n_zero = 50 - t1w_img.shape[-1]\n        t1w_img = np.concatenate((t1w_img, np.zeros((img_size, img_size, n_zero))), axis = -1)\n    #print(t1w_img.shape)\n    \n    \n    t1wce_img = np.array([dicom2array(a) for a in t1wce[len(t1wce)//2 - 25:len(t1wce)//2 + 25]]).T\n    if t1wce_img.shape[-1] < 50:\n        n_zero = 50 - t1wce_img.shape[-1]\n        t1wce_img = np.concatenate((t1wce_img, np.zeros((img_size, img_size, n_zero))), axis = -1)\n    #print(t1wce_img.shape)\n    \n    \n    t2w_img = np.array([dicom2array(a) for a in t2w[len(t2w)//2 - 25:len(t2w)//2 + 25]]).T\n    if t2w_img.shape[-1] < 50:\n        n_zero = 50 - t2w_img.shape[-1]\n        t2w_img = np.concatenate((t2w_img, np.zeros((img_size, img_size, n_zero))), axis = -1)\n    #print(t2w_img.shape)\n    \n    return np.concatenate((flair_img, t1w_img, t1wce_img, t2w_img), axis = -1)","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:15.856786Z","iopub.execute_input":"2021-10-12T03:49:15.857219Z","iopub.status.idle":"2021-10-12T03:49:15.8745Z","shell.execute_reply.started":"2021-10-12T03:49:15.857182Z","shell.execute_reply":"2021-10-12T03:49:15.873781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices = load_3d_dicom_images(\"00000\")\nprint(slices.shape)","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:15.875569Z","iopub.execute_input":"2021-10-12T03:49:15.875833Z","iopub.status.idle":"2021-10-12T03:49:17.881284Z","shell.execute_reply.started":"2021-10-12T03:49:15.875797Z","shell.execute_reply":"2021-10-12T03:49:17.880498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class BrainTumor(Dataset):\n    def __init__(self, path = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification', split = \"train\", validation_split = 0.0):\n        # labels\n        train_data = pd.read_csv(os.path.join(path, 'train_labels.csv'))\n        self.labels = {}\n        brats = list(train_data[\"BraTS21ID\"])\n        mgmt = list(train_data[\"MGMT_value\"])\n        for b, m in zip(brats, mgmt):\n            self.labels[str(b).zfill(5)] = m\n            \n        if split == \"valid\":\n            self.split = split\n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{split}/\" + \"/*\"))]\n            self.ids = self.ids[:int(len(self.ids)* validation_split)] # first 20% as validation\n        elif split == \"train\":\n            self.split = split\n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{split}/\" + \"/*\"))]\n            self.ids = self.ids[int(len(self.ids)* validation_split):] # last 80% as train\n        else:\n            self.split = split\n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{split}/\" + \"/*\"))]\n            \n    \n    def __len__(self):\n        return len(self.ids)\n    \n    def __getitem__(self, idx):\n        imgs = load_3d_dicom_images(self.ids[idx], self.split)\n        transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize((0.5,) * 200, (0.5,) * 200)])\n        imgs = transform(imgs)\n        \n        if self.split != \"test\":\n            label = self.labels[self.ids[idx]]\n            return torch.tensor(imgs, dtype = torch.float32), torch.tensor(label, dtype = torch.long)\n        else:\n            return torch.tensor(imgs, dtype = torch.float32)","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:17.882405Z","iopub.execute_input":"2021-10-12T03:49:17.883106Z","iopub.status.idle":"2021-10-12T03:49:17.897782Z","shell.execute_reply.started":"2021-10-12T03:49:17.883052Z","shell.execute_reply":"2021-10-12T03:49:17.896946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dataset = BrainTumor()\ntrain_loader = DataLoader(train_dataset, batch_size=2, shuffle=True, num_workers=4)","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:17.899028Z","iopub.execute_input":"2021-10-12T03:49:17.899292Z","iopub.status.idle":"2021-10-12T03:49:17.995339Z","shell.execute_reply.started":"2021-10-12T03:49:17.899261Z","shell.execute_reply":"2021-10-12T03:49:17.994607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for img, label in train_loader:\n    print(img.shape)\n    print(label.shape)\n    break","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:17.997201Z","iopub.execute_input":"2021-10-12T03:49:17.997386Z","iopub.status.idle":"2021-10-12T03:49:25.855376Z","shell.execute_reply.started":"2021-10-12T03:49:17.997364Z","shell.execute_reply":"2021-10-12T03:49:25.854473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MODEL","metadata":{}},{"cell_type":"code","source":"#https://github.com/MontaEllis/Pytorch-Medical-Classification/blob/main/models/three_d/densenet3d.py\n\nimport math\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom collections import OrderedDict\n\n\nclass _DenseLayer(nn.Sequential):\n\n    def __init__(self, num_input_features, growth_rate, bn_size, drop_rate):\n        super().__init__()\n        self.add_module('norm1', nn.BatchNorm3d(num_input_features))\n        self.add_module('relu1', nn.ReLU(inplace=True))\n        self.add_module(\n            'conv1',\n            nn.Conv3d(num_input_features,\n                      bn_size * growth_rate,\n                      kernel_size=1,\n                      stride=1,\n                      bias=False))\n        self.add_module('norm2', nn.BatchNorm3d(bn_size * growth_rate))\n        self.add_module('relu2', nn.ReLU(inplace=True))\n        self.add_module(\n            'conv2',\n            nn.Conv3d(bn_size * growth_rate,\n                      growth_rate,\n                      kernel_size=3,\n                      stride=1,\n                      padding=1,\n                      bias=False))\n        self.drop_rate = drop_rate\n\n    def forward(self, x):\n        new_features = super().forward(x)\n        if self.drop_rate > 0:\n            new_features = F.dropout(new_features,\n                                     p=self.drop_rate,\n                                     training=self.training)\n        return torch.cat([x, new_features], 1)\n\n\nclass _DenseBlock(nn.Sequential):\n\n    def __init__(self, num_layers, num_input_features, bn_size, growth_rate,\n                 drop_rate):\n        super().__init__()\n        for i in range(num_layers):\n            layer = _DenseLayer(num_input_features + i * growth_rate,\n                                growth_rate, bn_size, drop_rate)\n            self.add_module('denselayer{}'.format(i + 1), layer)\n\n\nclass _Transition(nn.Sequential):\n\n    def __init__(self, num_input_features, num_output_features):\n        super().__init__()\n        self.add_module('norm', nn.BatchNorm3d(num_input_features))\n        self.add_module('relu', nn.ReLU(inplace=True))\n        self.add_module(\n            'conv',\n            nn.Conv3d(num_input_features,\n                      num_output_features,\n                      kernel_size=1,\n                      stride=1,\n                      bias=False))\n        self.add_module('pool', nn.AvgPool3d(kernel_size=2, stride=2))\n\n\nclass DenseNet(nn.Module):\n    \"\"\"Densenet-BC model class\n    Args:\n        growth_rate (int) - how many filters to add each layer (k in paper)\n        block_config (list of 4 ints) - how many layers in each pooling block\n        num_init_features (int) - the number of filters to learn in the first convolution layer\n        bn_size (int) - multiplicative factor for number of bottle neck layers\n          (i.e. bn_size * k features in the bottleneck layer)\n        drop_rate (float) - dropout rate after each dense layer\n        num_classes (int) - number of classification classes\n    \"\"\"\n\n    def __init__(self,\n                 n_input_channels=1,\n                 conv1_t_size=7,\n                 conv1_t_stride=1,\n                 no_max_pool=False,\n                 growth_rate=32,\n                 block_config=(6, 12, 24, 16),\n                 num_init_features=64,\n                 bn_size=4,\n                 drop_rate=0,\n                 num_classes=2):\n\n        super().__init__()\n\n        # First convolution\n        self.features = [('conv1',\n                          nn.Conv3d(n_input_channels,\n                                    num_init_features,\n                                    kernel_size=(conv1_t_size, 7, 7),\n                                    stride=(conv1_t_stride, 2, 2),\n                                    padding=(conv1_t_size // 2, 3, 3),\n                                    bias=False)),\n                         ('norm1', nn.BatchNorm3d(num_init_features)),\n                         ('relu1', nn.ReLU(inplace=True))]\n        if not no_max_pool:\n            self.features.append(\n                ('pool1', nn.MaxPool3d(kernel_size=3, stride=2, padding=1)))\n        self.features = nn.Sequential(OrderedDict(self.features))\n\n        # Each denseblock\n        num_features = num_init_features\n        for i, num_layers in enumerate(block_config):\n            block = _DenseBlock(num_layers=num_layers,\n                                num_input_features=num_features,\n                                bn_size=bn_size,\n                                growth_rate=growth_rate,\n                                drop_rate=drop_rate)\n            self.features.add_module('denseblock{}'.format(i + 1), block)\n            num_features = num_features + num_layers * growth_rate\n            if i != len(block_config) - 1:\n                trans = _Transition(num_input_features=num_features,\n                                    num_output_features=num_features // 2)\n                self.features.add_module('transition{}'.format(i + 1), trans)\n                num_features = num_features // 2\n\n        # Final batch norm\n        self.features.add_module('norm5', nn.BatchNorm3d(num_features))\n\n        for m in self.modules():\n            if isinstance(m, nn.Conv3d):\n                m.weight = nn.init.kaiming_normal(m.weight, mode='fan_out')\n            elif isinstance(m, nn.BatchNorm3d) or isinstance(m, nn.BatchNorm2d):\n                m.weight.data.fill_(1)\n                m.bias.data.zero_()\n\n        # Linear layer\n        self.classifier = nn.Linear(num_features, num_classes)\n\n        for m in self.modules():\n            if isinstance(m, nn.Conv3d):\n                nn.init.kaiming_normal_(m.weight,\n                                        mode='fan_out',\n                                        nonlinearity='relu')\n            elif isinstance(m, nn.BatchNorm3d):\n                nn.init.constant_(m.weight, 1)\n                nn.init.constant_(m.bias, 0)\n            elif isinstance(m, nn.Linear):\n                nn.init.constant_(m.bias, 0)\n\n    def forward(self, x):\n        features = self.features(x)\n        out = F.relu(features, inplace=True)\n        out = F.adaptive_avg_pool3d(out,\n                                    output_size=(1, 1,\n                                                 1)).view(features.size(0), -1)\n        out = self.classifier(out)\n        return out\n\n\ndef generate_model(model_depth, **kwargs):\n    assert model_depth in [121, 169, 201, 264]\n\n    if model_depth == 121:\n        model = DenseNet(num_init_features=64,\n                         growth_rate=32,\n                         block_config=(6, 12, 24, 16),\n                         **kwargs)\n    elif model_depth == 169:\n        model = DenseNet(num_init_features=64,\n                         growth_rate=32,\n                         block_config=(6, 12, 32, 32),\n                         **kwargs)\n    elif model_depth == 201:\n        model = DenseNet(num_init_features=64,\n                         growth_rate=32,\n                         block_config=(6, 12, 48, 32),\n                         **kwargs)\n    elif model_depth == 264:\n        model = DenseNet(num_init_features=64,\n                         growth_rate=32,\n                         block_config=(6, 12, 64, 48),\n                         **kwargs)\n\n    return model\n\n\nif __name__ == \"__main__\":\n\n\n    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n    image_size = 128\n    x = torch.Tensor(1, 1, image_size, image_size, image_size)\n    x = x.to(device)\n    print(\"x size: {}\".format(x.size()))\n    \n    model = generate_model(201,n_input_channels=1,num_classes=2).to(device)\n    \n\n    out1 = model(x)\n    print(\"out size: {}\".format(out1.size()))","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:25.857477Z","iopub.execute_input":"2021-10-12T03:49:25.857759Z","iopub.status.idle":"2021-10-12T03:49:32.469936Z","shell.execute_reply.started":"2021-10-12T03:49:25.857723Z","shell.execute_reply":"2021-10-12T03:49:32.469135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class FocalLoss(nn.Module):\n    def __init__(self, alpha=1., gamma=1.):\n        super().__init__()\n        self.alpha = alpha\n        self.gamma = gamma\n\n    def forward(self, inputs, targets, **kwargs):\n        CE_loss = nn.CrossEntropyLoss(reduction='none')(inputs, targets)\n        pt = torch.exp(-CE_loss)\n        F_loss = self.alpha * ((1-pt)**self.gamma) * CE_loss\n        return F_loss.mean()\n\n    \ndef roc_score(inp, target):\n    _, indices = inp.max(1)\n    return torch.Tensor([roc_auc_score(target, indices)])[0]","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:32.472493Z","iopub.execute_input":"2021-10-12T03:49:32.473213Z","iopub.status.idle":"2021-10-12T03:49:32.481508Z","shell.execute_reply.started":"2021-10-12T03:49:32.473176Z","shell.execute_reply":"2021-10-12T03:49:32.480577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = DenseNet()\ncriterion = FocalLoss()\noptimizer = torch.optim.Adam(model.parameters(),lr = 0.0001)\nn_epochs = 10","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:32.482966Z","iopub.execute_input":"2021-10-12T03:49:32.483822Z","iopub.status.idle":"2021-10-12T03:49:32.978094Z","shell.execute_reply.started":"2021-10-12T03:49:32.483783Z","shell.execute_reply":"2021-10-12T03:49:32.977177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## TRAIN","metadata":{}},{"cell_type":"code","source":"gpu = torch.device(f\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nmodel.to(gpu)\n\nfor epoch in range(n_epochs):  # loop over the dataset multiple times\n\n    train_loss = []\n    best_pres = 10000\n    model.train()\n    for i, data in tqdm(enumerate(train_loader, 0)):\n        x, y = data\n        \n        x = torch.unsqueeze(x, dim = 1)\n        x = x.to(gpu)\n        y = y.to(gpu)\n\n        # zero the parameter gradients\n        optimizer.zero_grad()\n\n        # forward + backward + optimize\n        outputs = model(x)\n        loss = criterion(outputs, y)\n        loss.backward()\n        optimizer.step()\n\n        # print statistics\n        train_loss.append(loss.item())\n    avg_train = sum(train_loss) / len(train_loss)\n    print(f\"epoch {epoch+1} train: {avg_train}\")\n\n    if avg_train < best_pres:\n        print('save model...')\n        best_pres = avg_train\n        torch.save(model.state_dict(),'best_Densenet_201_loss.pt')","metadata":{"execution":{"iopub.status.busy":"2021-10-12T03:49:32.982427Z","iopub.execute_input":"2021-10-12T03:49:32.982701Z","iopub.status.idle":"2021-10-12T03:56:24.977861Z","shell.execute_reply.started":"2021-10-12T03:49:32.982668Z","shell.execute_reply":"2021-10-12T03:56:24.975483Z"},"trusted":true},"execution_count":null,"outputs":[]}]}