{"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":"## **Problem Description**:\n\nThere are structural multi-parametric MRI (mpMRI) scans for different subjects, in DICOM format. The exact mpMRI scans included are:\n\n* Fluid Attenuated Inversion Recovery (FLAIR)\n* T1-weighted pre-contrast (T1w)\n* T1-weighted post-contrast (T1Gd)\n* T2-weighted (T2)\n\n`train_labels.csv` - file contains the target MGMT_value for each subject in the training data **(e.g. the presence of MGMT promoter methylation)**.\n\nSo, it's a binary classification problem.\n\n## **MRI Augmentation**\n\n* https://arxiv.org/ftp/arxiv/papers/2006/2006.01693.pdf\n* https://arxiv.org/pdf/2106.14947.pdf\n* https://arxiv.org/pdf/1910.08112.pdf\n* https://arxiv.org/pdf/1902.09383.pdf\n* https://link.springer.com/chapter/10.1007/978-3-030-00536-8_1\n* https://arxiv.org/pdf/1805.09501v1.pdf (AutoAugment)\n","metadata":{}},{"cell_type":"code","source":"import sys\nsys.path.append('../input/autoaugment')\nsys.path.append('../input/imgaug/imgaug-master')","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:24:24.239364Z","iopub.execute_input":"2021-07-24T22:24:24.239861Z","iopub.status.idle":"2021-07-24T22:24:24.243959Z","shell.execute_reply.started":"2021-07-24T22:24:24.239831Z","shell.execute_reply":"2021-07-24T22:24:24.243320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **Importing libraries**","metadata":{}},{"cell_type":"code","source":"import os\nimport glob\nfrom tqdm import tqdm_notebook as tqdm\nimport math\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\nfrom sklearn.metrics import roc_auc_score\nimport imgaug as ia\nimport imgaug.augmenters as iaa\nfrom autoaugment import ImageNetPolicy\nfrom PIL import Image\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:55:46.740428Z","iopub.execute_input":"2021-07-24T22:55:46.740817Z","iopub.status.idle":"2021-07-24T22:55:46.748763Z","shell.execute_reply.started":"2021-07-24T22:55:46.740785Z","shell.execute_reply":"2021-07-24T22:55:46.747753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **Dicom Image Preprocessing**","metadata":{}},{"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))","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:32:25.544699Z","iopub.execute_input":"2021-07-24T22:32:25.545059Z","iopub.status.idle":"2021-07-24T22:32:25.563177Z","shell.execute_reply.started":"2021-07-24T22:32:25.545030Z","shell.execute_reply":"2021-07-24T22:32:25.561892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def dicom2array(path, voi_lut=True, fix_monochrome=True, remove_black_boundary=True, augmentation=False):\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    if remove_black_boundary: # we get slightly more details\n        (x, y) = np.where(data > 0)\n        if len(x) > 0 and len(y) > 0:\n            x_mn = np.min(x)\n            x_mx = np.max(x)\n            y_mn = np.min(y)\n            y_mx = np.max(y)\n            if (x_mx - x_mn) > 10 and (y_mx - y_mn) > 10:\n                data = data[:,np.min(y):np.max(y)]\n    data = cv2.resize(data, (512, 512))\n    return data\n\ndef load_rand_dicom_images(scan_id, split = \"train\"):\n    \"\"\"\n    send 4 random slices of each modality\n    \"\"\"\n    flair = sorted(glob.glob(f\"{path}/{split}/{scan_id}/FLAIR/*.dcm\"))\n    flair_img = dicom2array(random.sample(flair, 1)[0])\n    t1w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1w/*.dcm\"))\n    t1w_img = dicom2array(random.sample(t1w, 1)[0])\n    t1wce = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1wCE/*.dcm\"))\n    t1wce_img = dicom2array(random.sample(t1wce, 1)[0])\n    t2w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T2w/*.dcm\"))\n    t2w_img = dicom2array(random.sample(t2w, 1)[0])\n    \n    return np.array((flair_img, t1w_img, t1wce_img, t2w_img)).T","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:32:25.924499Z","iopub.execute_input":"2021-07-24T22:32:25.924837Z","iopub.status.idle":"2021-07-24T22:32:25.938437Z","shell.execute_reply.started":"2021-07-24T22:32:25.924809Z","shell.execute_reply":"2021-07-24T22:32:25.937349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"load_rand_dicom_images(\"00000\").shape","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:32:27.076777Z","iopub.execute_input":"2021-07-24T22:32:27.077159Z","iopub.status.idle":"2021-07-24T22:32:27.256404Z","shell.execute_reply.started":"2021-07-24T22:32:27.077122Z","shell.execute_reply":"2021-07-24T22:32:27.255456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_imgs(imgs, cols=4, size=7, is_rgb=True, title=\"\", cmap='gray', img_size=(512,512)):\n    rows = len(imgs)//cols + 1\n    fig = plt.figure(figsize=(cols*size, rows*size))\n    for i in range(4):\n        img = imgs[:,:,i]\n        fig.add_subplot(rows, cols, i+1)\n        plt.imshow(img, cmap=cmap)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:38:25.559411Z","iopub.execute_input":"2021-07-24T22:38:25.559762Z","iopub.status.idle":"2021-07-24T22:38:25.566495Z","shell.execute_reply.started":"2021-07-24T22:38:25.559731Z","shell.execute_reply":"2021-07-24T22:38:25.565465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices = load_rand_dicom_images(\"00000\")\nplot_imgs(slices, title=\"Original MRI Slices\")","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:45:09.895192Z","iopub.execute_input":"2021-07-24T22:45:09.895579Z","iopub.status.idle":"2021-07-24T22:45:10.920742Z","shell.execute_reply.started":"2021-07-24T22:45:09.895544Z","shell.execute_reply":"2021-07-24T22:45:10.919489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **Augmentation**","metadata":{}},{"cell_type":"code","source":"slices.T.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:48:07.291772Z","iopub.execute_input":"2021-07-24T22:48:07.292261Z","iopub.status.idle":"2021-07-24T22:48:07.297697Z","shell.execute_reply.started":"2021-07-24T22:48:07.292228Z","shell.execute_reply":"2021-07-24T22:48:07.296891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nsometimes = lambda aug: iaa.Sometimes(0.2, aug)\n\nseq = iaa.Sequential(\n    [\n        # apply the following augmenters to most images\n        iaa.Fliplr(0.5), # horizontally flip 50% of all images\n        iaa.Flipud(0.2), # vertically flip 20% of all images\n        # crop images by -5% to 10% of their height/width\n        sometimes(iaa.CropAndPad(\n            percent=(-0.05, 0.05),\n            pad_mode=ia.ALL,\n            pad_cval=(0, 255)\n        )),\n        sometimes(iaa.Affine(\n            scale={\"x\": (0.8, 1.2), \"y\": (0.8, 1.2)}, # scale images to 80-120% of their size, individually per axis\n            translate_percent={\"x\": (-0.2, 0.2), \"y\": (-0.2, 0.2)}, # translate by -20 to +20 percent (per axis)\n            rotate=(-45, 45), # rotate by -45 to +45 degrees\n            shear=(-16, 16), # shear by -16 to +16 degrees\n            order=[0, 1], # use nearest neighbour or bilinear interpolation (fast)\n            cval=(0, 255), # if mode is constant, use a cval between 0 and 255\n            mode=ia.ALL # use any of scikit-image's warping modes (see 2nd image from the top for examples)\n        )),\n        # execute 0 to 5 of the following (less important) augmenters per image\n        # don't execute all of them, as that would often be way too strong\n        iaa.SomeOf((0, 5),\n            [\n                sometimes(iaa.Superpixels(p_replace=(0, 1.0), n_segments=(20, 200))), # convert images into their superpixel representation\n                iaa.OneOf([\n                    iaa.GaussianBlur((0, 3.0)), # blur images with a sigma between 0 and 3.0\n                    iaa.AverageBlur(k=(2, 7)), # blur image using local means with kernel sizes between 2 and 7\n                    iaa.MedianBlur(k=(3, 11)), # blur image using local medians with kernel sizes between 2 and 7\n                ]),\n                iaa.Sharpen(alpha=(0, 1.0), lightness=(0.75, 1.5)), # sharpen images\n                iaa.Emboss(alpha=(0, 1.0), strength=(0, 2.0)), # emboss images\n                # search either for all edges or for directed edges,\n                # blend the result with the original image using a blobby mask\n                iaa.SimplexNoiseAlpha(iaa.OneOf([\n                    iaa.EdgeDetect(alpha=(0.5, 1.0)),\n                    iaa.DirectedEdgeDetect(alpha=(0.5, 1.0), direction=(0.0, 1.0)),\n                ])),\n                iaa.AdditiveGaussianNoise(loc=0, scale=(0.0, 0.05*255), per_channel=0.5), # add gaussian noise to images\n                iaa.OneOf([\n                    iaa.Dropout((0.01, 0.1), per_channel=0.5), # randomly remove up to 10% of the pixels\n                    iaa.CoarseDropout((0.03, 0.15), size_percent=(0.02, 0.05), per_channel=0.2),\n                ]),\n                iaa.Invert(0.05, per_channel=True), # invert color channels\n                iaa.Add((-10, 10), per_channel=0.5), # change brightness of images (by -10 to 10 of original value)\n                \n                # either change the brightness of the whole image (sometimes\n                # per channel) or change the brightness of subareas\n                iaa.OneOf([\n                    iaa.Multiply((0.5, 1.5), per_channel=0.5),\n                    iaa.FrequencyNoiseAlpha(\n                        exponent=(-4, 0),\n                        first=iaa.Multiply((0.5, 1.5), per_channel=True),\n                        second=iaa.LinearContrast((0.5, 2.0))\n                    )\n                ]),\n                iaa.LinearContrast((0.5, 2.0), per_channel=0.5), # improve or worsen the contrast\n                sometimes(iaa.ElasticTransformation(alpha=(0.5, 3.5), sigma=0.25)), # move pixels locally around (with random strengths)\n                sometimes(iaa.PiecewiseAffine(scale=(0.01, 0.05))), # sometimes move parts of the image around\n                sometimes(iaa.PerspectiveTransform(scale=(0.01, 0.1)))\n            ],\n            random_order=True\n        )\n    ],\n    random_order=True\n)\nslices_aug = seq(images=slices.T)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:51:30.049928Z","iopub.execute_input":"2021-07-24T22:51:30.050478Z","iopub.status.idle":"2021-07-24T22:51:30.138783Z","shell.execute_reply.started":"2021-07-24T22:51:30.050429Z","shell.execute_reply":"2021-07-24T22:51:30.137977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices_aug.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:51:34.974009Z","iopub.execute_input":"2021-07-24T22:51:34.974580Z","iopub.status.idle":"2021-07-24T22:51:34.980026Z","shell.execute_reply.started":"2021-07-24T22:51:34.974532Z","shell.execute_reply":"2021-07-24T22:51:34.979210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_imgs(slices_aug.T, title=\"Augmented MRI Slices\")","metadata":{"execution":{"iopub.status.busy":"2021-07-24T22:52:01.734656Z","iopub.execute_input":"2021-07-24T22:52:01.735018Z","iopub.status.idle":"2021-07-24T22:52:02.750676Z","shell.execute_reply.started":"2021-07-24T22:52:01.734987Z","shell.execute_reply":"2021-07-24T22:52:02.746821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **AutoAugment**","metadata":{}},{"cell_type":"code","source":"def autoaugment(np_arr):\n    out_np_arr = []\n    for s in np_arr: # D X H X W\n        slice_pil = Image.fromarray(np.uint8(s)).convert('RGB')\n        policy = ImageNetPolicy()\n        transformed = policy(slice_pil)\n        out_np_arr.append(np.array(transformed.convert('L')))\n    return np.array(out_np_arr)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T23:06:21.657895Z","iopub.execute_input":"2021-07-24T23:06:21.658262Z","iopub.status.idle":"2021-07-24T23:06:21.663654Z","shell.execute_reply.started":"2021-07-24T23:06:21.658213Z","shell.execute_reply":"2021-07-24T23:06:21.662651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices_autoaug = autoaugment(slices.T)\nprint(slices_autoaug.shape)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T23:06:22.954713Z","iopub.execute_input":"2021-07-24T23:06:22.955057Z","iopub.status.idle":"2021-07-24T23:06:23.026646Z","shell.execute_reply.started":"2021-07-24T23:06:22.955027Z","shell.execute_reply":"2021-07-24T23:06:23.025919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_imgs(slices_autoaug.T, title=\"AutoAugmented MRI Slices\")","metadata":{"execution":{"iopub.status.busy":"2021-07-24T23:06:31.774907Z","iopub.execute_input":"2021-07-24T23:06:31.775426Z","iopub.status.idle":"2021-07-24T23:06:32.726084Z","shell.execute_reply.started":"2021-07-24T23:06:31.775393Z","shell.execute_reply":"2021-07-24T23:06:32.725240Z"},"trusted":true},"execution_count":null,"outputs":[]}]}