{"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-MICCAI Brain Tumor Radiogenomic Classification","metadata":{}},{"cell_type":"code","source":"import os\nimport json\nimport glob\nimport random\nimport collections\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport random\nfrom skimage.transform import resize\n\n\ntrain_df = pd.read_csv(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv\")\ntest_df = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv')\n\n\nTYPES = [\"FLAIR\",]\nWHITE_THRESHOLD = 10 # out of 255\nEXCLUDE = [109, 123, 709]","metadata":{"execution":{"iopub.status.busy":"2021-08-17T20:54:04.722339Z","iopub.execute_input":"2021-08-17T20:54:04.722848Z","iopub.status.idle":"2021-08-17T20:54:06.707757Z","shell.execute_reply.started":"2021-08-17T20:54:04.722743Z","shell.execute_reply":"2021-08-17T20:54:06.706558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom(path, size = 64):\n    ''' \n    Reads a DICOM image, standardizes so that the pixel values are between 0 and 1, then rescales to 0 and 255\n    \n    Note super sure if this kind of scaling is appropriate, but everyone seems to do it. \n    '''\n    dicom = pydicom.read_file(path)\n    data = dicom.pixel_array\n    if np.max(data) != 0:\n        data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return cv2.resize(data, (size, size))\n\n\ndef get_all_image_paths(brats21id, image_type, folder='train'): \n    '''\n    Returns an arry of all the images of a particular type for a particular patient ID\n    '''\n    assert(image_type in TYPES)\n    \n    patient_path = os.path.join(\n        \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/%s/\" % folder, \n        str(brats21id).zfill(5),\n    )\n\n    return sorted(\n        glob.glob(os.path.join(patient_path, image_type, \"*\")), \n        key=lambda x: int(x[:-4].split(\"-\")[-1]),\n    )\n    \ndef get_white_area(image, threshold): \n    '''\n    Given an image, this function computes the fraction of pixels for which the greyscale value is at least threshold. \n    '''\n    return np.sum(image > threshold) / image.shape[0] / image.shape[1]\n\ndef filter_images(images, quantiles):\n    '''\n    Filters a list of images based on the quantiles given. \n    \n    For example, if quantiles = [0.33, 0.66], then we would find the photos with 0.33 and 0.66 of the maximum white area.\n    \n    Note that the images get bigger (up to the maximum) and then get smaller. Thus, we would extract 5 photos\n    corresponding to [0.33, 0.66, 1, 0.66, 0.33]\n    '''\n    quantiles.sort()\n    white_areas = [get_white_area(im, WHITE_THRESHOLD) for im in images]\n    max_white_area = max(white_areas)\n    \n    middle_index = white_areas.index(max_white_area)\n    \n    indices_to_keep = []\n    \n    q = 0\n    for i in range(middle_index): \n        if q >= len(quantiles): \n            break\n        if white_areas[i] >= quantiles[q] * max_white_area: \n            indices_to_keep.append(i)\n            q += 1\n    \n    indices_to_keep.append(middle_index)\n    \n    q = len(quantiles) - 1\n    for i in range(middle_index, len(white_areas)): \n        if q < 0: \n            break\n        if white_areas[i] <= quantiles[q] * max_white_area: \n            indices_to_keep.append(i)\n            q -= 1\n    \n    # expected # of images vs. actual number of images\n    difference = len(quantiles) * 2 + 1 - len(indices_to_keep)\n    \n    if difference > 0: \n        indices_to_keep += [i] * difference\n    \n    return indices_to_keep\n\n# Plotting Functions\n\ndef plot_image_white(brats21id, image_type):\n    images_paths = get_all_image_paths(brats21id, image_type)\n    images = [load_dicom(im) for im in images_paths]\n    plt.xlabel('Image Number')\n    plt.ylabel('White Area')\n    areas = [get_white_area(x, WHITE_THRESHOLD) for x in images]\n    plt.plot(range(len(images)), areas)\n    plt.show()\n    return max(areas)\n\ndef plot_images_at_quantiles(images, quantile_labels):  \n    assert(len(images) == len(quantile_labels))\n    \n    plt.figure(figsize=((30 // len(quantile_labels)) * len(images), 10))\n\n    for i in range(len(images)):\n        plt.subplot(1, len(images), i + 1)\n        plt.imshow(images[i], cmap=\"gray\")\n        plt.title(f\"{quantile_labels[i]}\", fontsize=16)\n        plt.axis(\"off\")\n    plt.show()\n    return \n\ndef center_images_for_patient(images, quantiles, size):\n    #select the image with the most white space\n    image_main = images_filtered[quantile_labels.index(1)]\n    \n    #find binding box of the image\n\n    return resized_img","metadata":{"execution":{"iopub.status.busy":"2021-08-16T20:48:17.708414Z","iopub.execute_input":"2021-08-16T20:48:17.708803Z","iopub.status.idle":"2021-08-16T20:48:17.731546Z","shell.execute_reply.started":"2021-08-16T20:48:17.708768Z","shell.execute_reply":"2021-08-16T20:48:17.730497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# image_paths = np.array(get_all_image_paths(100, 'FLAIR', 'train'))\n# images = np.array([load_dicom(im) for im in image_paths])\n\n# quantiles = [0.33, 0.67]\n# quantiles.sort()\n# quantile_labels = quantiles + [1] + quantiles[::-1]\n\n# images_filtered = images[filter_images(images, quantiles)]\n\n# plot_images_at_quantiles(images_filtered, quantile_labels)","metadata":{"execution":{"iopub.status.busy":"2021-08-13T19:25:22.685733Z","iopub.execute_input":"2021-08-13T19:25:22.686128Z","iopub.status.idle":"2021-08-13T19:25:23.067522Z","shell.execute_reply.started":"2021-08-13T19:25:22.686096Z","shell.execute_reply":"2021-08-13T19:25:23.066706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def filter_images_for_patient(brats21id, quantiles, folder='train'):\n    images_filtered = {}\n    for t in TYPES:\n        image_paths = np.array(get_all_image_paths(brats21id, t, folder))\n        \n        if len(image_paths) > 20: \n            image_paths = image_paths[:: len(image_paths) // 20]\n        images = np.array([load_dicom(im) for im in image_paths])\n        images = images[filter_images(images, quantiles)]\n        \n        image_main = images[len(quantiles)]\n        \n        col_sum = np.where(np.sum(image_main, axis = 0)>0)\n        row_sum = np.where(np.sum(image_main, axis = 1)>0)\n        y1, y2 = row_sum[0][0], row_sum[0][-1]\n        x1, x2 = col_sum[0][0], col_sum[0][-1]\n\n        images_filtered[t] = [resize(im[y1:y2, x1:x2], (64, 64), anti_aliasing=True) for im in images]        \n\n    return images_filtered","metadata":{"execution":{"iopub.status.busy":"2021-08-16T20:49:59.343111Z","iopub.execute_input":"2021-08-16T20:49:59.343594Z","iopub.status.idle":"2021-08-16T20:49:59.352621Z","shell.execute_reply.started":"2021-08-16T20:49:59.343564Z","shell.execute_reply":"2021-08-16T20:49:59.351789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_data_for_patients(patient_ids, folder='train'):\n    output = np.array([]) \n    count = 0 \n    for patient_id in patient_ids:\n        if patient_id in EXCLUDE: \n            continue\n        images = filter_images_for_patient(patient_id, [0.33, 0.66], folder)\n        output = np.append(output, images)\n        count += 1\n        if count % 100 == 99: \n            print('Done with %d out of %d' % (count + 1, len(patient_ids)))\n    \n    return output\n\ndef flatten_data_for_individual(data):\n    all_types = [data[t] for t in TYPES]\n    return np.array([item for sublist in all_types for item in sublist])\n\ndef format_data_for_keras(all_data):\n    result =  np.array([flatten_data_for_individual(data) for data in all_data])\n    result = np.swapaxes(result, 1, 2)\n    result = np.swapaxes(result, 2, 3)\n    return result\n\n","metadata":{"execution":{"iopub.status.busy":"2021-08-16T20:50:00.094462Z","iopub.execute_input":"2021-08-16T20:50:00.094847Z","iopub.status.idle":"2021-08-16T20:50:00.103619Z","shell.execute_reply.started":"2021-08-16T20:50:00.094813Z","shell.execute_reply":"2021-08-16T20:50:00.102606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_data_raw = get_data_for_patients(train_df.BraTS21ID)\ntesting_data_raw = get_data_for_patients(test_df.BraTS21ID, folder='test')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T20:50:00.502876Z","iopub.execute_input":"2021-08-16T20:50:00.503245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = format_data_for_keras(training_data_raw)\ny = train_df.MGMT_value[~train_df.BraTS21ID.isin(EXCLUDE)]\nX_test = format_data_for_keras(testing_data_raw)\n\nX.shape, X_test.shape, y.shape \n","metadata":{"execution":{"iopub.status.busy":"2021-08-13T19:37:37.644918Z","iopub.execute_input":"2021-08-13T19:37:37.645263Z","iopub.status.idle":"2021-08-13T19:37:37.761599Z","shell.execute_reply.started":"2021-08-13T19:37:37.645235Z","shell.execute_reply":"2021-08-13T19:37:37.760599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport pydicom\nimport ast\nimport cv2\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.model_selection import train_test_split\n\nX_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.25, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2021-08-13T19:37:43.022437Z","iopub.execute_input":"2021-08-13T19:37:43.022981Z","iopub.status.idle":"2021-08-13T19:37:43.050726Z","shell.execute_reply.started":"2021-08-13T19:37:43.022926Z","shell.execute_reply":"2021-08-13T19:37:43.049664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inpt = keras.Input(shape=X_train.shape[1:])\nh = keras.layers.experimental.preprocessing.Rescaling(1./255)(inpt)\n\n# convolutional layer!\nh = keras.layers.Conv2D(32, kernel_size=(4, 4), activation=\"relu\", name=\"Conv_1\")(h) \n# pooling layer\nh = keras.layers.MaxPool2D()(h) \n\n# convolutional layer!\nh = keras.layers.Conv2D(32, kernel_size=(4, 4), activation=\"relu\", name=\"Conv_2\")(h) \n# pooling layer\nh = keras.layers.MaxPool2D()(h)\n\nh = keras.layers.Flatten()(h)   \noutput = keras.layers.Dense(1, activation=\"sigmoid\")(h)\n\nmodel = keras.Model(inpt, output)","metadata":{"execution":{"iopub.status.busy":"2021-08-13T19:56:29.208686Z","iopub.execute_input":"2021-08-13T19:56:29.209077Z","iopub.status.idle":"2021-08-13T19:56:29.25975Z","shell.execute_reply.started":"2021-08-13T19:56:29.209042Z","shell.execute_reply":"2021-08-13T19:56:29.258789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.compile(loss='binary_crossentropy',\n             optimizer='adam',\n             metrics=['accuracy'])\nhistory = model.fit(x=X_train, y = y_train, \n                    epochs=20,\n                    validation_data= (X_valid, y_valid))\ny_pred = model.predict(X_valid)\ny_pred = np.reshape(y_pred, (y_pred.shape[0], ))\ny_pred.shape\nprint(roc_auc_score(y_valid, y_pred))","metadata":{"execution":{"iopub.status.busy":"2021-08-13T19:56:30.076661Z","iopub.execute_input":"2021-08-13T19:56:30.077032Z","iopub.status.idle":"2021-08-13T19:57:20.054815Z","shell.execute_reply.started":"2021-08-13T19:56:30.076999Z","shell.execute_reply":"2021-08-13T19:57:20.053758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = model.predict(X_test)\nprint(predictions)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv', converters={'BraTS21ID': lambda x: str(x)})\nsubmission['MGMT_value'] = predictions\n\nsubmission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T20:56:06.136792Z","iopub.execute_input":"2021-08-17T20:56:06.137200Z","iopub.status.idle":"2021-08-17T20:56:06.149057Z","shell.execute_reply.started":"2021-08-17T20:56:06.137165Z","shell.execute_reply":"2021-08-17T20:56:06.147791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}