{"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":"code","source":"!pip install ../input/nifti-converter/dicom2nifti-2.3.0/dicom2nifti-2.3.0\n#file:///srv/pkg/mypackage\n#!pip install dicom2nifti","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:20.873858Z","iopub.execute_input":"2021-12-23T06:25:20.874183Z","iopub.status.idle":"2021-12-23T06:25:52.277918Z","shell.execute_reply.started":"2021-12-23T06:25:20.874106Z","shell.execute_reply":"2021-12-23T06:25:52.277071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport glob\n\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nimport random\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\n\nimport pydicom # Handle MRI images\nimport dicom2nifti\nfrom dicom2nifti.exceptions import ConversionValidationError\n\nimport cv2  # OpenCV - https://docs.opencv.org/master/d6/d00/tutorial_py_root.html\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_auc_score\n\nfrom scipy import ndimage\n\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras.utils import to_categorical\nfrom tensorflow.keras import layers","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:52.279932Z","iopub.execute_input":"2021-12-23T06:25:52.280195Z","iopub.status.idle":"2021-12-23T06:25:58.467168Z","shell.execute_reply.started":"2021-12-23T06:25:52.280161Z","shell.execute_reply":"2021-12-23T06:25:58.466411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_dir = Path('../input/rsna-miccai-brain-tumor-radiogenomic-classification/')\n\nmri_types = [\"FLAIR\", \"T1w\", \"T2w\", \"T1wCE\"]\nexcluded_images = [109, 123, 709] # Bad images","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.468613Z","iopub.execute_input":"2021-12-23T06:25:58.46887Z","iopub.status.idle":"2021-12-23T06:25:58.47678Z","shell.execute_reply.started":"2021-12-23T06:25:58.468834Z","shell.execute_reply":"2021-12-23T06:25:58.475109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(data_dir / \"train_labels.csv\")\ntest_df = pd.read_csv(data_dir / \"sample_submission.csv\")\nsample_submission = pd.read_csv(data_dir / \"sample_submission.csv\")\n\ntrain_df = train_df[~train_df.BraTS21ID.isin(excluded_images)]\n\nprint(f\"train data: Rows={train_df.shape[0]}, Columns={train_df.shape[1]}\")\nprint(f\"test data : Rows={test_df.shape[0]}, Columns={test_df.shape[1]}\")","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.48065Z","iopub.execute_input":"2021-12-23T06:25:58.481383Z","iopub.status.idle":"2021-12-23T06:25:58.52504Z","shell.execute_reply.started":"2021-12-23T06:25:58.481312Z","shell.execute_reply":"2021-12-23T06:25:58.524299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resize_volume(img):\n    \"\"\"Resize across z-axis\"\"\"\n    # Set the desired depth\n    desired_depth = 64\n    desired_width = 128\n    desired_height = 128\n    # Get current depth\n    current_depth = img.shape[-1]\n    current_width = img.shape[0]\n    current_height = img.shape[1]\n    # Compute depth factor\n    depth = current_depth / desired_depth\n    width = current_width / desired_width\n    height = current_height / desired_height\n    depth_factor = 1 / depth\n    width_factor = 1 / width\n    height_factor = 1 / height\n    # Resize across z-axis\n    img = ndimage.zoom(img, (width_factor, height_factor, depth_factor), order=1)\n    return img","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.526333Z","iopub.execute_input":"2021-12-23T06:25:58.526642Z","iopub.status.idle":"2021-12-23T06:25:58.53429Z","shell.execute_reply.started":"2021-12-23T06:25:58.526603Z","shell.execute_reply":"2021-12-23T06:25:58.532582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom(path, size = 224):\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))","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.536281Z","iopub.execute_input":"2021-12-23T06:25:58.536641Z","iopub.status.idle":"2021-12-23T06:25:58.545405Z","shell.execute_reply.started":"2021-12-23T06:25:58.536596Z","shell.execute_reply":"2021-12-23T06:25:58.54448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom2(path):\n    data = np.concatenate([tf.expand_dims(pydicom.read_file(p).pixel_array, axis=-1) for p in path], axis=2)\n    if np.max(data) != 0:\n        data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return resize_volume(data)","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.547008Z","iopub.execute_input":"2021-12-23T06:25:58.547842Z","iopub.status.idle":"2021-12-23T06:25:58.55566Z","shell.execute_reply.started":"2021-12-23T06:25:58.547797Z","shell.execute_reply":"2021-12-23T06:25:58.554717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_all_image_paths(brats21id, image_type, folder='train'):\n    assert(image_type in mri_types)\n    \n    patient_path = os.path.join(\n        \"../input/rsna-miccai-brain-tumor-radiogenomic-classification\",\n        folder, \n        str(brats21id).zfill(5),\n    )\n\n    paths = sorted(\n        glob.glob(os.path.join(patient_path, image_type, \"*\")), \n        key=lambda x: int(x[:-4].split(\"-\")[-1]),\n    )\n    \n    num_images = len(paths)\n    \n    start = int(num_images * 0.25)\n    end = int(num_images * 0.75)\n\n    interval = 3\n    \n    if num_images < 10: \n        interval = 1\n    \n    return np.array(paths[start:end:interval])\n\ndef get_all_images2(brats21id, image_type, folder='train'):\n    return [load_dicom2(get_all_image_paths(brats21id, image_type, folder))]","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.557272Z","iopub.execute_input":"2021-12-23T06:25:58.557736Z","iopub.status.idle":"2021-12-23T06:25:58.568127Z","shell.execute_reply.started":"2021-12-23T06:25:58.557695Z","shell.execute_reply":"2021-12-23T06:25:58.567254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def first_last_true(arr):\n    first = 0\n    for i,ele in enumerate(arr):\n        if ele:\n            first = i\n            break\n    last = len(arr)\n    while i > 1:\n        if arr[last-1]:\n            break\n        last -= 1\n    return first, last\n\ndef remove_zeros(data):\n    axis0 = data.any(axis=(1,2))\n    axis0s, axis0e = first_last_true(axis0)\n    axis1 = data.any(axis=(0,2))\n    axis1s, axis1e = first_last_true(axis1)\n    axis2 = data.any(axis=(0,1))\n    axis2s, axis2e = first_last_true(axis2)\n    return data[axis0s:axis0e, axis1s:axis1e, axis2s:axis2e]\n\ndef load_nii_from_dicom_series(path):\n    data = dicom2nifti.dicom_series_to_nifti(path, output_file=\"test.nii\")\n    data = remove_zeros(data[\"NII\"].get_fdata())\n    data -= np.min(data)\n    max_val = np.max(data)\n    if max_val > 0:\n        data /= max_val\n    \n    num_images = data.shape[0]\n    \n    start = int(num_images * 0.25)\n    end = int(num_images * 0.75)\n\n    interval = 3\n    \n    if num_images < 10: \n        interval = 1\n        \n    return resize_volume((data[start:end:interval] * 255).astype(np.uint8))\n\ndef get_all_images3(brats21id, image_type, folder):\n    assert(image_type in mri_types)\n    \n    patient_path = os.path.join(\n        \"../input/rsna-miccai-brain-tumor-radiogenomic-classification\",\n        folder, str(brats21id).zfill(5), image_type\n    )\n    return [load_nii_from_dicom_series(patient_path)]","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.569757Z","iopub.execute_input":"2021-12-23T06:25:58.570483Z","iopub.status.idle":"2021-12-23T06:25:58.584952Z","shell.execute_reply.started":"2021-12-23T06:25:58.570439Z","shell.execute_reply":"2021-12-23T06:25:58.583965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_all_data_for_train(image_type):\n    global train_df\n    \n    X = []\n    y = []\n    train_ids = []\n\n    for i in tqdm(train_df.index):\n        x = train_df.loc[i]\n        try:\n            images = get_all_images3(int(x['BraTS21ID']), image_type, 'train')\n        except ConversionValidationError:\n            print(\"SLICE_INCREMENT_INCONSISTENT for\",x['BraTS21ID'])\n            images = get_all_images2(int(x['BraTS21ID']), image_type, 'train')\n        label = x['MGMT_value']\n\n        X += images\n        y += [label] * len(images)\n        train_ids += [int(x['BraTS21ID'])] * len(images)\n        assert(len(X) == len(y))\n    return np.array(X), np.array(y), np.array(train_ids)","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.588878Z","iopub.execute_input":"2021-12-23T06:25:58.589165Z","iopub.status.idle":"2021-12-23T06:25:58.611061Z","shell.execute_reply.started":"2021-12-23T06:25:58.58913Z","shell.execute_reply":"2021-12-23T06:25:58.609687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_all_data_for_test(image_type):\n    global test_df\n    \n    X = []\n    test_ids = []\n\n    for i in tqdm(test_df.index):\n        x = test_df.loc[i]\n        try:\n            images = get_all_images3(int(x['BraTS21ID']), image_type, 'test')\n        except ConversionValidationError:\n            print(\"SLICE_INCREMENT_INCONSISTENT for\",x['BraTS21ID'])\n            images = get_all_images2(int(x['BraTS21ID']), image_type, 'test')\n        X += images\n        test_ids += [int(x['BraTS21ID'])] * len(images)\n\n    return np.array(X), np.array(test_ids)","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.612972Z","iopub.execute_input":"2021-12-23T06:25:58.613518Z","iopub.status.idle":"2021-12-23T06:25:58.623635Z","shell.execute_reply.started":"2021-12-23T06:25:58.613474Z","shell.execute_reply":"2021-12-23T06:25:58.622683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X, y, trainidt = get_all_data_for_train('T2w')\nX_test, testidt = get_all_data_for_test('T2w')","metadata":{"execution":{"iopub.status.busy":"2021-12-23T06:25:58.625269Z","iopub.execute_input":"2021-12-23T06:25:58.625807Z","iopub.status.idle":"2021-12-23T07:02:56.507505Z","shell.execute_reply.started":"2021-12-23T06:25:58.625767Z","shell.execute_reply":"2021-12-23T07:02:56.506753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# source: https://keras.io/examples/vision/3D_image_classification/\ndef get_3DCNNmodel(width=128, height=128, depth=64, name='3dcnn'):\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n\n    inputs = tf.keras.Input((width, height, depth, 1))\n\n    x = tf.keras.layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(inputs)\n    x = tf.keras.layers.MaxPool3D(pool_size=2)(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n\n    x = tf.keras.layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = tf.keras.layers.MaxPool3D(pool_size=2)(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n\n    x = tf.keras.layers.Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n    x = tf.keras.layers.MaxPool3D(pool_size=2)(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n\n    x = tf.keras.layers.Conv3D(filters=256, kernel_size=3, activation=\"relu\")(x)\n    x = tf.keras.layers.MaxPool3D(pool_size=2)(x)\n    x = tf.keras.layers.BatchNormalization()(x)\n\n    x = tf.keras.layers.GlobalAveragePooling3D()(x)\n    x = tf.keras.layers.Dense(units=512, activation=\"relu\")(x)\n    x = tf.keras.layers.Dropout(0.3)(x)\n\n    #outputs = tf.keras.layers.Dense(units=1, activation=\"sigmoid\")(x)\n    output = keras.layers.Dense(2, activation=\"softmax\")(x)\n\n    #model = tf.keras.Model(inputs, outputs, name=name)\n\n    #initial_learning_rate = 0.0001\n    #lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(\n    #    initial_learning_rate, decay_steps=100000, decay_rate=0.96, staircase=True\n    #)\n    #model.compile(\n    #    loss=\"binary_crossentropy\",\n    #    optimizer=tf.keras.optimizers.Adam(learning_rate=lr_schedule),\n    #    metrics=[\"acc\"],\n    #)\n    model = keras.Model(inputs, output)\n    initial_learning_rate =  0.0001\n    lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(\n        initial_learning_rate,\n        decay_steps=100000,\n        decay_rate=0.96, \n        staircase=True\n    )\n  \n    roc_auc = tf.keras.metrics.AUC(name='roc_auc', curve='ROC')\n\n    model.compile(\n        loss=\"categorical_crossentropy\", \n        optimizer=keras.optimizers.Adam(),\n        metrics=[roc_auc],\n    )\n    return model","metadata":{"execution":{"iopub.status.busy":"2021-12-23T07:02:56.50876Z","iopub.execute_input":"2021-12-23T07:02:56.509227Z","iopub.status.idle":"2021-12-23T07:02:56.524109Z","shell.execute_reply.started":"2021-12-23T07:02:56.509189Z","shell.execute_reply":"2021-12-23T07:02:56.523301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"early_stopping_cb = tf.keras.callbacks.EarlyStopping(monitor=\"val_roc_auc\", mode='max', patience=10) #patience=10","metadata":{"execution":{"iopub.status.busy":"2021-12-23T07:18:59.967916Z","iopub.execute_input":"2021-12-23T07:18:59.96843Z","iopub.status.idle":"2021-12-23T07:18:59.972812Z","shell.execute_reply.started":"2021-12-23T07:18:59.968394Z","shell.execute_reply":"2021-12-23T07:18:59.972032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"auc_list = []\nbest_auc = float(\"-inf\")\n\nfor i in tqdm(range(20)):\n    checkpoint_filepath = \"best_model_\"+str(i)+\".h5\"\n\n    model_checkpoint_cb = tf.keras.callbacks.ModelCheckpoint(\n        filepath=checkpoint_filepath,\n        save_weights_only=False,\n        monitor=\"val_roc_auc\",\n        mode=\"max\",\n        save_best_only=True,\n        save_freq=\"epoch\",\n        verbose=0,\n    )\n    X_train, X_valid, y_train, y_valid, trainidt_train, trainidt_valid = train_test_split(X, y, trainidt, test_size=0.2, random_state=i)\n\n    X_train = tf.expand_dims(X_train, axis=-1)\n    X_valid = tf.expand_dims(X_valid, axis=-1)\n    y_train = to_categorical(y_train)\n    y_valid = to_categorical(y_valid)\n\n    model = get_3DCNNmodel()\n\n    history = model.fit(x=X_train, y = y_train, epochs=40, batch_size = 2,\n                        callbacks=[model_checkpoint_cb, early_stopping_cb],\n                        validation_data=(X_valid, y_valid), verbose=0)\n\n    model_best = tf.keras.models.load_model(filepath=checkpoint_filepath)\n    y_pred = model_best.predict(X_valid,batch_size = 2)\n\n    pred = np.argmax(y_pred, axis=1)\n\n    result = pd.DataFrame(trainidt_valid)\n    result[1] = pred\n\n    result.columns = [\"BraTS21ID\", \"MGMT_value\"]\n    result2 = result.groupby(\"BraTS21ID\", as_index=False).mean()\n\n    result2 = result2.merge(train_df, on=\"BraTS21ID\")\n    auc = roc_auc_score(\n        result2.MGMT_value_y,\n        result2.MGMT_value_x,\n    )\n    print(f\"Validation AUC={auc}\")\n    auc_list.append(auc)\n    if auc > best_auc:\n        best_i = i\n        best_auc = auc","metadata":{"execution":{"iopub.status.busy":"2021-12-23T07:19:00.072128Z","iopub.execute_input":"2021-12-23T07:19:00.072411Z","iopub.status.idle":"2021-12-23T09:46:49.496436Z","shell.execute_reply.started":"2021-12-23T07:19:00.072355Z","shell.execute_reply":"2021-12-23T09:46:49.495704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(auc_list)\nplt.xlabel(\"AUC\")\nplt.ylabel(\"No. of trials\")\nplt.title(f\"Mean AUC = {np.mean(auc_list)}\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-12-23T09:46:49.498499Z","iopub.execute_input":"2021-12-23T09:46:49.499001Z","iopub.status.idle":"2021-12-23T09:46:49.706043Z","shell.execute_reply.started":"2021-12-23T09:46:49.498943Z","shell.execute_reply":"2021-12-23T09:46:49.705395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"checkpoint_filepath = \"best_model_\"+str(best_i)+\".h5\"\nprint(f\"Using {checkpoint_filepath} with AUC = {best_auc}.\")\nmodel_best = tf.keras.models.load_model(filepath=checkpoint_filepath)\ny_pred = model_best.predict(X_test,batch_size = 2)\n\npred = np.argmax(y_pred, axis=1) #\n\nresult = pd.DataFrame(testidt)\nresult[1] = pred\npred","metadata":{"execution":{"iopub.status.busy":"2021-12-23T09:46:49.707627Z","iopub.execute_input":"2021-12-23T09:46:49.707877Z","iopub.status.idle":"2021-12-23T09:46:51.229195Z","shell.execute_reply.started":"2021-12-23T09:46:49.707842Z","shell.execute_reply":"2021-12-23T09:46:51.228537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result.columns=['BraTS21ID','MGMT_value']\n\nresult2 = result.groupby('BraTS21ID',as_index=False).mean()\nresult2['BraTS21ID'] = sample_submission['BraTS21ID']\n\nresult2['MGMT_value'] = result2['MGMT_value'].apply(lambda x:round(x*10)/10)\n\nresult2.to_csv('submission.csv',index=False)\nresult2","metadata":{"execution":{"iopub.status.busy":"2021-12-23T09:46:51.230542Z","iopub.execute_input":"2021-12-23T09:46:51.230804Z","iopub.status.idle":"2021-12-23T09:46:51.250339Z","shell.execute_reply.started":"2021-12-23T09:46:51.230768Z","shell.execute_reply":"2021-12-23T09:46:51.249493Z"},"trusted":true},"execution_count":null,"outputs":[]}]}