{"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":"import os, sys\nos.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport tensorflow as tf\nimport tensorflow_io as tfio\nimport pydicom\n\nfrom itertools import groupby\nprint(f\"Current Tensorflow version : {tf.__version__}\")\n\n#Turn TPU on\n# NEW on TPU in TensorFlow 24: shorter cross-compatible TPU/GPU/multi-GPU/cluster-GPU detection code\n\ntry: # detect TPUs\n    tpu = tf.distribute.cluster_resolver.TPUClusterResolver.connect() # TPU detection\n    strategy = tf.distribute.TPUStrategy(tpu)\n    # Step 1: Get the credential from the Cloud SDK\n    from kaggle_secrets import UserSecretsClient\n    user_secrets = UserSecretsClient()\n    user_credential = user_secrets.get_gcloud_credential()\n    \n    # Step 2: Set the credentials\n    user_secrets.set_tensorflow_credential(user_credential)\n\n    # Step 3: Use a familiar call to get the GCS path of the dataset\n    from kaggle_datasets import KaggleDatasets\n    GCS_DS_PATH = KaggleDatasets().get_gcs_path()\n    print(f\"TPU has been activated, the GCS dataset path : {GCS_DS_PATH}\", \"\\n\")\n    os.environ[\"TFHUB_MODEL_LOAD_FORMAT\"] = \"UNCOMPRESSED\" # tpu 환경에서는 이 옵션 켜야 tfhub 잘 됨\nexcept ValueError: # detect GPUs\n    tpu = False\n    strategy = tf.distribute.MirroredStrategy() # for GPU or multi-GPU machines\n    #strategy = tf.distribute.get_strategy() # default strategy that works on CPU and single GPU\n    #strategy = tf.distribute.experimental.MultiWorkerMirroredStrategy() # for clusters of multi-GPU machines\n\nimport tensorflow_addons as tfa\nfrom tensorflow import keras\nfrom tensorflow.keras import Input, Model\nfrom tensorflow.keras.models import load_model\nfrom tensorflow.keras.layers import Conv2D, Dense, Activation, BatchNormalization, LayerNormalization, MultiHeadAttention, Embedding, Subtract, Add, Multiply, GlobalAveragePooling2D, GlobalAveragePooling1D, LayerNormalization\nfrom tensorflow.keras.preprocessing.image import load_img, img_to_array\nimport albumentations as album\nimport cv2\nimport random\nfrom functools import partial\n\nfrom tqdm.notebook import tqdm\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nfrom sklearn.metrics import classification_report\nfrom sklearn.model_selection import train_test_split\nimport os\nimport matplotlib.pyplot as plt\n\n\n#Setting hyperparameters\n\nres = int(256*1.5) \nif tpu:\n    batch_size = 32\nelse:\n    batch_size = 16\nmaxlen = 128 \ntrain_N = int(40000/20)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-16T01:13:18.26262Z","iopub.execute_input":"2022-12-16T01:13:18.263745Z","iopub.status.idle":"2022-12-16T01:13:28.712726Z","shell.execute_reply.started":"2022-12-16T01:13:18.263641Z","shell.execute_reply":"2022-12-16T01:13:28.711736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> helper functions to read and visualize dcm files","metadata":{}},{"cell_type":"code","source":"dicom_dir = \"/kaggle/input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/stage_2_train\"\ndf_train = pd.read_csv(\"/kaggle/input/rsna-ich-detection-metadata/df_train_split.csv\")\ndf_val = pd.read_csv(\"/kaggle/input/rsna-ich-detection-metadata/df_val_splt.csv\")\nprint(df_train.head())\nexample_dicom = \"ID_01e7a9d10\"\ndicom_path = os.path.join(dicom_dir, example_dicom + \".dcm\")\n#Any, SAH, Intraventricular hem.\ndef dicom_to_tensor(dicom_path):\n    dataset = pydicom.dcmread(dicom_path)\n    tensor = np.array(dataset.pixel_array)\n    slope = dataset.RescaleSlope   # dicom header (Rescale slope)\n    intercept = dataset.RescaleIntercept   # dicom header (Rescale intercept)\n    center = dataset.WindowCenter   # dicom header (Window center)\n    width = dataset.WindowWidth   # dicom header (Window width)\n\n    if(type(dataset.WindowCenter) == pydicom.multival.MultiValue):\n            center = float(dataset.WindowCenter[0])\n            width = float(dataset.WindowWidth[0])       \n    else:    \n            center = float(dataset.WindowCenter)\n            width = float(dataset.WindowWidth)\n\n    tensor = slope*tensor + intercept\n    lbound, ubound = center - 0.5*width, center + 0.5*width\n    tensor[np.where(tensor < lbound)] = lbound\n    tensor[np.where(tensor > ubound)] = ubound\n    tensor = tf.image.resize(tensor[:,:,tf.newaxis], [res,res])\n    tensor = np.array(tensor)\n    try:\n        del dataset\n    except:\n        pass\n    return tensor\n\ntensor_example = dicom_to_tensor(dicom_path)\nplt.figure(figsize = (15,15))\nplt.title(\"Subarachnoid and Intraventricular hemorrhage\", fontsize = 20)\nplt.imshow(tensor_example, cmap = \"gray\")\nplt.colorbar()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-16T01:13:28.714287Z","iopub.execute_input":"2022-12-16T01:13:28.715564Z","iopub.status.idle":"2022-12-16T01:13:29.426418Z","shell.execute_reply.started":"2022-12-16T01:13:28.715523Z","shell.execute_reply":"2022-12-16T01:13:29.425523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class DicomDataGenerator(keras.utils.Sequence):\n    def __init__(self, dataframe, x_col, y_col, res, batch_size, shuffle = True):\n        self.df = dataframe\n        self.x_col = x_col\n        self.res = res\n        self.y_col = y_col\n        self.batch_size = batch_size\n        self.shuffle = shuffle\n        self.on_epoch_end()\n    \n    def dicom_to_tensor(self, dicom_path):\n        dataset = pydicom.dcmread(dicom_path)\n        tensor = np.array(dataset.pixel_array)\n        slope = dataset.RescaleSlope   # dicom header (Rescale slope)\n        intercept = dataset.RescaleIntercept   # dicom header (Rescale intercept)\n        center = dataset.WindowCenter   # dicom header (Window center)\n        width = dataset.WindowWidth   # dicom header (Window width)\n\n        if(type(dataset.WindowCenter) == pydicom.multival.MultiValue):\n                center = float(dataset.WindowCenter[0])\n                width = float(dataset.WindowWidth[0])       \n        else:    \n                center = float(dataset.WindowCenter)\n                width = float(dataset.WindowWidth)\n\n        tensor = slope*tensor + intercept\n        lbound, ubound = center - 0.5*width, center + 0.5*width\n        tensor[np.where(tensor < lbound)] = lbound\n        tensor[np.where(tensor > ubound)] = ubound\n        \n        tensor = tf.image.resize(tensor[:,:,tf.newaxis], [res,res])\n        tensor = np.array(tensor)\n        try:\n            del dataset\n        except:\n            pass\n        return tensor\n        \n    def on_epoch_end(self):\n        self.indexes = np.arange(len(self.df))\n        if self.shuffle:\n            np.random.shuffle(self.indexes)\n            \n    def __len__(self):\n        return int(np.floor(len(self.df) / self.batch_size))\n    \n    def __data_generation(self, img_name, y_list):\n        ## path를 받아 img화 및 token화 하여 실제로 Feeding할 데이터를 반환\n        X = []\n        y = []\n        for i, (fname, label) in enumerate(zip(img_name, y_list)):\n            dicom_path = os.path.join(dicom_dir, fname+\".dcm\")\n            img = self.dicom_to_tensor(dicom_path)\n            \n            X.append(img)\n            y.append(label)\n                \n        X = np.array(X)\n        y = np.array(y).astype(np.int32)\n                \n        return X, y\n        \n                \n    def __getitem__(self, index):\n        indexes = self.indexes[index * self.batch_size : (index + 1) * self.batch_size]\n        \n        img_name = [self.df[self.x_col].values[k] for k in indexes]\n        y_list = [self.df[self.y_col].values[k] for k in indexes]\n        \n        X, y = self.__data_generation(img_name, y_list)\n        \n        return X, y\n    \n\ndef get_train_gen():\n    return DicomDataGenerator(df_train, x_col = \"SOPInstanceUID\",\n                         res = res, y_col = \"any\", batch_size = batch_size)\n\ndef get_val_gen():\n    return DicomDataGenerator(df_val, x_col = \"SOPInstanceUID\",\n                         res = res, y_col = \"any\", batch_size = batch_size)\n\nval_ds = tf.data.Dataset.from_generator(get_val_gen, (tf.float32, tf.int32), output_shapes = ((None, res, res,1), (None, ) )).repeat()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-16T01:13:29.427908Z","iopub.execute_input":"2022-12-16T01:13:29.428545Z","iopub.status.idle":"2022-12-16T01:13:29.512503Z","shell.execute_reply.started":"2022-12-16T01:13:29.428508Z","shell.execute_reply":"2022-12-16T01:13:29.511304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clahe = album.Compose([ \n                     album.CLAHE(p = 1.0),\n                     \n                    ]\n                   )\n \ndef clahe_fn(image1):\n    if np.min(image1) < np.max(image1):\n        image1 = (image1 - np.min(image1)) / (np.max(image1) - np.min(image1))\n    else:\n        image1 = image1/255\n        image1 = tf.clip_by_value(image1, 0, 1)\n    image1 = np.array(image1*255.0).astype(np.uint8)\n    aug_data1 = clahe(image = image1)\n    aug_img1 = aug_data1['image']\n    aug_img1 = np.array(aug_img1).astype(np.uint8)\n    \n    return aug_img1\n\n\ndef preprocess_val_data(image, label):\n    aug_img = tf.numpy_function(func = clahe_fn, inp = [image], Tout = (tf.uint8))\n    return (aug_img, label)\n\nsup_val_ds = val_ds.unbatch().shuffle(buffer_size = 2560).map(preprocess_val_data, num_parallel_calls = tf.data.AUTOTUNE).batch(batch_size = batch_size).prefetch(tf.data.AUTOTUNE)\nvis = True\nif vis :\n    for img, lab in sup_val_ds.take(1):\n        sample_img_set = img\n        sample_labels = np.array(lab)\n    print(sample_labels)\n\n    fig, axes = plt.subplots(4,4, figsize = (20, 20))\n    axes = axes.flatten()\n    for idx, ax in enumerate(axes):\n        ax.imshow(sample_img_set[idx])\n        if sample_labels[idx] == 1:\n            ax.set_title(\"Hemorrhage\")\n        else:\n            ax.set_title(\"Normal/No hemorrhage\")\n        ax.axis(\"off\")\n    plt.suptitle(\"From Validation Dataset, after CLAHE\")\n    plt.show()\n    print(sample_labels)\n    print(sample_img_set.shape)\n    #del sample_img_set, sample_labels","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-16T01:13:29.514842Z","iopub.execute_input":"2022-12-16T01:13:29.515185Z","iopub.status.idle":"2022-12-16T01:14:49.439568Z","shell.execute_reply.started":"2022-12-16T01:13:29.515146Z","shell.execute_reply":"2022-12-16T01:14:49.438546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model import\nconv_sup = load_model(\"/kaggle/input/rsna-ich-detection-metadata/Fully_Supervised_models/Naive_Feature_Vector_Generator_effnet.h5\"\n                     )\nvit_sup = load_model(\"/kaggle/input/rsna-ich-detection-metadata/Fully_Supervised_models/Naive_Feature_Vector_Generator_vit.h5\")\nconv_sa_sup = load_model(\"/kaggle/input/rsna-ich-detection-metadata/Fully_Supervised_models/Naive_Feature_Vector_Generator_effnet_with_attention.h5\")\nswav_model = load_model(\"/kaggle/input/rsna-ich-detection-metadata/SwAV_TransferLearning.h5\")","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:14:49.440996Z","iopub.execute_input":"2022-12-16T01:14:49.441697Z","iopub.status.idle":"2022-12-16T01:15:18.102097Z","shell.execute_reply.started":"2022-12-16T01:14:49.441658Z","shell.execute_reply":"2022-12-16T01:15:18.100793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model evaluation with crosstables","metadata":{}},{"cell_type":"code","source":"conv_sup.summary(), vit_sup.summary(), conv_sa_sup.summary()","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:15:18.104667Z","iopub.execute_input":"2022-12-16T01:15:18.105464Z","iopub.status.idle":"2022-12-16T01:15:18.157516Z","shell.execute_reply.started":"2022-12-16T01:15:18.10542Z","shell.execute_reply":"2022-12-16T01:15:18.156278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"swav_model.summary()","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:15:18.159072Z","iopub.execute_input":"2022-12-16T01:15:18.159672Z","iopub.status.idle":"2022-12-16T01:15:18.187565Z","shell.execute_reply.started":"2022-12-16T01:15:18.159625Z","shell.execute_reply":"2022-12-16T01:15:18.186025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"eval_numeric = True\nif eval_numeric : \n    for model in [conv_sup, vit_sup, conv_sa_sup, swav_model]:\n        print(f\"Evaluation of model name : {model.name}\", \"\\n\")\n        tp = tf.keras.metrics.TruePositives()\n        tn = tf.keras.metrics.TrueNegatives()\n        fp = tf.keras.metrics.FalsePositives()\n        fn = tf.keras.metrics.FalseNegatives()\n        model.compile(loss = model.loss,\n                     metrics =  [\"accuracy\", \"AUC\", tfa.metrics.F1Score(num_classes = 1), tp, tn, fp, fn],\n                     optimizer = model.optimizer )\n        model.evaluate(val_ds, steps = len(df_val)//batch_size)\n        print(\"---------------------\", \"\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:15:18.189241Z","iopub.execute_input":"2022-12-16T01:15:18.18971Z","iopub.status.idle":"2022-12-16T01:15:18.198056Z","shell.execute_reply.started":"2022-12-16T01:15:18.189665Z","shell.execute_reply":"2022-12-16T01:15:18.196745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model evaluation with CAMs","metadata":{}},{"cell_type":"code","source":"def get_CAM(image, model, layer, \n           cam_type = \"gradcam\", flow_type = \"input\", pred_idx = None):\n    if flow_type == \"input\":\n        conv_layer = model.get_layer(layer).input\n    else:\n        conv_layer = model.get_layer(layer).output\n    \n    if len(image.shape) == 3:\n        image = image[tf.newaxis, :, :, :]\n        \n    gradient_model = Model(inputs=model.inputs,\n                          outputs = [conv_layer, model.outputs])\n    \n    with tf.GradientTape() as tape:\n        conv_tensor, total_proba = gradient_model(image)\n        total_proba = total_proba[0]\n        #print(total_proba)\n        if pred_idx is None:\n            pred_idx = np.argmax(total_proba)\n            #print(pred_idx)\n        proba = total_proba[pred_idx]\n    grads = tape.gradient(proba, conv_tensor)\n    if len(grads.shape) == 3:\n        _, prod, ndim = conv_layer.shape\n        w = tf.math.sqrt(float(prod))\n        h = tf.math.sqrt(float(prod))\n    else:\n        _, w, h, ndim = grads.shape\n    grads = tf.reshape(grads, [-1, w, h, ndim])\n    if cam_type == \"gradcam\":\n        grads = keras.layers.GlobalAveragePooling2D()(grads) #batch, ndim\n        grads = tf.reshape(grads, [-1,1])\n        conv_tensor = tf.reshape(conv_tensor, [1, w, h, ndim])\n        gcam = conv_tensor @ grads\n        gcam = tf.math.maximum(gcam, 0)\n        gcam = tf.squeeze(gcam)\n        return gcam\n    elif cam_type == \"layercam\":\n        grads = tf.math.maximum(grads, 0)\n        conv_tensor = tf.reshape(conv_tensor, [-1, w, h, ndim])\n        a_hat = tf.math.multiply(grads, conv_tensor)\n        a_hat = tf.squeeze(a_hat)\n        layercam = tf.reduce_mean(a_hat, axis = -1)\n        return layercam\n    else:\n        print(\"CAM type should be gradcam or layercam\")\n        return grads","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:37:18.249732Z","iopub.execute_input":"2022-12-16T01:37:18.250142Z","iopub.status.idle":"2022-12-16T01:37:18.2661Z","shell.execute_reply.started":"2022-12-16T01:37:18.250108Z","shell.execute_reply":"2022-12-16T01:37:18.26475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Superimpose\ndef get_cam_viz(image, model, feature_map, alpha = 0.2, ground_truth = None):\n    if len(image.shape) == 3:\n        image = image[tf.newaxis, :, :, :]\n    \n    feature_map = np.array(feature_map)\n    if (np.max(feature_map) - np.min(feature_map)) == 0:\n        feature_map /= 255\n    else:\n        feature_map = (feature_map-np.min(feature_map)) / (np.max(feature_map) - np.min(feature_map))\n    feature_map = cv2.resize(feature_map, [res,res])\n    x = 255*feature_map\n    x = np.array(x).astype(np.uint8)\n    heatmap = cv2.applyColorMap(x, cv2.COLORMAP_INFERNO)\n    heatmap = heatmap/255\n    superimposed = alpha*heatmap + (1-alpha)*np.array(image[0]).astype(np.float32)/255\n    #superimposed = (superimposed - np.min(superimposed)) / (np.max(superimposed) - np.min(superimposed))\n    \n    pred_proba = model.predict(image)\n    pred_idx = np.round(float(pred_proba[0]))\n    \n    fig, axes = plt.subplots(1, 3, figsize=(45,15))\n    axes = axes.flatten()\n    axes[0].imshow(tf.squeeze(image), cmap = 'gray')\n    if ground_truth:\n        axes[0].set_title(ground_truth, fontsize = 40)\n    axes[1].imshow(heatmap)\n    axes[1].set_title(\"Heatmap of given CT\", fontsize = 40)\n    axes[2].imshow(superimposed)\n    axes[2].set_title(f\"Prediction : {pred_idx}, probability : {np.round(pred_proba[0],2)}\", fontsize = 30)\n    plt.suptitle(f\"{model.name}_Evaluation\", fontsize = 35)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:39:21.214092Z","iopub.execute_input":"2022-12-16T01:39:21.214541Z","iopub.status.idle":"2022-12-16T01:39:22.010369Z","shell.execute_reply.started":"2022-12-16T01:39:21.214503Z","shell.execute_reply":"2022-12-16T01:39:22.008584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for image, label in zip(sample_img_set, sample_labels):\n    truth_title = str(label)+\"_GradCAM\"\n    gradcam_conv = get_CAM(image = image, model = swav_model, layer = \"GAP_layer_type_effnet\",\n                       cam_type = \"gradcam\")\n    get_cam_viz(image, swav_model, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"---------------------------------------------------------------------\", \"\\n\")\n    truth_title = str(label)+\"_LayerCAM\"\n    layercam_conv = get_CAM(image = image, model = swav_model, layer = \"GAP_layer_type_effnet\",\n                           cam_type = \"layercam\")\n    get_cam_viz(image, swav_model, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"=====================================================================\", \"\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:37:23.132364Z","iopub.execute_input":"2022-12-16T01:37:23.132775Z","iopub.status.idle":"2022-12-16T01:38:51.538631Z","shell.execute_reply.started":"2022-12-16T01:37:23.132743Z","shell.execute_reply":"2022-12-16T01:38:51.537616Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for image, label in zip(sample_img_set, sample_labels):\n    truth_title = str(label)+\"_GradCAM\"\n    gradcam_conv = get_CAM(image = image, model = conv_sup, layer = \"GAP_layer_type_effnet\",\n                       cam_type = \"gradcam\")\n    get_cam_viz(image, conv_sup, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"---------------------------------------------------------------------\", \"\\n\")\n    truth_title = str(label)+\"_LayerCAM\"\n    layercam_conv = get_CAM(image = image, model = conv_sup, layer = \"GAP_layer_type_effnet\",\n                           cam_type = \"layercam\")\n    get_cam_viz(image, conv_sup, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"=====================================================================\", \"\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:39:23.693419Z","iopub.execute_input":"2022-12-16T01:39:23.693823Z","iopub.status.idle":"2022-12-16T01:40:56.117536Z","shell.execute_reply.started":"2022-12-16T01:39:23.693792Z","shell.execute_reply":"2022-12-16T01:40:56.116727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for image, label in zip(sample_img_set, sample_labels):\n    truth_title = str(label)+\"_GradCAM\"\n    gradcam_conv = get_CAM(image = image, model = conv_sa_sup, layer = \"GAP_layer_type_effnet_with_attention\",\n                       cam_type = \"gradcam\")\n    get_cam_viz(image, conv_sup, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"---------------------------------------------------------------------\", \"\\n\")\n    truth_title = str(label)+\"_LayerCAM\"\n    layercam_conv = get_CAM(image = image, model = conv_sa_sup, layer = \"GAP_layer_type_effnet_with_attention\",\n                           cam_type = \"layercam\")\n    get_cam_viz(image, conv_sup, gradcam_conv, alpha = 0.2, ground_truth = truth_title)\n    print(\"=====================================================================\", \"\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-12-16T01:40:56.11918Z","iopub.execute_input":"2022-12-16T01:40:56.120176Z"},"trusted":true},"execution_count":null,"outputs":[]}]}