{"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":"What follows is **a** solution for the 2021 [RSNA-MICCAI Brain Tumor Radiogenomic Classification](https://www.kaggle.com/c/rsna-miccai-brain-tumor-radiogenomic-classification/overview) Kaggle competition. This is definitely **not** a winning solution, but is offered as an example of an end-to-end solution. It uses TorchIO for data manipulation, NiBabel for reading NifiTi images, and Tensorflow to train a 3d Convolutional Neural Network.\n\nThe only **prerequisite** to running *this* kernel, is to \"Add data\" from the \"Notebook Output File\" of \"pip-download-torchio\" (https://www.kaggle.com/ohbewise/pip-download-torchio) to this kernel.  This is needed to enable pip to install TorchIO offline.\n\nNote with `demo = True` only 10 patients are processed, and when `demo = False` all patients are processed but the notebook runs out of ram.  For this purpose I have provided link to each step of the solution that can process all the patients and not run out of ram!","metadata":{}},{"cell_type":"markdown","source":"## Install and import libraries","metadata":{}},{"cell_type":"code","source":"# If this line fails please see the prerequisite above\n!pip install --quiet --no-index --find-links ../input/pip-download-torchio/ --requirement ../input/pip-download-torchio/requirements.txt","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:15:45.712371Z","iopub.execute_input":"2021-09-27T20:15:45.71299Z","iopub.status.idle":"2021-09-27T20:15:54.366102Z","shell.execute_reply.started":"2021-09-27T20:15:45.712889Z","shell.execute_reply":"2021-09-27T20:15:54.365246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import libraries\nimport os\nimport csv\nimport pickle\nimport numpy as np\nimport pandas as pd\nimport nibabel as nib\nimport torchio as tio\nimport tensorflow as tf\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\n\n# Parameters to limit the processing power needed.\ndemo  = True # if True limits to 10 patients\nscan_types    = ['FLAIR','T1w','T1wCE','T2w'] # uses all scan types","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:15:54.369718Z","iopub.execute_input":"2021-09-27T20:15:54.369953Z","iopub.status.idle":"2021-09-27T20:16:04.183012Z","shell.execute_reply.started":"2021-09-27T20:15:54.369926Z","shell.execute_reply":"2021-09-27T20:16:04.182271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocess data: DICOM to normalized NIfTI with TorchIO\n Uses TorchIO to convert folders of DICOM images into a NIfTI file. More importantly it normalizes, resizes and rotates the MRI scans. For a stand alone kernel of this section see https://www.kaggle.com/ohbewise/dicom-to-normalized-nifti-with-torchio","metadata":{}},{"cell_type":"code","source":"# Preprocess data \ndata_dir   = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification/'\nout_dir    = '/kaggle/working/processed'\n\nfor dataset in ['train']:\n    dataset_dir = f'{data_dir}{dataset}'\n    patients = os.listdir(dataset_dir)\n    if demo:\n        patients = patients[:10]\n    \n    # Remove cases the competion host said to exclude \n    # https://www.kaggle.com/c/rsna-miccai-brain-tumor-radiogenomic-classification/discussion/262046\n    if '00109' in patients: patients.remove('00109')\n    if '00123' in patients: patients.remove('00123')\n    if '00709' in patients: patients.remove('00709')\n    \n    print(f'Total patients in {dataset} dataset: {len(patients)}')\n\n    count = 0\n    for patient in patients:\n        count = count + 1\n        print(f'{dataset}: {count}/{len(patients)}')\n\n        for scan_type in scan_types:\n            scan_src  = f'{dataset_dir}/{patient}/{scan_type}/'\n            scan_dest = f'{out_dir}/{dataset}/{patient}/{scan_type}/'\n            Path(scan_dest).mkdir(parents=True, exist_ok=True)\n            image = tio.ScalarImage(scan_src)\n            transforms = [\n                tio.ToCanonical(),\n                tio.Resample(1),\n                tio.ZNormalization(masking_method=tio.ZNormalization.mean),\n                tio.RescaleIntensity((-1, 1)),\n                tio.CropOrPad((128,128,64)),\n            ]\n            transform = tio.Compose(transforms)\n            preprocessed = transform(image)\n            preprocessed.save(f'{scan_dest}/{scan_type}.nii.gz')","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:16:04.184329Z","iopub.execute_input":"2021-09-27T20:16:04.184588Z","iopub.status.idle":"2021-09-27T20:18:52.289285Z","shell.execute_reply.started":"2021-09-27T20:16:04.184557Z","shell.execute_reply":"2021-09-27T20:18:52.28839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Build datasets: NIfTI to Split Dataset with NiBabel\nUses NiBabel to read the NIfTI files in the \"processed/train\" folder and split them into a Training and Validation dataset. For a stand alone kernel of this section see https://www.kaggle.com/ohbewise/nifti-to-split-dataset-with-nibabel","metadata":{}},{"cell_type":"code","source":"# build datasets\n\n# dataset processing functions\ndef read_nifti_file(filepath):\n    \"\"\"Read and load volume\"\"\"\n    # Read file\n    scan = nib.load(filepath)\n    # Get raw data\n    scan = scan.get_fdata()\n    return scan\n\ndef add_batch_channel(volume):\n    \"\"\"Process validation data by adding a channel.\"\"\"\n    volume = tf.expand_dims(volume, axis=-1)\n    volume = tf.expand_dims(volume, axis=0)\n    return volume\n\ndef process_scan(filepath):\n    scan = read_nifti_file(filepath)\n    volume = add_batch_channel(scan)\n    return volume\n\n# get labels\nlabels_df = pd.read_csv(data_dir+'train_labels.csv', index_col=0)\n\n# split patients\npatients = os.listdir(f'{out_dir}/train')\nfrom sklearn.model_selection import train_test_split\ntrain, validation = train_test_split(patients, test_size=0.3, random_state=42)\nprint(f'{len(patients)} total patients.\\n   {len(train)} in the train split.\\n   {len(validation)} in the validation split')\n\nsplits_dict = {'train':train, 'validation':validation}\n\nfor scan_type in scan_types:\n    print(f'{scan_type} start')\n    for split_name, split_list in splits_dict.items():\n        print(f'   {split_name} start')\n        label_list = []\n        filepaths = []\n        for patient in split_list:\n            label = labels_df._get_value(int(patient), 'MGMT_value')\n            label = add_batch_channel(label)\n            label_list.append(label)\n            filepath  = f'{out_dir}/train/{patient}/{scan_type}/{scan_type}.nii.gz'\n            filepaths.append(filepath)\n\n        features = np.array([process_scan(filepath) for filepath in filepaths if filepath])\n        labels = np.array(label_list, dtype=np.uint8)\n        dataset = tf.data.Dataset.from_tensor_slices((features, labels))\n        \n        # save dataset   \n        tf_data_path = f'./datasets/{scan_type}_{split_name}_dataset'\n        tf.data.experimental.save(dataset, tf_data_path, compression='GZIP')\n        with open(tf_data_path + '/element_spec', 'wb') as out_:  # also save the element_spec to disk for future loading\n            pickle.dump(dataset.element_spec, out_)\n        print(f'   {split_name} done')\n    print(f'{scan_type} done')","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:18:52.291164Z","iopub.execute_input":"2021-09-27T20:18:52.291408Z","iopub.status.idle":"2021-09-27T20:19:13.479719Z","shell.execute_reply.started":"2021-09-27T20:18:52.291382Z","shell.execute_reply":"2021-09-27T20:19:13.478989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define, train, and evaluate model:  Dataset to Model with Tensorflow\n\nKernel uses Tensorflow to define, train, and evaluate a model.  For a stand alone kernel of this section see https://www.kaggle.com/ohbewise/dataset-to-model-with-tensorflow","metadata":{}},{"cell_type":"code","source":"# Define, train, and evaluate model\n# source: https://keras.io/examples/vision/3D_image_classification/\ndef get_model(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\n    # Define the model.\n    model = tf.keras.Model(inputs, outputs, name=name)\n    \n    # Compile model.\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    \n    return model","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:19:13.481001Z","iopub.execute_input":"2021-09-27T20:19:13.481351Z","iopub.status.idle":"2021-09-27T20:19:13.494693Z","shell.execute_reply.started":"2021-09-27T20:19:13.481314Z","shell.execute_reply":"2021-09-27T20:19:13.493925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for scan_type in scan_types:\n    # load train_dataset dataset\n    tf_data_path = f'./datasets/{scan_type}_train_dataset'\n    with open(tf_data_path + '/element_spec', 'rb') as in_:\n        es = pickle.load(in_)\n    train_dataset = tf.data.experimental.load(tf_data_path, es, compression='GZIP')\n    \n    # load validation_dataset\n    tf_data_path = f'./datasets/{scan_type}_validation_dataset'\n    with open(tf_data_path + '/element_spec', 'rb') as in_:\n        es = pickle.load(in_)\n    validation_dataset = tf.data.experimental.load(tf_data_path, es, compression='GZIP')\n\n    # Get Model\n    model = get_model(width=128, height=128, depth=64,name=scan_type)\n    \n    # Define callbacks.\n    checkpoint_cb = tf.keras.callbacks.ModelCheckpoint(\n        f'{scan_type}_3d_image_classification.h5', save_best_only=True\n    )\n    early_stopping_cb = tf.keras.callbacks.EarlyStopping(monitor=\"val_acc\", patience=15)\n\n    epochs = 100\n    model.fit(\n        train_dataset,\n        validation_data=validation_dataset,\n        epochs=epochs,\n        shuffle=True,\n        verbose=2,\n        callbacks=[checkpoint_cb, early_stopping_cb],\n    )\n    \n    #save model\n    model.save(f'./models/{scan_type}')\n    \n    # show metrics\n    fig, ax = plt.subplots(1, 2, figsize=(20, 3))\n    ax = ax.ravel()\n\n    for i, metric in enumerate([\"acc\", \"loss\"]):\n        ax[i].plot(model.history.history[metric])\n        ax[i].plot(model.history.history[\"val_\" + metric])\n        ax[i].set_title(\"{} Model {}\".format(scan_type, metric))\n        ax[i].set_xlabel(\"epochs\")\n        ax[i].set_ylabel(metric)\n        ax[i].legend([\"train\", \"val\"])","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:19:13.49591Z","iopub.execute_input":"2021-09-27T20:19:13.496615Z","iopub.status.idle":"2021-09-27T20:20:21.121999Z","shell.execute_reply.started":"2021-09-27T20:19:13.496579Z","shell.execute_reply":"2021-09-27T20:20:21.121237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Write predictions to submission.csv: Model Prediction to Submission\nKernel uses the model to predict the test set and write results to submission.csv.  For a stand alone kernel of this section see https://www.kaggle.com/ohbewise/model-prediction-to-submission","metadata":{}},{"cell_type":"code","source":"# write predictions to submission.csv\n\n# Set up directories\ndata_dir   = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification/'\ntest_dir   = f'{data_dir}test'\npatients = os.listdir(test_dir)\nif demo:\n    patients = patients[:10]\nprint(f'Total patients: {len(patients)}\\n\\n')\n\nout_dir    = '/kaggle/working/processed'\n\nscan_types = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\nscan_types = ['T1wCE']\n\nfor scan_type in scan_types:\n    f = open(f'/kaggle/working/submission.csv', 'w')\n    writer = csv.writer(f)\n    writer.writerow(['BraTS21ID','MGMT_value'])\n    for patient in patients:\n        # dicom to nifiti\n        scan_src  = f'{test_dir}/{patient}/{scan_type}/'\n        scan_dest = f'{out_dir}/test/{patient}/{scan_type}/'\n        Path(scan_dest).mkdir(parents=True, exist_ok=True)\n        image = tio.ScalarImage(scan_src)  # subclass of Image\n        transforms = [\n            tio.ToCanonical(),\n            tio.Resample(1),\n            tio.ZNormalization(masking_method=tio.ZNormalization.mean),\n            tio.RescaleIntensity((-1, 1)),\n            tio.CropOrPad((128,128,64)),\n        ]\n        transform = tio.Compose(transforms)\n        preprocessed = transform(image)\n        filepath = f'{scan_dest}/{scan_type}.nii.gz'\n        preprocessed.save(filepath)\n        \n        # process_scan\n        case = process_scan(filepath)\n\n        # tf model\n        model = tf.keras.models.load_model(f'./models/{scan_type}')\n\n        # get prediction\n        prediction = model.predict(case)\n        \n        # write prediction\n        print(f'{patient},{prediction[0][0]}')\n        writer.writerow([patient, prediction[0][0]])\n\n    f.close()","metadata":{"execution":{"iopub.status.busy":"2021-09-27T20:20:21.123348Z","iopub.execute_input":"2021-09-27T20:20:21.123901Z","iopub.status.idle":"2021-09-27T20:21:00.899674Z","shell.execute_reply.started":"2021-09-27T20:20:21.123863Z","shell.execute_reply":"2021-09-27T20:21:00.898107Z"},"trusted":true},"execution_count":null,"outputs":[]}]}